1a5eb4965SSatish Balay #ifdef PETSC_RCS_HEADER 2*71c5468dSBarry Smith static char vcid[] = "$Id: dgefa.c,v 1.13 1998/12/24 02:56:18 bsmith Exp bsmith $"; 34431cf12SSatish Balay #endif 49fd41dd0SBarry Smith /* 59fd41dd0SBarry Smith This routine was converted by f2c from Linpack source 69fd41dd0SBarry Smith linpack. this version dated 08/14/78 79fd41dd0SBarry Smith cleve moler, university of new mexico, argonne national lab. 859539b86SBarry Smith 959539b86SBarry Smith Does an LU factorization with partial pivoting of a dense 109d8b60e5SBarry Smith n by n matrix. 119d8b60e5SBarry Smith 129d8b60e5SBarry Smith Used by the sparse factorization routines in 1359539b86SBarry Smith src/mat/impls/baij/seq and src/mat/impls/bdiag/seq 1459539b86SBarry Smith 15*71c5468dSBarry Smith See also src/inline/ilu.h 169fd41dd0SBarry Smith */ 179fd41dd0SBarry Smith #include "petsc.h" 189fd41dd0SBarry Smith 195615d1e5SSatish Balay #undef __FUNC__ 205615d1e5SSatish Balay #define __FUNC__ "Linpack_DGEFA" 213f1db9ecSBarry Smith int Linpack_DGEFA(MatScalar *a, int n, int *ipvt) 229fd41dd0SBarry Smith { 238d3e6ddaSBarry Smith int i__2, i__3, kp1, nm1, j, k, l,ll,kn,knp1,jn1; 243f1db9ecSBarry Smith MatScalar t,*ax,*ay,*aa; 253f1db9ecSBarry Smith MatFloat tmp,max; 269fd41dd0SBarry Smith 279fd41dd0SBarry Smith /* gaussian elimination with partial pivoting */ 289fd41dd0SBarry Smith 293a40ed3dSBarry Smith PetscFunctionBegin; 309fd41dd0SBarry Smith /* Parameter adjustments */ 319fd41dd0SBarry Smith --ipvt; 3239d66777SBarry Smith a -= n + 1; 339fd41dd0SBarry Smith 349fd41dd0SBarry Smith /* Function Body */ 359fd41dd0SBarry Smith nm1 = n - 1; 3639d66777SBarry Smith for (k = 1; k <= nm1; ++k) { 379fd41dd0SBarry Smith kp1 = k + 1; 3839d66777SBarry Smith kn = k*n; 3939d66777SBarry Smith knp1 = k*n + k; 409fd41dd0SBarry Smith 419fd41dd0SBarry Smith /* find l = pivot index */ 429fd41dd0SBarry Smith 439fd41dd0SBarry Smith i__2 = n - k + 1; 4439d66777SBarry Smith aa = &a[knp1]; 459fd41dd0SBarry Smith max = PetscAbsScalar(aa[0]); 469fd41dd0SBarry Smith l = 1; 479fd41dd0SBarry Smith for ( ll=1; ll<i__2; ll++ ) { 489fd41dd0SBarry Smith tmp = PetscAbsScalar(aa[ll]); 499fd41dd0SBarry Smith if (tmp > max) { max = tmp; l = ll+1;} 509fd41dd0SBarry Smith } 519fd41dd0SBarry Smith l += k - 1; 529fd41dd0SBarry Smith ipvt[k] = l; 539fd41dd0SBarry Smith 5439d66777SBarry Smith if (a[l + kn] == 0.) { 55e3372554SBarry Smith SETERRQ(k,0,"Zero pivot"); 569fd41dd0SBarry Smith } 579fd41dd0SBarry Smith 589fd41dd0SBarry Smith /* interchange if necessary */ 599fd41dd0SBarry Smith 609fd41dd0SBarry Smith if (l != k) { 6139d66777SBarry Smith t = a[l + kn]; 6239d66777SBarry Smith a[l + kn] = a[knp1]; 6339d66777SBarry Smith a[knp1] = t; 649fd41dd0SBarry Smith } 659fd41dd0SBarry Smith 669fd41dd0SBarry Smith /* compute multipliers */ 679fd41dd0SBarry Smith 6839d66777SBarry Smith t = -1. / a[knp1]; 699fd41dd0SBarry Smith i__2 = n - k; 7039d66777SBarry Smith aa = &a[1 + knp1]; 719fd41dd0SBarry Smith for ( ll=0; ll<i__2; ll++ ) { 729fd41dd0SBarry Smith aa[ll] *= t; 739fd41dd0SBarry Smith } 749fd41dd0SBarry Smith 759fd41dd0SBarry Smith /* row elimination with column indexing */ 769fd41dd0SBarry Smith 7739d66777SBarry Smith ax = aa; 789fd41dd0SBarry Smith for (j = kp1; j <= n; ++j) { 798d3e6ddaSBarry Smith jn1 = j*n; 808d3e6ddaSBarry Smith t = a[l + jn1]; 819fd41dd0SBarry Smith if (l != k) { 828d3e6ddaSBarry Smith a[l + jn1] = a[k + jn1]; 838d3e6ddaSBarry Smith a[k + jn1] = t; 849fd41dd0SBarry Smith } 859fd41dd0SBarry Smith 869fd41dd0SBarry Smith i__3 = n - k; 878d3e6ddaSBarry Smith ay = &a[1+k+jn1]; 889fd41dd0SBarry Smith for ( ll=0; ll<i__3; ll++ ) { 899fd41dd0SBarry Smith ay[ll] += t*ax[ll]; 909fd41dd0SBarry Smith } 919fd41dd0SBarry Smith } 929fd41dd0SBarry Smith } 939fd41dd0SBarry Smith ipvt[n] = n; 949fd41dd0SBarry Smith if (a[n + n * n] == 0.) { 95e3372554SBarry Smith SETERRQ(n,0,"Zero pivot,final row"); 969fd41dd0SBarry Smith } 973a40ed3dSBarry Smith PetscFunctionReturn(0); 989fd41dd0SBarry Smith } 999fd41dd0SBarry Smith 100