1 /*
2 Factorization code for BAIJ format.
3 */
4 #include <../src/mat/impls/baij/seq/baij.h>
5 #include <petsc/private/kernels/blockinvert.h>
6
MatILUFactorNumeric_SeqBAIJ_N_inplace(Mat C,Mat A,const MatFactorInfo * info)7 PetscErrorCode MatILUFactorNumeric_SeqBAIJ_N_inplace(Mat C, Mat A, const MatFactorInfo *info)
8 {
9 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
10 IS isrow = b->row, isicol = b->icol;
11 const PetscInt *r, *ic;
12 PetscInt i, j, n = a->mbs, *bi = b->i, *bj = b->j;
13 PetscInt *ajtmpold, *ajtmp, nz, row, *ai = a->i, *aj = a->j, k, flg;
14 const PetscInt *diag_offset;
15 PetscInt diag, bs = A->rmap->bs, bs2 = a->bs2, *pj, *v_pivots;
16 MatScalar *ba = b->a, *aa = a->a, *pv, *v, *rtmp, *multiplier, *v_work, *pc, *w;
17 PetscBool allowzeropivot, zeropivotdetected;
18
19 PetscFunctionBegin;
20 /* Since A is C and C is labeled as a factored matrix we need to lie to MatGetDiagonalMarkers_SeqBAIJ() to get it to compute the diagonals */
21 A->factortype = MAT_FACTOR_NONE;
22 PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, &diag_offset, NULL));
23 A->factortype = MAT_FACTOR_ILU;
24 PetscCall(ISGetIndices(isrow, &r));
25 PetscCall(ISGetIndices(isicol, &ic));
26 allowzeropivot = PetscNot(A->erroriffailure);
27
28 PetscCall(PetscCalloc1(bs2 * (n + 1), &rtmp));
29 /* generate work space needed by dense LU factorization */
30 PetscCall(PetscMalloc3(bs, &v_work, bs2, &multiplier, bs, &v_pivots));
31
32 for (i = 0; i < n; i++) {
33 nz = bi[i + 1] - bi[i];
34 ajtmp = bj + bi[i];
35 for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * ajtmp[j], bs2));
36 /* load in initial (unfactored row) */
37 nz = ai[r[i] + 1] - ai[r[i]];
38 ajtmpold = aj + ai[r[i]];
39 v = aa + bs2 * ai[r[i]];
40 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmpold[j]], v + bs2 * j, bs2));
41 row = *ajtmp++;
42 while (row < i) {
43 pc = rtmp + bs2 * row;
44 /* if (*pc) { */
45 for (flg = 0, k = 0; k < bs2; k++) {
46 if (pc[k] != 0.0) {
47 flg = 1;
48 break;
49 }
50 }
51 if (flg) {
52 pv = ba + bs2 * diag_offset[row];
53 pj = bj + diag_offset[row] + 1;
54 PetscKernel_A_gets_A_times_B(bs, pc, pv, multiplier);
55 nz = bi[row + 1] - diag_offset[row] - 1;
56 pv += bs2;
57 for (j = 0; j < nz; j++) PetscKernel_A_gets_A_minus_B_times_C(bs, rtmp + bs2 * pj[j], pc, pv + bs2 * j);
58 PetscCall(PetscLogFlops(2.0 * bs * bs2 * (nz + 1.0) - bs));
59 }
60 row = *ajtmp++;
61 }
62 /* finished row so stick it into b->a */
63 pv = ba + bs2 * bi[i];
64 pj = bj + bi[i];
65 nz = bi[i + 1] - bi[i];
66 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
67 diag = diag_offset[i] - bi[i];
68 /* invert diagonal block */
69 w = pv + bs2 * diag;
70
71 PetscCall(PetscKernel_A_gets_inverse_A(bs, w, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
72 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
73 }
74
75 PetscCall(PetscFree(rtmp));
76 PetscCall(PetscFree3(v_work, multiplier, v_pivots));
77 PetscCall(ISRestoreIndices(isicol, &ic));
78 PetscCall(ISRestoreIndices(isrow, &r));
79
80 C->ops->solve = MatSolve_SeqBAIJ_N_inplace;
81 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N_inplace;
82 C->assembled = PETSC_TRUE;
83
84 PetscCall(PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */
85 PetscFunctionReturn(PETSC_SUCCESS);
86 }
87