1*6d84be18SBarry Smith 2*6d84be18SBarry Smith #ifndef lint 3*6d84be18SBarry Smith static char vcid[] = "$Id: baijfact.c,v 1.6 1996/02/20 18:52:30 curfman Exp bsmith $"; 4*6d84be18SBarry Smith #endif 5*6d84be18SBarry Smith 6*6d84be18SBarry Smith /* 7*6d84be18SBarry Smith This file creating by running f2c 8*6d84be18SBarry Smith linpack. this version dated 08/14/78 9*6d84be18SBarry Smith cleve moler, university of new mexico, argonne national lab. 10*6d84be18SBarry Smith 11*6d84be18SBarry Smith Computes the inverse of a matrix given its factors and pivots 12*6d84be18SBarry Smith calculated by Linpack_DGEFA(). 132bb67c89SBarry Smith */ 142bb67c89SBarry Smith 15*6d84be18SBarry Smith #include "petsc.h" 162bb67c89SBarry Smith 17*6d84be18SBarry Smith int Linpack_DGEDI(Scalar *a,int n,int *ipvt,Scalar *work) 182bb67c89SBarry Smith { 19*6d84be18SBarry Smith int a_offset, i__2,kb, kp1, nm1,i, j, k, l, ll; 20*6d84be18SBarry Smith Scalar t, *aa,*ax,*ay,tmp; 212bb67c89SBarry Smith 22*6d84be18SBarry Smith --work; 23*6d84be18SBarry Smith --ipvt; 24*6d84be18SBarry Smith a_offset = n + 1; 25*6d84be18SBarry Smith a -= a_offset; 262bb67c89SBarry Smith 272bb67c89SBarry Smith /* compute inverse(u) */ 282bb67c89SBarry Smith 29*6d84be18SBarry Smith for (k = 1; k <= n; ++k) { 30*6d84be18SBarry Smith a[k + k * n] = 1. / a[k + k * n]; 31*6d84be18SBarry Smith t = -a[k + k * n]; 322bb67c89SBarry Smith i__2 = k - 1; 33*6d84be18SBarry Smith /* dscal_(&i__2, &t, &a[k * n + 1], &c__1); */ 34*6d84be18SBarry Smith aa = &a[k * n + 1]; 35*6d84be18SBarry Smith for ( ll=0; ll<i__2; ll++ ) aa[ll] *= t; 362bb67c89SBarry Smith kp1 = k + 1; 37*6d84be18SBarry Smith if (n < kp1) continue; 38*6d84be18SBarry Smith ax = aa; 39*6d84be18SBarry Smith for (j = kp1; j <= n; ++j) { 40*6d84be18SBarry Smith t = a[k + j * n]; 41*6d84be18SBarry Smith a[k + j * n] = 0.; 42*6d84be18SBarry Smith /* daxpy_(&k, &t, &a[k * n + 1], &c__1, &a[j * n + 1], &c__1);*/ 43*6d84be18SBarry Smith ay = &a[j * n + 1]; 44*6d84be18SBarry Smith for ( ll=0; ll<k; ll++ ) { 45*6d84be18SBarry Smith ay[ll] += t*ax[ll]; 462bb67c89SBarry Smith } 472bb67c89SBarry Smith } 482bb67c89SBarry Smith } 492bb67c89SBarry Smith 502bb67c89SBarry Smith /* form inverse(u)*inverse(l) */ 512bb67c89SBarry Smith 52*6d84be18SBarry Smith nm1 = n - 1; 532bb67c89SBarry Smith if (nm1 < 1) { 54*6d84be18SBarry Smith return 0; 552bb67c89SBarry Smith } 56*6d84be18SBarry Smith for (kb = 1; kb <= nm1; ++kb) { 57*6d84be18SBarry Smith k = n - kb; 582bb67c89SBarry Smith kp1 = k + 1; 59*6d84be18SBarry Smith aa = a + k * n; 60*6d84be18SBarry Smith for (i = kp1; i <= n; ++i) { 61*6d84be18SBarry Smith work[i] = aa[i]; 62*6d84be18SBarry Smith aa[i] = 0.; 632bb67c89SBarry Smith } 64*6d84be18SBarry Smith for (j = kp1; j <= n; ++j) { 652bb67c89SBarry Smith t = work[j]; 66*6d84be18SBarry Smith /* daxpy_(n, &t, &a[j * n + 1], &c__1, &a[k * n + 1], &c__1);*/ 67*6d84be18SBarry Smith ax = &a[j * n + 1]; 68*6d84be18SBarry Smith ay = &a[k * n + 1]; 69*6d84be18SBarry Smith for ( ll=0; ll<n; ll++ ) { 70*6d84be18SBarry Smith ay[ll] += t*ax[ll]; 71*6d84be18SBarry Smith } 722bb67c89SBarry Smith } 732bb67c89SBarry Smith l = ipvt[k]; 742bb67c89SBarry Smith if (l != k) { 75*6d84be18SBarry Smith /* dswap_(n, &a[k * n + 1], &c__1, &a[l * n + 1], &c__1); */ 76*6d84be18SBarry Smith ax = &a[k * n + 1]; 77*6d84be18SBarry Smith ay = &a[l * n + 1]; 78*6d84be18SBarry Smith for ( ll=0; ll<n; ll++ ) { 79*6d84be18SBarry Smith tmp = ax[ll]; 80*6d84be18SBarry Smith ax[ll] = ay[ll]; 81*6d84be18SBarry Smith ay[ll] = tmp; 822bb67c89SBarry Smith } 832bb67c89SBarry Smith } 84*6d84be18SBarry Smith } 852bb67c89SBarry Smith return 0; 86*6d84be18SBarry Smith } 872bb67c89SBarry Smith 88