1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dgedi.c,v 1.11 1998/12/24 04:04:06 bsmith Exp bsmith $"; 3 #endif 4 5 /* 6 This file creating by running f2c 7 linpack. this version dated 08/14/78 8 cleve moler, university of new mexico, argonne national lab. 9 10 Computes the inverse of a matrix given its factors and pivots 11 calculated by LINPACKdgefa(). Performed in-place for an n by n 12 dense matrix. 13 14 Used by the sparse factorization routines in 15 src/mat/impls/baij/seq and src/mat/impls/bdiag/seq 16 17 See also src/inline/ilu.h 18 */ 19 20 #include "petsc.h" 21 22 #undef __FUNC__ 23 #define __FUNC__ "LINPACKdgedi" 24 int LINPACKdgedi(MatScalar *a,int n,int *ipvt,MatScalar *work) 25 { 26 int i__2,kb, kp1, nm1,i, j, k, l, ll,kn,knp1,jn1; 27 MatScalar *aa,*ax,*ay,tmp; 28 MatScalar t; 29 30 PetscFunctionBegin; 31 --work; 32 --ipvt; 33 a -= n + 1; 34 35 /* compute inverse(u) */ 36 37 for (k = 1; k <= n; ++k) { 38 kn = k*n; 39 knp1 = kn + k; 40 a[knp1] = 1.0 / a[knp1]; 41 t = -a[knp1]; 42 i__2 = k - 1; 43 aa = &a[1 + kn]; 44 for ( ll=0; ll<i__2; ll++ ) aa[ll] *= t; 45 kp1 = k + 1; 46 if (n < kp1) continue; 47 ax = aa; 48 for (j = kp1; j <= n; ++j) { 49 jn1 = j*n; 50 t = a[k + jn1]; 51 a[k + jn1] = 0.; 52 ay = &a[1 + jn1]; 53 for ( ll=0; ll<k; ll++ ) { 54 ay[ll] += t*ax[ll]; 55 } 56 } 57 } 58 59 /* form inverse(u)*inverse(l) */ 60 61 nm1 = n - 1; 62 if (nm1 < 1) { 63 PetscFunctionReturn(0); 64 } 65 for (kb = 1; kb <= nm1; ++kb) { 66 k = n - kb; 67 kn = k*n; 68 kp1 = k + 1; 69 aa = a + kn; 70 for (i = kp1; i <= n; ++i) { 71 work[i] = aa[i]; 72 aa[i] = 0.; 73 } 74 for (j = kp1; j <= n; ++j) { 75 t = work[j]; 76 ax = &a[j * n + 1]; 77 ay = &a[kn + 1]; 78 for ( ll=0; ll<n; ll++ ) { 79 ay[ll] += t*ax[ll]; 80 } 81 } 82 l = ipvt[k]; 83 if (l != k) { 84 ax = &a[kn + 1]; 85 ay = &a[l * n + 1]; 86 for ( ll=0; ll<n; ll++ ) { 87 tmp = ax[ll]; 88 ax[ll] = ay[ll]; 89 ay[ll] = tmp; 90 } 91 } 92 } 93 PetscFunctionReturn(0); 94 } 95 96