16d84be18SBarry Smith /* 26d84be18SBarry Smith This file creating by running f2c 36d84be18SBarry Smith linpack. this version dated 08/14/78 46d84be18SBarry Smith cleve moler, university of new mexico, argonne national lab. 56d84be18SBarry Smith 66d84be18SBarry Smith Computes the inverse of a matrix given its factors and pivots 796b95a6bSBarry Smith calculated by PetscLINPACKgefa(). Performed in-place for an n by n 89d8b60e5SBarry Smith dense matrix. 99d8b60e5SBarry Smith 109d8b60e5SBarry Smith Used by the sparse factorization routines in 11dd882469SBarry Smith src/mat/impls/baij/seq 1271c5468dSBarry Smith 132bb67c89SBarry Smith */ 142bb67c89SBarry Smith 155f748d31SJed Brown #include <petsc/private/kernels/blockinvert.h> 162bb67c89SBarry Smith 17*9de2952eSStefano Zampini PetscErrorCode PetscLINPACKgedi(MatScalar *a, PetscInt n, PetscInt *ipvt, MatScalar *work) 18d71ae5a4SJacob Faibussowitsch { 19690b6cddSBarry Smith PetscInt i__2, kb, kp1, nm1, i, j, k, l, ll, kn, knp1, jn1; 203f1db9ecSBarry Smith MatScalar *aa, *ax, *ay, tmp; 213f1db9ecSBarry Smith MatScalar t; 222bb67c89SBarry Smith 233a40ed3dSBarry Smith PetscFunctionBegin; 246d84be18SBarry Smith --work; 256d84be18SBarry Smith --ipvt; 2639d66777SBarry Smith a -= n + 1; 272bb67c89SBarry Smith 282bb67c89SBarry Smith /* compute inverse(u) */ 292bb67c89SBarry Smith 306d84be18SBarry Smith for (k = 1; k <= n; ++k) { 3139d66777SBarry Smith kn = k * n; 3239d66777SBarry Smith knp1 = kn + k; 3339d66777SBarry Smith a[knp1] = 1.0 / a[knp1]; 3439d66777SBarry Smith t = -a[knp1]; 352bb67c89SBarry Smith i__2 = k - 1; 3639d66777SBarry Smith aa = &a[1 + kn]; 376d84be18SBarry Smith for (ll = 0; ll < i__2; ll++) aa[ll] *= t; 382bb67c89SBarry Smith kp1 = k + 1; 396d84be18SBarry Smith if (n < kp1) continue; 406d84be18SBarry Smith ax = aa; 416d84be18SBarry Smith for (j = kp1; j <= n; ++j) { 428d3e6ddaSBarry Smith jn1 = j * n; 438d3e6ddaSBarry Smith t = a[k + jn1]; 448d3e6ddaSBarry Smith a[k + jn1] = 0.; 458d3e6ddaSBarry Smith ay = &a[1 + jn1]; 4626fbe8dcSKarl Rupp for (ll = 0; ll < k; ll++) ay[ll] += t * ax[ll]; 472bb67c89SBarry Smith } 482bb67c89SBarry Smith } 492bb67c89SBarry Smith 502bb67c89SBarry Smith /* form inverse(u)*inverse(l) */ 512bb67c89SBarry Smith 526d84be18SBarry Smith nm1 = n - 1; 533ba16761SJacob Faibussowitsch if (nm1 < 1) PetscFunctionReturn(PETSC_SUCCESS); 546d84be18SBarry Smith for (kb = 1; kb <= nm1; ++kb) { 556d84be18SBarry Smith k = n - kb; 5639d66777SBarry Smith kn = k * n; 572bb67c89SBarry Smith kp1 = k + 1; 5839d66777SBarry Smith aa = a + kn; 596d84be18SBarry Smith for (i = kp1; i <= n; ++i) { 606d84be18SBarry Smith work[i] = aa[i]; 616d84be18SBarry Smith aa[i] = 0.; 622bb67c89SBarry Smith } 636d84be18SBarry Smith for (j = kp1; j <= n; ++j) { 642bb67c89SBarry Smith t = work[j]; 656d84be18SBarry Smith ax = &a[j * n + 1]; 6639d66777SBarry Smith ay = &a[kn + 1]; 6726fbe8dcSKarl Rupp for (ll = 0; ll < n; ll++) ay[ll] += t * ax[ll]; 682bb67c89SBarry Smith } 692bb67c89SBarry Smith l = ipvt[k]; 702bb67c89SBarry Smith if (l != k) { 7139d66777SBarry Smith ax = &a[kn + 1]; 726d84be18SBarry Smith ay = &a[l * n + 1]; 736d84be18SBarry Smith for (ll = 0; ll < n; ll++) { 746d84be18SBarry Smith tmp = ax[ll]; 756d84be18SBarry Smith ax[ll] = ay[ll]; 766d84be18SBarry Smith ay[ll] = tmp; 772bb67c89SBarry Smith } 782bb67c89SBarry Smith } 796d84be18SBarry Smith } 803ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 816d84be18SBarry Smith } 82