19fd41dd0SBarry Smith /* 29fd41dd0SBarry Smith This routine was converted by f2c from Linpack source 39fd41dd0SBarry Smith linpack. this version dated 08/14/78 49fd41dd0SBarry Smith cleve moler, university of new mexico, argonne national lab. 559539b86SBarry Smith 659539b86SBarry Smith Does an LU factorization with partial pivoting of a dense 79d8b60e5SBarry Smith n by n matrix. 89d8b60e5SBarry Smith 99d8b60e5SBarry Smith Used by the sparse factorization routines in 1059539b86SBarry Smith src/mat/impls/baij/seq and src/mat/impls/bdiag/seq 1159539b86SBarry Smith 1271c5468dSBarry Smith See also src/inline/ilu.h 139fd41dd0SBarry Smith */ 149fd41dd0SBarry Smith #include "petsc.h" 159fd41dd0SBarry Smith 164a2ae208SSatish Balay #undef __FUNCT__ 174a2ae208SSatish Balay #define __FUNCT__ "LINPACKdgefa" 18596552b5SBarry Smith int LINPACKdgefa(MatScalar *a,int n,int *ipvt) 199fd41dd0SBarry Smith { 208d3e6ddaSBarry Smith int i__2,i__3,kp1,nm1,j,k,l,ll,kn,knp1,jn1; 213f1db9ecSBarry Smith MatScalar t,*ax,*ay,*aa; 22329f5518SBarry Smith MatReal tmp,max; 239fd41dd0SBarry Smith 249fd41dd0SBarry Smith /* gaussian elimination with partial pivoting */ 259fd41dd0SBarry Smith 263a40ed3dSBarry Smith PetscFunctionBegin; 279fd41dd0SBarry Smith /* Parameter adjustments */ 289fd41dd0SBarry Smith --ipvt; 2939d66777SBarry Smith a -= n + 1; 309fd41dd0SBarry Smith 319fd41dd0SBarry Smith /* Function Body */ 329fd41dd0SBarry Smith nm1 = n - 1; 3339d66777SBarry Smith for (k = 1; k <= nm1; ++k) { 349fd41dd0SBarry Smith kp1 = k + 1; 3539d66777SBarry Smith kn = k*n; 3639d66777SBarry Smith knp1 = k*n + k; 379fd41dd0SBarry Smith 389fd41dd0SBarry Smith /* find l = pivot index */ 399fd41dd0SBarry Smith 409fd41dd0SBarry Smith i__2 = n - k + 1; 4139d66777SBarry Smith aa = &a[knp1]; 429fd41dd0SBarry Smith max = PetscAbsScalar(aa[0]); 439fd41dd0SBarry Smith l = 1; 449fd41dd0SBarry Smith for (ll=1; ll<i__2; ll++) { 459fd41dd0SBarry Smith tmp = PetscAbsScalar(aa[ll]); 469fd41dd0SBarry Smith if (tmp > max) { max = tmp; l = ll+1;} 479fd41dd0SBarry Smith } 489fd41dd0SBarry Smith l += k - 1; 499fd41dd0SBarry Smith ipvt[k] = l; 509fd41dd0SBarry Smith 51*5b8514ebSBarry Smith if (a[l + kn] == 0.0) { 5229bbc08cSBarry Smith SETERRQ(k,"Zero pivot"); 539fd41dd0SBarry Smith } 549fd41dd0SBarry Smith 559fd41dd0SBarry Smith /* interchange if necessary */ 569fd41dd0SBarry Smith 579fd41dd0SBarry Smith if (l != k) { 5839d66777SBarry Smith t = a[l + kn]; 5939d66777SBarry Smith a[l + kn] = a[knp1]; 6039d66777SBarry Smith a[knp1] = t; 619fd41dd0SBarry Smith } 629fd41dd0SBarry Smith 639fd41dd0SBarry Smith /* compute multipliers */ 649fd41dd0SBarry Smith 6539d66777SBarry Smith t = -1. / a[knp1]; 669fd41dd0SBarry Smith i__2 = n - k; 6739d66777SBarry Smith aa = &a[1 + knp1]; 689fd41dd0SBarry Smith for (ll=0; ll<i__2; ll++) { 699fd41dd0SBarry Smith aa[ll] *= t; 709fd41dd0SBarry Smith } 719fd41dd0SBarry Smith 729fd41dd0SBarry Smith /* row elimination with column indexing */ 739fd41dd0SBarry Smith 7439d66777SBarry Smith ax = aa; 759fd41dd0SBarry Smith for (j = kp1; j <= n; ++j) { 768d3e6ddaSBarry Smith jn1 = j*n; 778d3e6ddaSBarry Smith t = a[l + jn1]; 789fd41dd0SBarry Smith if (l != k) { 798d3e6ddaSBarry Smith a[l + jn1] = a[k + jn1]; 808d3e6ddaSBarry Smith a[k + jn1] = t; 819fd41dd0SBarry Smith } 829fd41dd0SBarry Smith 839fd41dd0SBarry Smith i__3 = n - k; 848d3e6ddaSBarry Smith ay = &a[1+k+jn1]; 859fd41dd0SBarry Smith for (ll=0; ll<i__3; ll++) { 869fd41dd0SBarry Smith ay[ll] += t*ax[ll]; 879fd41dd0SBarry Smith } 889fd41dd0SBarry Smith } 899fd41dd0SBarry Smith } 909fd41dd0SBarry Smith ipvt[n] = n; 91*5b8514ebSBarry Smith if (a[n + n * n] == 0.0) { 9229bbc08cSBarry Smith SETERRQ(n,"Zero pivot,final row"); 939fd41dd0SBarry Smith } 943a40ed3dSBarry Smith PetscFunctionReturn(0); 959fd41dd0SBarry Smith } 969fd41dd0SBarry Smith 97