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