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++) {
 32:       PetscArrayzero(rtmp+bs2*ajtmp[j],bs2);
 33:     }
 34:     /* load in initial (unfactored row) */
 35:     nz       = ai[r[i]+1] - ai[r[i]];
 36:     ajtmpold = aj + ai[r[i]];
 37:     v        = aa + bs2*ai[r[i]];
 38:     for (j=0; j<nz; j++) {
 39:       PetscArraycpy(rtmp+bs2*ic[ajtmpold[j]],v+bs2*j,bs2);
 40:     }
 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++) {
 58:           PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j);
 59:         }
 60:         PetscLogFlops(2.0*bs*bs2*(nz+1.0)-bs);
 61:       }
 62:       row = *ajtmp++;
 63:     }
 64:     /* finished row so stick it into b->a */
 65:     pv = ba + bs2*bi[i];
 66:     pj = bj + bi[i];
 67:     nz = bi[i+1] - bi[i];
 68:     for (j=0; j<nz; j++) {
 69:       PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);
 70:     }
 71:     diag = diag_offset[i] - bi[i];
 72:     /* invert diagonal block */
 73:     w    = pv + bs2*diag;

 75:     PetscKernel_A_gets_inverse_A(bs,w,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
 76:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
 77:   }

 79:   PetscFree(rtmp);
 80:   PetscFree3(v_work,multiplier,v_pivots);
 81:   ISRestoreIndices(isicol,&ic);
 82:   ISRestoreIndices(isrow,&r);

 84:   C->ops->solve          = MatSolve_SeqBAIJ_N_inplace;
 85:   C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N_inplace;
 86:   C->assembled           = PETSC_TRUE;

 88:   PetscLogFlops(1.333333333333*bs*bs2*b->mbs); /* from inverting diagonal blocks */
 89:   return 0;
 90: }