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 a_offset, i__1, i__2, i__3, kp1, nm1, j, k, l,ll; 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_offset = n + 1; 20 a -= a_offset; 21 22 /* Function Body */ 23 nm1 = n - 1; 24 if (nm1 < 1) { 25 goto L70; 26 } 27 i__1 = nm1; 28 for (k = 1; k <= i__1; ++k) { 29 kp1 = k + 1; 30 31 /* find l = pivot index */ 32 33 i__2 = n - k + 1; 34 /* l = idamax_(&i__2, &a[k + k * n], &c__1) + k - 1; */ 35 aa = &a[k + k * n]; 36 max = PetscAbsScalar(aa[0]); 37 l = 1; 38 for ( ll=1; ll<i__2; ll++ ) { 39 tmp = PetscAbsScalar(aa[ll]); 40 if (tmp > max) { max = tmp; l = ll+1;} 41 } 42 l += k - 1; 43 ipvt[k] = l; 44 45 /* zero pivot implies this column already triangularized */ 46 47 if (a[l + k * n] == 0.) { 48 SETERRQ(k,"Linpack_DGEFA:Zero pivot"); 49 } 50 51 /* interchange if necessary */ 52 53 if (l != k) { 54 t = a[l + k * n]; 55 a[l + k * n] = a[k + k * n]; 56 a[k + k * n] = t; 57 } 58 59 /* compute multipliers */ 60 61 t = -1. / a[k + k * n]; 62 i__2 = n - k; 63 /* dscal_(&i__2, &t, &a[k + 1 + k * n], &c__1); */ 64 aa = &a[k + 1 + k * n]; 65 for ( ll=0; ll<i__2; ll++ ) { 66 aa[ll] *= t; 67 } 68 69 70 /* row elimination with column indexing */ 71 72 ax = &a[k+1+k*n]; 73 for (j = kp1; j <= n; ++j) { 74 t = a[l + j * n]; 75 if (l != k) { 76 a[l + j * n] = a[k + j * n]; 77 a[k + j * n] = t; 78 } 79 80 i__3 = n - k; 81 /* daxpy_(&i__3, &t, &a[k + 1 + k * n], &c__1, &a[k + 1 + j * 82 n], &c__1); */ 83 ay = &a[k+1+j*n]; 84 for ( ll=0; ll<i__3; ll++ ) { 85 ay[ll] += t*ax[ll]; 86 } 87 } 88 } 89 L70: 90 ipvt[n] = n; 91 if (a[n + n * n] == 0.) { 92 SETERRQ(n,"Linpack_DGEFA:Zero pivot,final row"); 93 } 94 return 0; 95 } 96 97