1 2 /* 3 Factorization code for BAIJ format. 4 */ 5 #include <../src/mat/impls/baij/seq/baij.h> 6 #include <petsc/private/kernels/blockinvert.h> 7 8 /* ----------------------------------------------------------- */ 9 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N_inplace(Mat C, Mat A, const MatFactorInfo *info) { 10 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data; 11 IS isrow = b->row, isicol = b->icol; 12 const PetscInt *r, *ic; 13 PetscInt i, j, n = a->mbs, *bi = b->i, *bj = b->j; 14 PetscInt *ajtmpold, *ajtmp, nz, row, *ai = a->i, *aj = a->j, k, flg; 15 PetscInt *diag_offset = b->diag, 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 PetscCall(ISGetIndices(isrow, &r)); 21 PetscCall(ISGetIndices(isicol, &ic)); 22 allowzeropivot = PetscNot(A->erroriffailure); 23 24 PetscCall(PetscCalloc1(bs2 * (n + 1), &rtmp)); 25 /* generate work space needed by dense LU factorization */ 26 PetscCall(PetscMalloc3(bs, &v_work, bs2, &multiplier, bs, &v_pivots)); 27 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++) { PetscCall(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++) { PetscCall(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 PetscCall(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++) { PetscCall(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; 66 67 PetscCall(PetscKernel_A_gets_inverse_A(bs, w, v_pivots, v_work, allowzeropivot, &zeropivotdetected)); 68 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 69 } 70 71 PetscCall(PetscFree(rtmp)); 72 PetscCall(PetscFree3(v_work, multiplier, v_pivots)); 73 PetscCall(ISRestoreIndices(isicol, &ic)); 74 PetscCall(ISRestoreIndices(isrow, &r)); 75 76 C->ops->solve = MatSolve_SeqBAIJ_N_inplace; 77 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N_inplace; 78 C->assembled = PETSC_TRUE; 79 80 PetscCall(PetscLogFlops(1.333333333333 * bs * bs2 * b->mbs)); /* from inverting diagonal blocks */ 81 PetscFunctionReturn(0); 82 } 83