/* dgesl.c -- double precision solution of general matrix by * Gaussian elimination. * * Solves the system A * X = B or TRANS(A) * X = B. * * Operates on the LU decomposition (factorization) * generated by dgefa.c or dgeco.c. * * From "Linpack Users' Guide", Dongarra, et al. * Translated from FORTRAN by C. Bond, 1994. */ #include #include #include "func.h" #define max(a,b) (((a) > (b)) ? (a) : (b)) #define min(a,b) (((a) < (b)) ? (a) : (b)) void dgesl(double **a,int n,int ipvt[],double b[],int job) { double t; int i,k,kb,l,nm1; nm1 = n - 1; if (job == 0) { /* job = 0, solve a * x = b. * First solve l * y = b. */ if (n > 0) { for (k = 0; k < nm1; k++) { l = ipvt[k]; t = b[l]; if (l != k) { b[l] = b[k]; b[k] = t; } /* saxpy(n-k,t,a(k+1,k),1,b(k+1),1); */ for (i=k+1;i < n;i++) b[i] += (t * a[i][k]); } } /* Now solve u * x = y. */ for (kb = 0; kb < n; kb++) { k = n - kb-1; b[k] /= a[k][k]; t = -b[k]; /* saxpy(k-1,t,a(1,k),1,b(1),1); */ for (i = 0; i < k ; i++) b[i] += (t * a[i][k]); } return; } /* job != 0, solve trans(a) * x = b. * First solve trans(u) * x = y. */ for (k = 0; k < n; k++) { /* t = ddot(k-1,a(1,k),1,b(1),1); */ t = 0; for (i = 0; i < k; i++) t += (a[i][k] * b[i]); b[k] = (b[k] - t)/a[k][k]; } /* Now solve trans(l) * x = y. */ if (n > 0) { for (kb = 0; kb < nm1; kb++) { k = n - 2 - kb; /* b[k] = b[k] + ddot(n-k,a(k+1,k),1,b(k+1),1); */ t = 0; for (i = k+1;i < n; i++) t += (a[i][k] * b[i]); b[k] += t; l = ipvt[k]; if (l != k) { t = b[l]; b[l] = b[k]; b[k] = t; } } } }