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: }