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 { 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; 19 20 PetscFunctionBegin; 21 PetscCall(ISGetIndices(isrow,&r)); 22 PetscCall(ISGetIndices(isicol,&ic)); 23 allowzeropivot = PetscNot(A->erroriffailure); 24 25 PetscCall(PetscCalloc1(bs2*(n+1),&rtmp)); 26 /* generate work space needed by dense LU factorization */ 27 PetscCall(PetscMalloc3(bs,&v_work,bs2,&multiplier,bs,&v_pivots)); 28 29 for (i=0; i<n; i++) { 30 nz = bi[i+1] - bi[i]; 31 ajtmp = bj + bi[i]; 32 for (j=0; j<nz; j++) { 33 PetscCall(PetscArrayzero(rtmp+bs2*ajtmp[j],bs2)); 34 } 35 /* load in initial (unfactored row) */ 36 nz = ai[r[i]+1] - ai[r[i]]; 37 ajtmpold = aj + ai[r[i]]; 38 v = aa + bs2*ai[r[i]]; 39 for (j=0; j<nz; j++) { 40 PetscCall(PetscArraycpy(rtmp+bs2*ic[ajtmpold[j]],v+bs2*j,bs2)); 41 } 42 row = *ajtmp++; 43 while (row < i) { 44 pc = rtmp + bs2*row; 45 /* if (*pc) { */ 46 for (flg=0,k=0; k<bs2; k++) { 47 if (pc[k]!=0.0) { 48 flg = 1; 49 break; 50 } 51 } 52 if (flg) { 53 pv = ba + bs2*diag_offset[row]; 54 pj = bj + diag_offset[row] + 1; 55 PetscKernel_A_gets_A_times_B(bs,pc,pv,multiplier); 56 nz = bi[row+1] - diag_offset[row] - 1; 57 pv += bs2; 58 for (j=0; j<nz; j++) { 59 PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); 60 } 61 PetscCall(PetscLogFlops(2.0*bs*bs2*(nz+1.0)-bs)); 62 } 63 row = *ajtmp++; 64 } 65 /* finished row so stick it into b->a */ 66 pv = ba + bs2*bi[i]; 67 pj = bj + bi[i]; 68 nz = bi[i+1] - bi[i]; 69 for (j=0; j<nz; j++) { 70 PetscCall(PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2)); 71 } 72 diag = diag_offset[i] - bi[i]; 73 /* invert diagonal block */ 74 w = pv + bs2*diag; 75 76 PetscCall(PetscKernel_A_gets_inverse_A(bs,w,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 77 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 78 } 79 80 PetscCall(PetscFree(rtmp)); 81 PetscCall(PetscFree3(v_work,multiplier,v_pivots)); 82 PetscCall(ISRestoreIndices(isicol,&ic)); 83 PetscCall(ISRestoreIndices(isrow,&r)); 84 85 C->ops->solve = MatSolve_SeqBAIJ_N_inplace; 86 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N_inplace; 87 C->assembled = PETSC_TRUE; 88 89 PetscCall(PetscLogFlops(1.333333333333*bs*bs2*b->mbs)); /* from inverting diagonal blocks */ 90 PetscFunctionReturn(0); 91 } 92