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 <petscsys.h> 14 15 PETSC_INTERN 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 if (allowzeropivot) { 53 PetscCall(PetscInfo(NULL, "Zero pivot, row %" PetscInt_FMT "\n", k - 1)); 54 if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE; 55 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot, row %" PetscInt_FMT, k - 1); 56 } 57 58 /* interchange if necessary */ 59 if (l != k) { 60 t = a[l + kn]; 61 a[l + kn] = a[knp1]; 62 a[knp1] = t; 63 } 64 65 /* compute multipliers */ 66 t = -1. / a[knp1]; 67 i__2 = n - k; 68 aa = &a[1 + knp1]; 69 for (ll = 0; ll < i__2; ll++) aa[ll] *= t; 70 71 /* row elimination with column indexing */ 72 ax = aa; 73 for (j = kp1; j <= n; ++j) { 74 jn1 = j * n; 75 t = a[l + jn1]; 76 if (l != k) { 77 a[l + jn1] = a[k + jn1]; 78 a[k + jn1] = t; 79 } 80 81 i__3 = n - k; 82 ay = &a[1 + k + jn1]; 83 for (ll = 0; ll < i__3; ll++) ay[ll] += t * ax[ll]; 84 } 85 } 86 ipvt[n] = n; 87 if (a[n + n * n] == 0.0) { 88 if (allowzeropivot) { 89 PetscCall(PetscInfo(NULL, "Zero pivot, row %" PetscInt_FMT "\n", n - 1)); 90 if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE; 91 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot, row %" PetscInt_FMT, n - 1); 92 } 93 PetscFunctionReturn(PETSC_SUCCESS); 94 } 95