16d84be18SBarry Smith 26d84be18SBarry Smith #ifndef lint 3*39d66777SBarry Smith static char vcid[] = "$Id: dgedi.c,v 1.2 1996/03/04 05:16:18 bsmith Exp bsmith $"; 46d84be18SBarry Smith #endif 56d84be18SBarry Smith 66d84be18SBarry Smith /* 76d84be18SBarry Smith This file creating by running f2c 86d84be18SBarry Smith linpack. this version dated 08/14/78 96d84be18SBarry Smith cleve moler, university of new mexico, argonne national lab. 106d84be18SBarry Smith 116d84be18SBarry Smith Computes the inverse of a matrix given its factors and pivots 126d84be18SBarry Smith calculated by Linpack_DGEFA(). 132bb67c89SBarry Smith */ 142bb67c89SBarry Smith 156d84be18SBarry Smith #include "petsc.h" 162bb67c89SBarry Smith 176d84be18SBarry Smith int Linpack_DGEDI(Scalar *a,int n,int *ipvt,Scalar *work) 182bb67c89SBarry Smith { 19*39d66777SBarry Smith int i__2,kb, kp1, nm1,i, j, k, l, ll,kn,knp1,jn; 206d84be18SBarry Smith Scalar t, *aa,*ax,*ay,tmp; 212bb67c89SBarry Smith 226d84be18SBarry Smith --work; 236d84be18SBarry Smith --ipvt; 24*39d66777SBarry Smith a -= n + 1; 252bb67c89SBarry Smith 262bb67c89SBarry Smith /* compute inverse(u) */ 272bb67c89SBarry Smith 286d84be18SBarry Smith for (k = 1; k <= n; ++k) { 29*39d66777SBarry Smith kn = k*n; 30*39d66777SBarry Smith knp1 = kn + k; 31*39d66777SBarry Smith a[knp1] = 1.0 / a[knp1]; 32*39d66777SBarry Smith t = -a[knp1]; 332bb67c89SBarry Smith i__2 = k - 1; 34*39d66777SBarry Smith aa = &a[1 + kn]; 356d84be18SBarry Smith for ( ll=0; ll<i__2; ll++ ) aa[ll] *= t; 362bb67c89SBarry Smith kp1 = k + 1; 376d84be18SBarry Smith if (n < kp1) continue; 386d84be18SBarry Smith ax = aa; 396d84be18SBarry Smith for (j = kp1; j <= n; ++j) { 40*39d66777SBarry Smith jn = j*n; 41*39d66777SBarry Smith t = a[k + jn]; 42*39d66777SBarry Smith a[k + jn] = 0.; 43*39d66777SBarry Smith ay = &a[1 + jn]; 446d84be18SBarry Smith for ( ll=0; ll<k; ll++ ) { 456d84be18SBarry Smith ay[ll] += t*ax[ll]; 462bb67c89SBarry Smith } 472bb67c89SBarry Smith } 482bb67c89SBarry Smith } 492bb67c89SBarry Smith 502bb67c89SBarry Smith /* form inverse(u)*inverse(l) */ 512bb67c89SBarry Smith 526d84be18SBarry Smith nm1 = n - 1; 532bb67c89SBarry Smith if (nm1 < 1) { 546d84be18SBarry Smith return 0; 552bb67c89SBarry Smith } 566d84be18SBarry Smith for (kb = 1; kb <= nm1; ++kb) { 576d84be18SBarry Smith k = n - kb; 58*39d66777SBarry Smith kn = k*n; 592bb67c89SBarry Smith kp1 = k + 1; 60*39d66777SBarry Smith aa = a + kn; 616d84be18SBarry Smith for (i = kp1; i <= n; ++i) { 626d84be18SBarry Smith work[i] = aa[i]; 636d84be18SBarry Smith aa[i] = 0.; 642bb67c89SBarry Smith } 656d84be18SBarry Smith for (j = kp1; j <= n; ++j) { 662bb67c89SBarry Smith t = work[j]; 676d84be18SBarry Smith ax = &a[j * n + 1]; 68*39d66777SBarry Smith ay = &a[kn + 1]; 696d84be18SBarry Smith for ( ll=0; ll<n; ll++ ) { 706d84be18SBarry Smith ay[ll] += t*ax[ll]; 716d84be18SBarry Smith } 722bb67c89SBarry Smith } 732bb67c89SBarry Smith l = ipvt[k]; 742bb67c89SBarry Smith if (l != k) { 75*39d66777SBarry Smith ax = &a[kn + 1]; 766d84be18SBarry Smith ay = &a[l * n + 1]; 776d84be18SBarry Smith for ( ll=0; ll<n; ll++ ) { 786d84be18SBarry Smith tmp = ax[ll]; 796d84be18SBarry Smith ax[ll] = ay[ll]; 806d84be18SBarry Smith ay[ll] = tmp; 812bb67c89SBarry Smith } 822bb67c89SBarry Smith } 836d84be18SBarry Smith } 842bb67c89SBarry Smith return 0; 856d84be18SBarry Smith } 862bb67c89SBarry Smith 87