1*be1d678aSKris Buschelman #define PETSCMAT_DLL 2*be1d678aSKris Buschelman 36d84be18SBarry Smith /* 46d84be18SBarry Smith This file creating by running f2c 56d84be18SBarry Smith linpack. this version dated 08/14/78 66d84be18SBarry Smith cleve moler, university of new mexico, argonne national lab. 76d84be18SBarry Smith 86d84be18SBarry Smith Computes the inverse of a matrix given its factors and pivots 9596552b5SBarry Smith calculated by LINPACKdgefa(). Performed in-place for an n by n 109d8b60e5SBarry Smith dense matrix. 119d8b60e5SBarry Smith 129d8b60e5SBarry Smith Used by the sparse factorization routines in 139d8b60e5SBarry Smith src/mat/impls/baij/seq and src/mat/impls/bdiag/seq 1471c5468dSBarry Smith 1571c5468dSBarry Smith See also src/inline/ilu.h 162bb67c89SBarry Smith */ 172bb67c89SBarry Smith 186d84be18SBarry Smith #include "petsc.h" 192bb67c89SBarry Smith 204a2ae208SSatish Balay #undef __FUNCT__ 214a2ae208SSatish Balay #define __FUNCT__ "LINPACKdgedi" 22690b6cddSBarry Smith PetscErrorCode LINPACKdgedi(MatScalar *a,PetscInt n,PetscInt *ipvt,MatScalar *work) 232bb67c89SBarry Smith { 24690b6cddSBarry Smith PetscInt i__2,kb,kp1,nm1,i,j,k,l,ll,kn,knp1,jn1; 253f1db9ecSBarry Smith MatScalar *aa,*ax,*ay,tmp; 263f1db9ecSBarry Smith MatScalar t; 272bb67c89SBarry Smith 283a40ed3dSBarry Smith PetscFunctionBegin; 296d84be18SBarry Smith --work; 306d84be18SBarry Smith --ipvt; 3139d66777SBarry Smith a -= n + 1; 322bb67c89SBarry Smith 332bb67c89SBarry Smith /* compute inverse(u) */ 342bb67c89SBarry Smith 356d84be18SBarry Smith for (k = 1; k <= n; ++k) { 3639d66777SBarry Smith kn = k*n; 3739d66777SBarry Smith knp1 = kn + k; 3839d66777SBarry Smith a[knp1] = 1.0 / a[knp1]; 3939d66777SBarry Smith t = -a[knp1]; 402bb67c89SBarry Smith i__2 = k - 1; 4139d66777SBarry Smith aa = &a[1 + kn]; 426d84be18SBarry Smith for (ll=0; ll<i__2; ll++) aa[ll] *= t; 432bb67c89SBarry Smith kp1 = k + 1; 446d84be18SBarry Smith if (n < kp1) continue; 456d84be18SBarry Smith ax = aa; 466d84be18SBarry Smith for (j = kp1; j <= n; ++j) { 478d3e6ddaSBarry Smith jn1 = j*n; 488d3e6ddaSBarry Smith t = a[k + jn1]; 498d3e6ddaSBarry Smith a[k + jn1] = 0.; 508d3e6ddaSBarry Smith ay = &a[1 + jn1]; 516d84be18SBarry Smith for (ll=0; ll<k; ll++) { 526d84be18SBarry Smith ay[ll] += t*ax[ll]; 532bb67c89SBarry Smith } 542bb67c89SBarry Smith } 552bb67c89SBarry Smith } 562bb67c89SBarry Smith 572bb67c89SBarry Smith /* form inverse(u)*inverse(l) */ 582bb67c89SBarry Smith 596d84be18SBarry Smith nm1 = n - 1; 602bb67c89SBarry Smith if (nm1 < 1) { 613a40ed3dSBarry Smith PetscFunctionReturn(0); 622bb67c89SBarry Smith } 636d84be18SBarry Smith for (kb = 1; kb <= nm1; ++kb) { 646d84be18SBarry Smith k = n - kb; 6539d66777SBarry Smith kn = k*n; 662bb67c89SBarry Smith kp1 = k + 1; 6739d66777SBarry Smith aa = a + kn; 686d84be18SBarry Smith for (i = kp1; i <= n; ++i) { 696d84be18SBarry Smith work[i] = aa[i]; 706d84be18SBarry Smith aa[i] = 0.; 712bb67c89SBarry Smith } 726d84be18SBarry Smith for (j = kp1; j <= n; ++j) { 732bb67c89SBarry Smith t = work[j]; 746d84be18SBarry Smith ax = &a[j * n + 1]; 7539d66777SBarry Smith ay = &a[kn + 1]; 766d84be18SBarry Smith for (ll=0; ll<n; ll++) { 776d84be18SBarry Smith ay[ll] += t*ax[ll]; 786d84be18SBarry Smith } 792bb67c89SBarry Smith } 802bb67c89SBarry Smith l = ipvt[k]; 812bb67c89SBarry Smith if (l != k) { 8239d66777SBarry Smith ax = &a[kn + 1]; 836d84be18SBarry Smith ay = &a[l * n + 1]; 846d84be18SBarry Smith for (ll=0; ll<n; ll++) { 856d84be18SBarry Smith tmp = ax[ll]; 866d84be18SBarry Smith ax[ll] = ay[ll]; 876d84be18SBarry Smith ay[ll] = tmp; 882bb67c89SBarry Smith } 892bb67c89SBarry Smith } 906d84be18SBarry Smith } 913a40ed3dSBarry Smith PetscFunctionReturn(0); 926d84be18SBarry Smith } 932bb67c89SBarry Smith 94