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