home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / biology / gsrc208a.zip / DGESL.C < prev    next >
C/C++ Source or Header  |  1992-05-01  |  2KB  |  102 lines
  1. extern double ddot();
  2. extern void daxpy();
  3.  
  4. void dgesl( a, n, ipvt, b, job )
  5.  
  6. double **a, *b;
  7. int n, *ipvt, job;
  8.  
  9. /*
  10.    Purpose : dgesl solves the linear system
  11.    a * x = b or Transpose(a) * x = b
  12.    using the factors computed by dgeco or degfa.
  13.  
  14.  
  15.    On Entry :
  16.  
  17.       a    : double matrix of dimension ( n+1, n+1 ),
  18.              the output from dgeco or dgefa.
  19.              The 0-th row and column are not used.
  20.       n    : the row dimension of a.
  21.       ipvt : the pivot vector from degco or dgefa.
  22.       b    : the right hand side vector.
  23.       job  : = 0       to solve a * x = b,
  24.              = nonzero to solve Transpose(a) * x = b.
  25.  
  26.  
  27.    On Return :
  28.  
  29.       b : the solution vector x.
  30.  
  31.  
  32.    Error Condition :
  33.  
  34.       A division by zero will occur if the input factor contains
  35.       a zero on the diagonal.  Technically this indicates
  36.       singularity but it is often caused by improper argments or
  37.       improper setting of the pointers of a.  It will not occur
  38.       if the subroutines are called correctly and if dgeco has
  39.       set rcond > 0 or dgefa has set info = 0.
  40.  
  41.  
  42.    BLAS : daxpy, ddot
  43. */
  44.  
  45. {
  46.    int nm1, k, j;
  47.    double t;
  48.  
  49.    nm1 = n - 1;
  50.  
  51. /*
  52.    Job = 0, solve a * x = b.
  53. */
  54.    if ( job == 0 ) {
  55. /*
  56.    First solve L * y = b.
  57. */
  58.       for ( k = 1 ; k <= n ; k++ ) {
  59.          t = ddot( k-1, a[k], 1, b, 1 );
  60.          b[k] = ( b[k] - t ) / a[k][k];
  61.       }
  62. /*
  63.    Now solve U * x = y.
  64. */
  65.       for ( k = n - 1 ; k >= 1 ; k-- ) {
  66.          b[k] = b[k] + ddot( n-k, a[k]+k, 1, b+k, 1 );
  67.          j = ipvt[k];
  68.          if ( j != k ) {
  69.             t = b[j];
  70.             b[j] = b[k];
  71.             b[k] = t;
  72.          }
  73.       }
  74.       return;
  75.    }
  76.  
  77. /*
  78.    Job = nonzero, solve Transpose(a) * x = b.
  79.  
  80.    First solve Transpose(U) * y = b.
  81. */
  82.    for ( k = 1 ; k <= n - 1 ; k++ ) {
  83.       j = ipvt[k];
  84.       t = b[j];
  85.       if ( j != k ) {
  86.          b[j] = b[k];
  87.          b[k] = t;
  88.       }
  89.       daxpy( n-k, t, a[k]+k, 1, b+k, 1 );
  90.    }
  91. /*
  92.    Now solve Transpose(L) * x = y.
  93. */
  94.    for ( k = n ; k >= 1 ; k-- ) {
  95.       b[k] = b[k] / a[k][k];
  96.       t = -b[k];
  97.       daxpy( k-1, t, a[k], 1, b, 1 );
  98.    }
  99.  
  100. }
  101.  
  102.