Actual source code: dgefa.c

      1: #define PETSCMAT_DLL

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

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

     11:        Used by the sparse factorization routines in
     12:      src/mat/impls/baij/seq and src/mat/impls/bdiag/seq

     14:        See also src/inline/ilu.h
     15: */
     16:  #include petsc.h

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

     26: /*     gaussian elimination with partial pivoting */

     29:     /* Parameter adjustments */
     30:     --ipvt;
     31:     a       -= n + 1;

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

     40: /*        find l = pivot index */

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

     53:         if (a[l + kn] == 0.0) {
     54:           SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D",k-1);
     55:         }

     57: /*           interchange if necessary */

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

     65: /*           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++) {
     71:           aa[ll] *= t;
     72:         }

     74: /*           row elimination with column indexing */

     76:         ax = aa;
     77:         for (j = kp1; j <= n; ++j) {
     78:             jn1 = j*n;
     79:             t = a[l + jn1];
     80:             if (l != k) {
     81:               a[l + jn1] = a[k + jn1];
     82:               a[k + jn1] = t;
     83:             }

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