19fd41dd0SBarry Smith 29fd41dd0SBarry Smith /* 39fd41dd0SBarry Smith This routine was converted by f2c from Linpack source 49fd41dd0SBarry Smith linpack. this version dated 08/14/78 59fd41dd0SBarry Smith cleve moler, university of new mexico, argonne national lab. 69fd41dd0SBarry Smith */ 79fd41dd0SBarry Smith #include "petsc.h" 89fd41dd0SBarry Smith 99fd41dd0SBarry Smith int Linpack_DGEFA(Scalar *a, int n, int *ipvt) 109fd41dd0SBarry Smith { 11*39d66777SBarry Smith int i__2, i__3, kp1, nm1, j, k, l,ll,kn,knp1,jn; 126d84be18SBarry Smith Scalar t,*aa,*ax,*ay; 136d84be18SBarry Smith double tmp,max; 149fd41dd0SBarry Smith 159fd41dd0SBarry Smith /* gaussian elimination with partial pivoting */ 169fd41dd0SBarry Smith 179fd41dd0SBarry Smith /* Parameter adjustments */ 189fd41dd0SBarry Smith --ipvt; 19*39d66777SBarry Smith a -= n + 1; 209fd41dd0SBarry Smith 219fd41dd0SBarry Smith /* Function Body */ 229fd41dd0SBarry Smith nm1 = n - 1; 23*39d66777SBarry Smith for (k = 1; k <= nm1; ++k) { 249fd41dd0SBarry Smith kp1 = k + 1; 25*39d66777SBarry Smith kn = k*n; 26*39d66777SBarry Smith knp1 = k*n + k; 279fd41dd0SBarry Smith 289fd41dd0SBarry Smith /* find l = pivot index */ 299fd41dd0SBarry Smith 309fd41dd0SBarry Smith i__2 = n - k + 1; 31*39d66777SBarry Smith aa = &a[knp1]; 329fd41dd0SBarry Smith max = PetscAbsScalar(aa[0]); 339fd41dd0SBarry Smith l = 1; 349fd41dd0SBarry Smith for ( ll=1; ll<i__2; ll++ ) { 359fd41dd0SBarry Smith tmp = PetscAbsScalar(aa[ll]); 369fd41dd0SBarry Smith if (tmp > max) { max = tmp; l = ll+1;} 379fd41dd0SBarry Smith } 389fd41dd0SBarry Smith l += k - 1; 399fd41dd0SBarry Smith ipvt[k] = l; 409fd41dd0SBarry Smith 41*39d66777SBarry Smith if (a[l + kn] == 0.) { 429fd41dd0SBarry Smith SETERRQ(k,"Linpack_DGEFA:Zero pivot"); 439fd41dd0SBarry Smith } 449fd41dd0SBarry Smith 459fd41dd0SBarry Smith /* interchange if necessary */ 469fd41dd0SBarry Smith 479fd41dd0SBarry Smith if (l != k) { 48*39d66777SBarry Smith t = a[l + kn]; 49*39d66777SBarry Smith a[l + kn] = a[knp1]; 50*39d66777SBarry Smith a[knp1] = t; 519fd41dd0SBarry Smith } 529fd41dd0SBarry Smith 539fd41dd0SBarry Smith /* compute multipliers */ 549fd41dd0SBarry Smith 55*39d66777SBarry Smith t = -1. / a[knp1]; 569fd41dd0SBarry Smith i__2 = n - k; 57*39d66777SBarry Smith aa = &a[1 + knp1]; 589fd41dd0SBarry Smith for ( ll=0; ll<i__2; ll++ ) { 599fd41dd0SBarry Smith aa[ll] *= t; 609fd41dd0SBarry Smith } 619fd41dd0SBarry Smith 629fd41dd0SBarry Smith /* row elimination with column indexing */ 639fd41dd0SBarry Smith 64*39d66777SBarry Smith ax = aa; 659fd41dd0SBarry Smith for (j = kp1; j <= n; ++j) { 66*39d66777SBarry Smith jn = j*n; 67*39d66777SBarry Smith t = a[l + jn]; 689fd41dd0SBarry Smith if (l != k) { 69*39d66777SBarry Smith a[l + jn] = a[k + jn]; 70*39d66777SBarry Smith a[k + jn] = t; 719fd41dd0SBarry Smith } 729fd41dd0SBarry Smith 739fd41dd0SBarry Smith i__3 = n - k; 74*39d66777SBarry Smith ay = &a[1+k+jn]; 759fd41dd0SBarry Smith for ( ll=0; ll<i__3; ll++ ) { 769fd41dd0SBarry Smith ay[ll] += t*ax[ll]; 779fd41dd0SBarry Smith } 789fd41dd0SBarry Smith } 799fd41dd0SBarry Smith } 809fd41dd0SBarry Smith ipvt[n] = n; 819fd41dd0SBarry Smith if (a[n + n * n] == 0.) { 829fd41dd0SBarry Smith SETERRQ(n,"Linpack_DGEFA:Zero pivot,final row"); 839fd41dd0SBarry Smith } 849fd41dd0SBarry Smith return 0; 859fd41dd0SBarry Smith } 869fd41dd0SBarry Smith 87