1 /* 2 This routine was converted by f2c from Linpack source 3 linpack. this version dated 08/14/78 4 cleve moler, university of new mexico, argonne national lab. 5 6 Does an LU factorization with partial pivoting of a dense 7 n by n matrix. 8 9 Used by the sparse factorization routines in 10 src/mat/impls/baij/seq 11 12 */ 13 #include <petsc/private/kernels/blockinvert.h> 14 15 PetscErrorCode PetscLINPACKgefa(MatScalar *a, PetscInt n, PetscInt *ipvt, PetscBool allowzeropivot, PetscBool *zeropivotdetected) 16 { 17 PetscInt i__2, i__3, kp1, nm1, j, k, l, ll, kn, knp1, jn1; 18 MatScalar t, *ax, *ay, *aa; 19 MatReal tmp, max; 20 21 PetscFunctionBegin; 22 if (zeropivotdetected) *zeropivotdetected = PETSC_FALSE; 23 24 /* Parameter adjustments */ 25 --ipvt; 26 a -= n + 1; 27 28 /* Function Body */ 29 nm1 = n - 1; 30 for (k = 1; k <= nm1; ++k) { 31 kp1 = k + 1; 32 kn = k * n; 33 knp1 = k * n + k; 34 35 /* find l = pivot index */ 36 37 i__2 = n - k + 1; 38 aa = &a[knp1]; 39 max = PetscAbsScalar(aa[0]); 40 l = 1; 41 for (ll = 1; ll < i__2; ll++) { 42 tmp = PetscAbsScalar(aa[ll]); 43 if (tmp > max) { 44 max = tmp; 45 l = ll + 1; 46 } 47 } 48 l += k - 1; 49 ipvt[k] = l; 50 51 if (a[l + kn] == 0.0) { 52 PetscCheck(allowzeropivot, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot, row %" PetscInt_FMT, k - 1); 53 PetscCall(PetscInfo(NULL, "Zero pivot, row %" PetscInt_FMT "\n", k - 1)); 54 if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE; 55 } 56 57 /* interchange if necessary */ 58 if (l != k) { 59 t = a[l + kn]; 60 a[l + kn] = a[knp1]; 61 a[knp1] = t; 62 } 63 64 /* compute multipliers */ 65 t = -1. / a[knp1]; 66 i__2 = n - k; 67 aa = &a[1 + knp1]; 68 for (ll = 0; ll < i__2; ll++) aa[ll] *= t; 69 70 /* row elimination with column indexing */ 71 ax = aa; 72 for (j = kp1; j <= n; ++j) { 73 jn1 = j * n; 74 t = a[l + jn1]; 75 if (l != k) { 76 a[l + jn1] = a[k + jn1]; 77 a[k + jn1] = t; 78 } 79 80 i__3 = n - k; 81 ay = &a[1 + k + jn1]; 82 for (ll = 0; ll < i__3; ll++) ay[ll] += t * ax[ll]; 83 } 84 } 85 ipvt[n] = n; 86 if (a[n + n * n] == 0.0) { 87 PetscCheck(allowzeropivot, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot, row %" PetscInt_FMT, n - 1); 88 PetscCall(PetscInfo(NULL, "Zero pivot, row %" PetscInt_FMT "\n", n - 1)); 89 if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE; 90 } 91 PetscFunctionReturn(PETSC_SUCCESS); 92 } 93