2: /*
  3:        This routine was converted by f2c from Linpack source
  4:              linpack. this version dated 08/14/78 
  5:       cleve moler, university of new mexico, argonne national lab.

  7:         Does an LU factorization with partial pivoting of a dense
  8:      n by n matrix.

 10:        Used by the sparse factorization routines in 
 11:      src/mat/impls/baij/seq

 13: */
 14: #include <petscsys.h>

 18: PetscErrorCode PetscLINPACKgefa(MatScalar *a,PetscInt n,PetscInt *ipvt)
 19: {
 20:     PetscInt   i__2,i__3,kp1,nm1,j,k,l,ll,kn,knp1,jn1;
 21:     MatScalar  t,*ax,*ay,*aa;
 22:     MatReal    tmp,max;

 24: /*     gaussian elimination with partial pivoting */

 27:     /* Parameter adjustments */
 28:     --ipvt;
 29:     a       -= n + 1;

 31:     /* Function Body */
 32:     nm1 = n - 1;
 33:     for (k = 1; k <= nm1; ++k) {
 34:         kp1  = k + 1;
 35:         kn   = k*n;
 36:         knp1 = k*n + k;

 38: /*        find l = pivot index */

 40:         i__2 = n - k + 1;
 41:         aa = &a[knp1];
 42:         max = PetscAbsScalar(aa[0]);
 43:         l = 1;
 44:         for (ll=1; ll<i__2; ll++) {
 45:           tmp = PetscAbsScalar(aa[ll]);
 46:           if (tmp > max) { max = tmp; l = ll+1;}
 47:         }
 48:         l += k - 1;
 49:         ipvt[k] = l;

 51:         if (a[l + kn] == 0.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",k-1);

 53: /*           interchange if necessary */

 55:         if (l != k) {
 56:           t = a[l + kn];
 57:           a[l + kn] = a[knp1];
 58:           a[knp1] = t;
 59:         }

 61: /*           compute multipliers */

 63:         t = -1. / a[knp1];
 64:         i__2 = n - k;
 65:         aa = &a[1 + knp1];
 66:         for (ll=0; ll<i__2; ll++) {
 67:           aa[ll] *= t;
 68:         }

 70: /*           row elimination with column indexing */

 72:         ax = aa;
 73:         for (j = kp1; j <= n; ++j) {
 74:             jn1 = j*n;
 75:             t = a[l + jn1];
 76:             if (l != k) {
 77:               a[l + jn1] = a[k + jn1];
 78:               a[k + jn1] = t;
 79:             }

 81:             i__3 = n - k;
 82:             ay = &a[1+k+jn1];
 83:             for (ll=0; ll<i__3; ll++) {
 84:               ay[ll] += t*ax[ll];
 85:             }
 86:         }
 87:     }
 88:     ipvt[n] = n;
 89:     if (a[n + n * n] == 0.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",n-1);
 90:     return(0);
 91: }