1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dgefa.c,v 1.8 1997/07/09 20:55:07 balay Exp bsmith $"; 3 #endif 4 /* 5 This routine was converted by f2c from Linpack source 6 linpack. this version dated 08/14/78 7 cleve moler, university of new mexico, argonne national lab. 8 */ 9 #include "petsc.h" 10 11 #undef __FUNC__ 12 #define __FUNC__ "Linpack_DGEFA" 13 int Linpack_DGEFA(Scalar *a, int n, int *ipvt) 14 { 15 int i__2, i__3, kp1, nm1, j, k, l,ll,kn,knp1,jn; 16 Scalar t,*aa,*ax,*ay; 17 double tmp,max; 18 19 /* gaussian elimination with partial pivoting */ 20 21 PetscFunctionBegin; 22 /* Parameter adjustments */ 23 --ipvt; 24 a -= n + 1; 25 26 /* Function Body */ 27 nm1 = n - 1; 28 for (k = 1; k <= nm1; ++k) { 29 kp1 = k + 1; 30 kn = k*n; 31 knp1 = k*n + k; 32 33 /* find l = pivot index */ 34 35 i__2 = n - k + 1; 36 aa = &a[knp1]; 37 max = PetscAbsScalar(aa[0]); 38 l = 1; 39 for ( ll=1; ll<i__2; ll++ ) { 40 tmp = PetscAbsScalar(aa[ll]); 41 if (tmp > max) { max = tmp; l = ll+1;} 42 } 43 l += k - 1; 44 ipvt[k] = l; 45 46 if (a[l + kn] == 0.) { 47 SETERRQ(k,0,"Zero pivot"); 48 } 49 50 /* interchange if necessary */ 51 52 if (l != k) { 53 t = a[l + kn]; 54 a[l + kn] = a[knp1]; 55 a[knp1] = t; 56 } 57 58 /* compute multipliers */ 59 60 t = -1. / a[knp1]; 61 i__2 = n - k; 62 aa = &a[1 + knp1]; 63 for ( ll=0; ll<i__2; ll++ ) { 64 aa[ll] *= t; 65 } 66 67 /* row elimination with column indexing */ 68 69 ax = aa; 70 for (j = kp1; j <= n; ++j) { 71 jn = j*n; 72 t = a[l + jn]; 73 if (l != k) { 74 a[l + jn] = a[k + jn]; 75 a[k + jn] = t; 76 } 77 78 i__3 = n - k; 79 ay = &a[1+k+jn]; 80 for ( ll=0; ll<i__3; ll++ ) { 81 ay[ll] += t*ax[ll]; 82 } 83 } 84 } 85 ipvt[n] = n; 86 if (a[n + n * n] == 0.) { 87 SETERRQ(n,0,"Zero pivot,final row"); 88 } 89 PetscFunctionReturn(0); 90 } 91 92