1 #ifndef lint 2 static char vcid[] = "$Id: dpause.c,v 1.7 1996/12/18 20:55:58 balay Exp $"; 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 __FUNCTION__ 12 #define __FUNCTION__ "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 /* Parameter adjustments */ 22 --ipvt; 23 a -= n + 1; 24 25 /* Function Body */ 26 nm1 = n - 1; 27 for (k = 1; k <= nm1; ++k) { 28 kp1 = k + 1; 29 kn = k*n; 30 knp1 = k*n + k; 31 32 /* find l = pivot index */ 33 34 i__2 = n - k + 1; 35 aa = &a[knp1]; 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 if (a[l + kn] == 0.) { 46 SETERRQ(k,"Linpack_DGEFA:Zero pivot"); 47 } 48 49 /* interchange if necessary */ 50 51 if (l != k) { 52 t = a[l + kn]; 53 a[l + kn] = a[knp1]; 54 a[knp1] = t; 55 } 56 57 /* compute multipliers */ 58 59 t = -1. / a[knp1]; 60 i__2 = n - k; 61 aa = &a[1 + knp1]; 62 for ( ll=0; ll<i__2; ll++ ) { 63 aa[ll] *= t; 64 } 65 66 /* row elimination with column indexing */ 67 68 ax = aa; 69 for (j = kp1; j <= n; ++j) { 70 jn = j*n; 71 t = a[l + jn]; 72 if (l != k) { 73 a[l + jn] = a[k + jn]; 74 a[k + jn] = t; 75 } 76 77 i__3 = n - k; 78 ay = &a[1+k+jn]; 79 for ( ll=0; ll<i__3; ll++ ) { 80 ay[ll] += t*ax[ll]; 81 } 82 } 83 } 84 ipvt[n] = n; 85 if (a[n + n * n] == 0.) { 86 SETERRQ(n,"Linpack_DGEFA:Zero pivot,final row"); 87 } 88 return 0; 89 } 90 91