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