home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_09_05 / 9n05086a < prev    next >
Text File  |  1991-03-21  |  11KB  |  253 lines
  1.  
  2. /**********************************************************
  3.  *
  4.  * Copyright (c) 1991, John Forkosh.  All rights reserved.
  5.  * Released to the public domain only for use that is both
  6.  * (a) by an individual, and (b) not for profit.
  7.  * --------------------------------------------------------
  8.  *
  9.  * Function:    simq ( a, b, n )
  10.  *
  11.  * Purpose:     Solves the set of n simultaneous linear
  12.  *              equations ax=b.
  13.  *
  14.  * Arguments:
  15.  *              a (I)   Address of n-by-n array of doubles
  16.  *                      stored columnwise (see notes)
  17.  *                      containing the matrix of
  18.  *                      coefficients (which are destroyed
  19.  *                      during the computation).
  20.  *              b (I/O) Address of vector of n doubles,
  21.  *                      containing the original constants,
  22.  *                      and returning the final solution.
  23.  *              n (I)   Int containing the number of
  24.  *                      equations to be solved (must be >1)
  25.  *
  26.  * Returns:     (int)   0 for a normal solution, or
  27.  *                      1 for a singular set of equations.
  28.  *
  29.  * Notes:     o Simq() is derived from the Fortran routine
  30.  *              of the same name in IBM's System/360
  31.  *              Scientific Subroutine Package, Publication
  32.  *              GH20-0205-4, Fifth Edition (August 1970),
  33.  *              Page 120.  See the discussion of Gauss-
  34.  *              Jordan elimination in Chapter 2 of 
  35.  *              Numerical Recipes in C for a similar
  36.  *              treatment.
  37.  *            o The matrix of coefficients, a[], must be
  38.  *              stored columnwise in a singly-dimensioned
  39.  *              array, e.g., for a 3x3 matrix the nine
  40.  *              elements of a[] are
  41.  *                      a[0]  a[3]  a[6]
  42.  *                      a[1]  a[4]  a[7]
  43.  *                      a[2]  a[5]  a[8]
  44.  *              In other words, the index for row i and
  45.  *              column j (i,j=0,...,n-1) is a[n*j+i].
  46.  *              Note that many index calculations are
  47.  *              buried in initialization steps of loops.
  48.  *            o The method of solution is by elimination
  49.  *              using largest pivotal divisor.  Each stage
  50.  *              of elimination consists of interchanging
  51.  *              rows when necessary to avoid division by
  52.  *              zero or small elements.  First the forward
  53.  *              solution to obtain the n-th variable is
  54.  *              done in n stages.  Then the back solution
  55.  *              for the other variables is done by
  56.  *              successive substitution.  Final solution
  57.  *              values are returned in vector b[], with
  58.  *              variable 1 in b[0], variable 2 in b[1],
  59.  *              etc.  If no pivot can be found exceeding
  60.  *              the #defined TOLerance, the matrix is
  61.  *              considered singular and an error returned.
  62.  *
  63.  * Source:      SIMQ.C
  64.  *
  65.  * --------------------------------------------------------
  66.  * Revision History:
  67.  *
  68.  * 01/07/91     J.Forkosh       Installation
  69.  *
  70.  *********************************************************/
  71.  
  72. /* --- standard headers --- */
  73. #include <stdio.h>
  74. #define dabs(x) ((x)>=0.0?(x):(-(x))) /* absolute value */
  75.  
  76. /* --- user-adjustable tolerance --- */
  77. #define TOL 0.0
  78.  
  79. /* --- set testdrive to compile (or not) test main() --- */
  80. #define TESTDRIVE 0             /* 1=compile it (0=not) */
  81.  
  82. /* --------------------------------------------------------
  83. Entry point
  84. -------------------------------------------------------- */
  85. int     simq ( a, b, n )        /* returns 0=ok, 1=error */
  86. double  *a;                     /* coefficient matrix */
  87. double  *b;                     /* constant vector */
  88. int     n;                      /* number of equations */
  89. {
  90. /* --- local allocations and declarations --- */
  91. int     i,j, ix,jx;             /* row,col indexes */
  92. int     it,k;                   /* temp indexes */
  93.  
  94. /* --------------------------------------------------------
  95. Forward Solution
  96. -------------------------------------------------------- */
  97. for ( j=0; j<n; j++ )           /* cols */
  98.   {
  99.   /* --- local allocations and declarations --- */
  100.   int   imax;                   /* pivot row */
  101.   double biga=0.0;              /*largest element in col*/
  102.   double save;                  /*save value during pivot*/
  103.  
  104.   /* ------------------------------------------------------
  105.   Find maximum coefficient (pivot row) in column
  106.   ------------------------------------------------------ */
  107.   for ( it=j*n,i=j; i<n; i++ )  /*rows in lower diagonal*/
  108.     if ( dabs(a[it+i]) > dabs(biga) ) /* got bigger biga */
  109.       { biga = a[it+i];         /* so store biggest */
  110.         imax = i; }             /* and its col offset */
  111.   /* --- error if pivot less than tolerance --- */
  112.   if ( dabs(biga) <= TOL )      /* pivot too small so... */
  113.     return ( 1 );               /*back to caller with err*/
  114.  
  115.   /* ------------------------------------------------------
  116.   Interchange rows as necessary and divide by leading coeff
  117.   ------------------------------------------------------ */
  118.   /* --- easy to interchange constant vector --- */
  119.   save = b[imax];               /* save b[imax] */
  120.   b[imax] = b[j];               /* replace it with b[j] */
  121.   b[j] = save/biga;             /* swap and rescale */
  122.   /* --- must interchange row one col at a time --- */
  123.   for ( i=j*(n+1),it=imax-j,k=j; k<n; k++,i+=n )
  124.     {
  125.     save = a[it+i];             /* save swap value */
  126.     a[it+i] = a[i];             /* replace it with a[i] */
  127.     a[i] = save/biga;           /* swap and rescale */
  128.     } /* --- end-of-for(k=j...n-1) --- */
  129.  
  130.   /* ------------------------------------------------------
  131.   Eliminate next variable (Note: the index calculations
  132.   required beyond this point get much more complicated.)
  133.   ------------------------------------------------------ */
  134.   if ( j < (n-1) )              /* except on last j */
  135.    for ( ix=j+1; ix<n; ix++ )   /* lower diagonal */
  136.     {
  137.     int ixj = (n*j)+ix;         /*lowr diag rows in col j*/
  138.     for( it=j-ix,jx=j+1; jx<n; jx++ )
  139.       { int ixjx = (n*jx)+ix;
  140.         a[ixjx] -= a[ixj]*a[it+ixjx]; }
  141.     b[ix] -= b[j]*a[ixj];
  142.     } /* --- end-of-for(ix=j+1...n-1) --- */
  143.   } /* --- end-of-for(j=0...n-1) --- */
  144.  
  145. /* --------------------------------------------------------
  146. Back Solution
  147. -------------------------------------------------------- */
  148. for ( it=n*n,i=2; i<=n; i++ )   /* all rows except last */
  149.   {
  150.   int   ia=it-i, ib=n-i, ic=n-1;
  151.   for ( k=1; k<i; k++,ia-=n,ic-- )
  152.     b[ib] -= a[ia]*b[ic];
  153.   } /* --- end-of-for(i=2...n) --- */
  154. return ( 0 );                   /*back to caller with b[]*/
  155. } /* --- end-of-function simq() --- */
  156.  
  157.  
  158. #if TESTDRIVE
  159. /**********************************************************
  160.  *
  161.  * Purpose:     Test driver for simq().
  162.  *              Solves a set of simultaneous equations
  163.  *              of the form ax=b.
  164.  *
  165.  * Arguments:
  166.  *              N (I)   Int containing the number of
  167.  *                      simultaneous equations to be
  168.  *                      solved (must be <=8 for printf
  169.  *                      output to fit on 80-col screen).
  170.  *              a (I)   Double array containing matrix
  171.  *                      of NxN coefficients.  Note that
  172.  *                      a[] is singly subscripted.
  173.  *                      Values are initialized rowwise
  174.  *                      for easy reading, and transposed
  175.  *                      before calling simq() which
  176.  *                      requires columnwise input.
  177.  *              b (I)   Double array containing vector
  178.  *                      of constants.
  179.  *
  180.  * Notes:     o All input is supplied through the
  181.  *              statements immediately below.  Since
  182.  *              this is only a test/demo of simq(),
  183.  *              no real user interface is supplied.
  184.  *
  185.  *********************************************************/
  186.  
  187. /* --------------------------------------------------------
  188. Input data to test simq()
  189. -------------------------------------------------------- */
  190. /* --- The user may change the N,a[],b[] test data
  191.        to solve any set of simultaneous equations
  192.        (but don't try more than 8x8 or the results
  193.        won't printf nicely on an 80-column screen). --- */
  194. #define N 7                     /* N-by-N test data */
  195. /* --- Note: a[] is entered rowwise for easy reading,
  196.        and then transposed before calling simq(). --- */
  197. static double a[] =             /*matrix of coefficients*/
  198.  { 1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  199.    0.0, 1.0, 2.0, 0.0, 0.0, 0.0, 0.0,
  200.    1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
  201.    0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0,
  202.    0.0, 0.0, 0.0, 2.0, 1.0, 0.0, 0.0,
  203.    0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0,
  204.    0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 1.0,
  205.  };
  206. static double b[] =             /* vector of constants */
  207.  { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
  208.  };
  209.  
  210. /* --------------------------------------------------------
  211. Entry point
  212. -------------------------------------------------------- */
  213. main ( argc, argv )
  214. int argc; char *argv[];         /* args not used */
  215. {
  216. /* --- local allocations and declarations --- */
  217. double  ahold[N*N], bhold[N];   /* copy of original data */
  218. int     i,j, n=N;               /* loop indexes and n */
  219.  
  220. /* --- copy input arrays for display later --- */
  221. for ( i=0; i<n; i++ )           /*row# (before transpose)*/
  222.   { int k=n*i;                  /* 1st col of row */
  223.     bhold[i] = b[i];            /* copy constant vector */
  224.     for ( j=0; j<n; j++ )       /* each col of row */
  225.         ahold[k+j] = a[k+j]; }  /*copy coefficient matrix*/
  226.  
  227. /* --- transpose a[] for columnwise input to simq() --- */
  228. for ( i=1; i<n; i++ )           /* each row except 1st */
  229.   for ( j=0; j<i; j++ )         /* upper diagonal */
  230.     {
  231.     int ij = n*i + j;           /* i,j-th element */
  232.     int ji = n*j + i;           /* and its transpose */
  233.     double swap = a[ij];        /* save ... */
  234.     a[ij] = a[ji];  a[ji] = swap; /* ... and swap */
  235.     } /* --- end-of-for(i,j) --- */
  236.  
  237. /* --- call simq to solve --- */
  238. if ( simq(a,b,n) )              /* error if non-0 return */
  239.   { printf("simq() returned error\n");
  240.     exit(0); }
  241.  
  242. /* --- output results --- */
  243. printf("Simq() test results (Ax=b)...\n");
  244. for ( i=0; i<n; i++ )           /* each row */
  245.   {
  246.   /* --- print original coefficients (A) for row --- */
  247.   for ( j=0; j<n; j++ ) printf("%6.2lf",ahold[n*i+j]);
  248.   /* --- print result (x) and constant for row (b) --- */
  249.   printf("  %8.2lf   %8.2lf\n", b[i],bhold[i]);
  250.   } /* --- end-of-for(i) --- */
  251. } /* --- end-of-main() --- */
  252. #endif
  253.