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: PETSC_INTERN PetscErrorCode PetscLINPACKgefa(MatScalar *a,PetscInt n,PetscInt *ipvt,PetscBool allowzeropivot,PetscBool *zeropivotdetected)
 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;

 25:   if (zeropivotdetected) *zeropivotdetected = PETSC_FALSE;

 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) {
 52:       if (allowzeropivot) {
 54:         PetscInfo1(NULL,"Zero pivot, row %D\n",k-1);
 55:         if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE;
 56:       } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",k-1);
 57:     }

 59:     /* interchange if necessary */
 60:     if (l != k) {
 61:       t         = a[l + kn];
 62:       a[l + kn] = a[knp1];
 63:       a[knp1]   = t;
 64:     }

 66:     /* compute multipliers */
 67:     t    = -1. / a[knp1];
 68:     i__2 = n - k;
 69:     aa   = &a[1 + knp1];
 70:     for (ll=0; ll<i__2; ll++) aa[ll] *= t;

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

 82:       i__3 = n - k;
 83:       ay   = &a[1+k+jn1];
 84:       for (ll=0; ll<i__3; ll++) ay[ll] += t*ax[ll];
 85:     }
 86:   }
 87:   ipvt[n] = n;
 88:   if (a[n + n * n] == 0.0) {
 89:     if (allowzeropivot) {
 91:       PetscInfo1(NULL,"Zero pivot, row %D\n",n-1);
 92:       if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE;
 93:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",n-1);
 94:   }
 95:   return(0);
 96: }