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: if (zeropivotdetected) *zeropivotdetected = PETSC_FALSE;
24: /* Parameter adjustments */
25: --ipvt;
26: a -= n + 1;
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;
35: /* find l = pivot index */
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) { max = tmp; l = ll+1;}
44: }
45: l += k - 1;
46: ipvt[k] = l;
48: if (a[l + kn] == 0.0) {
49: if (allowzeropivot) {
50: PetscInfo(NULL,"Zero pivot, row %" PetscInt_FMT "\n",k-1);
51: if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE;
52: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %" PetscInt_FMT,k-1);
53: }
55: /* interchange if necessary */
56: if (l != k) {
57: t = a[l + kn];
58: a[l + kn] = a[knp1];
59: a[knp1] = t;
60: }
62: /* compute multipliers */
63: t = -1. / a[knp1];
64: i__2 = n - k;
65: aa = &a[1 + knp1];
66: for (ll=0; ll<i__2; ll++) aa[ll] *= t;
68: /* row elimination with column indexing */
69: ax = aa;
70: for (j = kp1; j <= n; ++j) {
71: jn1 = j*n;
72: t = a[l + jn1];
73: if (l != k) {
74: a[l + jn1] = a[k + jn1];
75: a[k + jn1] = t;
76: }
78: i__3 = n - k;
79: ay = &a[1+k+jn1];
80: for (ll=0; ll<i__3; ll++) ay[ll] += t*ax[ll];
81: }
82: }
83: ipvt[n] = n;
84: if (a[n + n * n] == 0.0) {
85: if (allowzeropivot) {
86: PetscInfo(NULL,"Zero pivot, row %" PetscInt_FMT "\n",n-1);
87: if (zeropivotdetected) *zeropivotdetected = PETSC_TRUE;
88: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %" PetscInt_FMT,n-1);
89: }
90: return 0;
91: }