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