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 { 119fd41dd0SBarry Smith int a_offset, i__1, i__2, i__3, kp1, nm1, j, k, l,ll; 12*6d84be18SBarry Smith Scalar t,*aa,*ax,*ay; 13*6d84be18SBarry Smith double tmp,max; 149fd41dd0SBarry Smith 159fd41dd0SBarry Smith /* gaussian elimination with partial pivoting */ 169fd41dd0SBarry Smith 179fd41dd0SBarry Smith /* Parameter adjustments */ 189fd41dd0SBarry Smith --ipvt; 199fd41dd0SBarry Smith a_offset = n + 1; 209fd41dd0SBarry Smith a -= a_offset; 219fd41dd0SBarry Smith 229fd41dd0SBarry Smith /* Function Body */ 239fd41dd0SBarry Smith nm1 = n - 1; 249fd41dd0SBarry Smith if (nm1 < 1) { 259fd41dd0SBarry Smith goto L70; 269fd41dd0SBarry Smith } 279fd41dd0SBarry Smith i__1 = nm1; 289fd41dd0SBarry Smith for (k = 1; k <= i__1; ++k) { 299fd41dd0SBarry Smith kp1 = k + 1; 309fd41dd0SBarry Smith 319fd41dd0SBarry Smith /* find l = pivot index */ 329fd41dd0SBarry Smith 339fd41dd0SBarry Smith i__2 = n - k + 1; 349fd41dd0SBarry Smith /* l = idamax_(&i__2, &a[k + k * n], &c__1) + k - 1; */ 359fd41dd0SBarry Smith aa = &a[k + k * n]; 369fd41dd0SBarry Smith max = PetscAbsScalar(aa[0]); 379fd41dd0SBarry Smith l = 1; 389fd41dd0SBarry Smith for ( ll=1; ll<i__2; ll++ ) { 399fd41dd0SBarry Smith tmp = PetscAbsScalar(aa[ll]); 409fd41dd0SBarry Smith if (tmp > max) { max = tmp; l = ll+1;} 419fd41dd0SBarry Smith } 429fd41dd0SBarry Smith l += k - 1; 439fd41dd0SBarry Smith ipvt[k] = l; 449fd41dd0SBarry Smith 459fd41dd0SBarry Smith /* zero pivot implies this column already triangularized */ 469fd41dd0SBarry Smith 479fd41dd0SBarry Smith if (a[l + k * n] == 0.) { 489fd41dd0SBarry Smith SETERRQ(k,"Linpack_DGEFA:Zero pivot"); 499fd41dd0SBarry Smith } 509fd41dd0SBarry Smith 519fd41dd0SBarry Smith /* interchange if necessary */ 529fd41dd0SBarry Smith 539fd41dd0SBarry Smith if (l != k) { 549fd41dd0SBarry Smith t = a[l + k * n]; 559fd41dd0SBarry Smith a[l + k * n] = a[k + k * n]; 569fd41dd0SBarry Smith a[k + k * n] = t; 579fd41dd0SBarry Smith } 589fd41dd0SBarry Smith 599fd41dd0SBarry Smith /* compute multipliers */ 609fd41dd0SBarry Smith 619fd41dd0SBarry Smith t = -1. / a[k + k * n]; 629fd41dd0SBarry Smith i__2 = n - k; 639fd41dd0SBarry Smith /* dscal_(&i__2, &t, &a[k + 1 + k * n], &c__1); */ 649fd41dd0SBarry Smith aa = &a[k + 1 + k * n]; 659fd41dd0SBarry Smith for ( ll=0; ll<i__2; ll++ ) { 669fd41dd0SBarry Smith aa[ll] *= t; 679fd41dd0SBarry Smith } 689fd41dd0SBarry Smith 699fd41dd0SBarry Smith 709fd41dd0SBarry Smith /* row elimination with column indexing */ 719fd41dd0SBarry Smith 729fd41dd0SBarry Smith ax = &a[k+1+k*n]; 739fd41dd0SBarry Smith for (j = kp1; j <= n; ++j) { 749fd41dd0SBarry Smith t = a[l + j * n]; 759fd41dd0SBarry Smith if (l != k) { 769fd41dd0SBarry Smith a[l + j * n] = a[k + j * n]; 779fd41dd0SBarry Smith a[k + j * n] = t; 789fd41dd0SBarry Smith } 799fd41dd0SBarry Smith 809fd41dd0SBarry Smith i__3 = n - k; 819fd41dd0SBarry Smith /* daxpy_(&i__3, &t, &a[k + 1 + k * n], &c__1, &a[k + 1 + j * 829fd41dd0SBarry Smith n], &c__1); */ 839fd41dd0SBarry Smith ay = &a[k+1+j*n]; 849fd41dd0SBarry Smith for ( ll=0; ll<i__3; ll++ ) { 859fd41dd0SBarry Smith ay[ll] += t*ax[ll]; 869fd41dd0SBarry Smith } 879fd41dd0SBarry Smith } 889fd41dd0SBarry Smith } 899fd41dd0SBarry Smith L70: 909fd41dd0SBarry Smith ipvt[n] = n; 919fd41dd0SBarry Smith if (a[n + n * n] == 0.) { 929fd41dd0SBarry Smith SETERRQ(n,"Linpack_DGEFA:Zero pivot,final row"); 939fd41dd0SBarry Smith } 949fd41dd0SBarry Smith return 0; 959fd41dd0SBarry Smith } 969fd41dd0SBarry Smith 97