/* dgedi.c -- double precision solution of general matrix by * Gaussian elimination. computes determinant and * inverse of matrix (following dgeco or dgefa). * * From "Linpack Users' Guide", Dongarra, et al., */ #include #include #include "func.h" #define max(a,b) (((a) > (b)) ? (a) : (b)) #define min(a,b) (((a) < (b)) ? (a) : (b)) void dgedi(double **a,int n,int ipvt[],double det[],double work[],int job) { double t,ten; int i,j,k,kb,kp1,l,nm1; /* Compute determinant. */ if ((job / 10) != 0) { det[0] = 1.0; det[1] = 0.0; ten = 10.0; for (i = 0; i < n; i++) { if (ipvt[i] != i) det[0] = -det[0]; det[0] *= a[i][i]; if (det[0] != 0.0) { while (fabs(det[0]) < 1.0) { det[0] *= ten; det[1] -= 1.0; } while (fabs(det[0]) >= ten) { det[0] /= ten; det[1] += 1.0; } } } } /* Compute inverse. */ if (job == 10) return; for (k = 0; k < n; k++) { a[k][k] = 1.0 / a[k][k]; t = -a[k][k]; /* dscal(k-1,t,a(1,k),1) */ for (i = 0; i < k; i++) a[i][k] *= t; kp1 = k + 1; if ( n > kp1) { for (j = kp1; j < n; j++) { t = a[k][j]; a[k][j] = 0.0; /* daxpy(k,t,a(1,k),1,a(1,j),1) */ for (i = 0; i < k+1; i++) a[i][j] += (t * a[i][k]); } } } /* Form inverse(u) * inverse(l) */ nm1 = n - 1; if (n > 0) { for (kb = 0; kb < nm1; kb++) { k = n - kb - 2; kp1 = k + 1; for (i = kp1; i < n; i++) { work[i] = a[i][k]; a[i][k] = 0.0; } for (j = kp1; j < n; j++) { t = work[j]; /* daxpy(n,t,a(1,j),1,a(1,k),1) */ for (i = 0; i < n; i++) a[i][k] += (t * a[i][j]); } l = ipvt[k]; if (l != k) { /* dswap(n,a(1,k),1,a(1,l),1) */ for (i = 0; i < n; i++) { t = a[i][k]; a[i][k] = a[i][l]; a[i][l] = t; } } } } }