/* 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 (nm1 < 0) goto _70; 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) goto _40; /* Interchange if necessary. */ if (l == k) goto _10; t = a[l][k]; a[l][k] = a[k][k]; a[k][k] = t; _10: /* Compute multipliers. */ if (a[k][k] == 0.0) printf("\n!ERROR. Singular matrix.\m"); 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) goto _20; a[l][j] = a[k][j]; a[k][j] = t; _20: for (i = k+1; i < n; i++ ) a[i][j] += (t * a[i][k]); } goto _50; _40: *info = k; _50: } _70: ipvt[n-1] = n-1; if (a[n-1][n-1] == 0.0) *info = n-1; }