Actual source code: dgedi.c

      1: #define PETSCMAT_DLL

      3: /*
      4:               This file creating by running f2c
      5:             linpack. this version dated 08/14/78
      6:       cleve moler, university of new mexico, argonne national lab.

      8:       Computes the inverse of a matrix given its factors and pivots
      9:     calculated by LINPACKdgefa(). Performed in-place for an n by n
     10:     dense matrix.

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

     15:        See also src/inline/ilu.h
     16: */

     18:  #include petsc.h

     22: PetscErrorCode LINPACKdgedi(MatScalar *a,PetscInt n,PetscInt *ipvt,MatScalar *work)
     23: {
     24:     PetscInt   i__2,kb,kp1,nm1,i,j,k,l,ll,kn,knp1,jn1;
     25:     MatScalar  *aa,*ax,*ay,tmp;
     26:     MatScalar  t;

     29:     --work;
     30:     --ipvt;
     31:     a       -= n + 1;

     33:    /*     compute inverse(u) */

     35:     for (k = 1; k <= n; ++k) {
     36:         kn           = k*n;
     37:         knp1         = kn + k;
     38:         a[knp1]      = 1.0 / a[knp1];
     39:         t            = -a[knp1];
     40:         i__2         = k - 1;
     41:         aa           = &a[1 + kn];
     42:         for (ll=0; ll<i__2; ll++) aa[ll] *= t;
     43:         kp1 = k + 1;
     44:         if (n < kp1) continue;
     45:         ax = aa;
     46:         for (j = kp1; j <= n; ++j) {
     47:             jn1 = j*n;
     48:             t = a[k + jn1];
     49:             a[k + jn1] = 0.;
     50:             ay = &a[1 + jn1];
     51:             for (ll=0; ll<k; ll++) {
     52:               ay[ll] += t*ax[ll];
     53:             }
     54:         }
     55:     }

     57:    /*    form inverse(u)*inverse(l) */

     59:     nm1 = n - 1;
     60:     if (nm1 < 1) {
     61:         return(0);
     62:     }
     63:     for (kb = 1; kb <= nm1; ++kb) {
     64:         k   = n - kb;
     65:         kn  = k*n;
     66:         kp1 = k + 1;
     67:         aa  = a + kn;
     68:         for (i = kp1; i <= n; ++i) {
     69:             work[i] = aa[i];
     70:             aa[i]   = 0.;
     71:         }
     72:         for (j = kp1; j <= n; ++j) {
     73:             t = work[j];
     74:             ax = &a[j * n + 1];
     75:             ay = &a[kn + 1];
     76:             for (ll=0; ll<n; ll++) {
     77:               ay[ll] += t*ax[ll];
     78:             }
     79:         }
     80:         l = ipvt[k];
     81:         if (l != k) {
     82:             ax = &a[kn + 1];
     83:             ay = &a[l * n + 1];
     84:             for (ll=0; ll<n; ll++) {
     85:               tmp    = ax[ll];
     86:               ax[ll] = ay[ll];
     87:               ay[ll] = tmp;
     88:             }
     89:         }
     90:     }
     91:     return(0);
     92: }