home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / biology / gsrc208a.zip / DGEFA.C < prev    next >
C/C++ Source or Header  |  1992-05-01  |  3KB  |  101 lines
  1. extern void daxpy(), dscal();
  2. extern int idamax();
  3.  
  4. void dgefa( a, n, ipvt, info )
  5.  
  6. double **a;
  7. int n, *ipvt, *info;
  8.  
  9. /*
  10.    Purpose : dgefa factors a double matrix by Gaussian elimination.
  11.  
  12.    dgefa is usually called by dgeco, but it can be called directly
  13.    with a saving in time if rcond is not needed.
  14.    (Time for dgeco) = (1+9/n)*(time for dgefa).
  15.  
  16.    This c version uses algorithm kji rather than the kij in dgefa.f.
  17.    Note that the fortran version input variable lda is not needed.
  18.  
  19.  
  20.    On Entry :
  21.  
  22.       a   : double matrix of dimension ( n+1, n+1 ),
  23.             the 0-th row and column are not used.
  24.             a is created using NewDoubleMatrix, hence
  25.             lda is unnecessary.
  26.       n   : the row dimension of a.
  27.  
  28.    On Return :
  29.  
  30.       a     : a lower triangular matrix and the multipliers
  31.               which were used to obtain it.  The factorization
  32.               can be written a = L * U where U is a product of
  33.               permutation and unit upper triangular matrices
  34.               and L is lower triangular.
  35.       ipvt  : an n+1 integer vector of pivot indices.
  36.       *info : = 0 normal value,
  37.               = k if U[k][k] == 0.  This is not an error
  38.                 condition for this subroutine, but it does
  39.                 indicate that dgesl or dgedi will divide by
  40.                 zero if called.  Use rcond in dgeco for
  41.                 a reliable indication of singularity.
  42.  
  43.                 Notice that the calling program must use &info.
  44.  
  45.    BLAS : daxpy, dscal, idamax
  46. */
  47.  
  48. {
  49.    int j, k, i;
  50.    double t;
  51.  
  52. /* Gaussian elimination with partial pivoting.   */
  53.  
  54.    *info = 0;
  55.    for ( k = 1 ; k <= n - 1 ; k++ ) {
  56. /*
  57.    Find j = pivot index.  Note that a[k]+k-1 is the address of
  58.    the 0-th element of the row vector whose 1st element is a[k][k].
  59. */
  60.       j = idamax( n-k+1, a[k]+k-1, 1 ) + k - 1;      
  61.       ipvt[k] = j;
  62. /*
  63.    Zero pivot implies this row already triangularized.
  64. */
  65.       if ( a[k][j] == 0. ) {
  66.          *info = k;
  67.          continue;
  68.       }      
  69. /*
  70.    Interchange if necessary.
  71. */
  72.       if ( j != k ) {
  73.          t = a[k][j];
  74.          a[k][j] = a[k][k];
  75.          a[k][k] = t;
  76.       }
  77. /*
  78.    Compute multipliers.
  79. */
  80.       t = -1. / a[k][k];
  81.       dscal( n-k, t, a[k]+k, 1 );
  82. /*
  83.    Column elimination with row indexing.
  84. */
  85.       for ( i = k + 1 ; i <= n ; i++ ) {
  86.          t = a[i][j];
  87.          if ( j != k ) {
  88.             a[i][j] = a[i][k];
  89.             a[i][k] = t;
  90.          }
  91.          daxpy( n-k, t, a[k]+k, 1, a[i]+k, 1 );
  92.       }
  93.    }                     /*  end k-loop  */
  94.  
  95.    ipvt[n] = n;
  96.    if ( a[n][n] == 0. )
  97.       *info = n;
  98.  
  99. }
  100.  
  101.