/* dgefa.c -- double precision solution of general matrix by * Gaussian elimination. * * This routine performs the LU factorization. (You should * use dgeco.c if you need a condition estimate.) * * Follow this with a call to dgesl.c, which is the linear * equation solver, or dgedi.c which computes the determinant * and inverse. * * From "Linpack Users' Guide", Dongarra, et al. * Translated from FORTRAN by C. Bond, 1994. */ #include #include #define max(a,b) (((a) > (b)) ? (a) : (b)) #define min(a,b) (((a) < (b)) ? (a) : (b)) void dgefa(double **a,int n, int ipvt[],int *info) { double dmax,t; int i,j,k,kp1,l,nm1; *info = 0; nm1 = n - 1; if (n > 0) { for (k = 0; k < nm1; k++) { kp1 = k + 1; /* Find l = pivot index. */ dmax = fabs(a[k][k]); l = k; for (i = k+1; i < n; i++) { if (fabs(a[i][k]) <= dmax) continue; l = i; } ipvt[k] = l; /* Zero pivot implies this column already triangularized. */ if (a[l][k] == 0.0) { *info = k; continue; } /* Interchange if necessary. */ if (l != k) { t = a[l][k]; a[l][k] = a[k][k]; a[k][k] = t; } /* Compute multipliers. */ if (a[k][k] == 0.0) printf("\n!ERROR. Singular matrix.\n"); t = -1.0/a[k][k]; for (i = k+1; i < n; i++) a[i][k] *= t; /* Row elimination with column indexing. */ for (j = kp1; j < n; j++) { t = a[l][j]; if (l != k) { a[l][j] = a[k][j]; a[k][j] = t; } for (i = k+1; i < n; i++ ) a[i][j] += (t * a[i][k]); } } } ipvt[n-1] = n-1; if (a[n-1][n-1] == 0.0) *info = n-1; }