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 #undef __FUNCT__ 10 #define __FUNCT__ "MatLUFactorNumeric_SeqBAIJ_N_inplace" 11 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N_inplace(Mat C,Mat A,const MatFactorInfo *info) 12 { 13 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data; 14 IS isrow = b->row,isicol = b->icol; 15 PetscErrorCode ierr; 16 const PetscInt *r,*ic; 17 PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j; 18 PetscInt *ajtmpold,*ajtmp,nz,row,bslog,*ai=a->i,*aj=a->j,k,flg; 19 PetscInt *diag_offset=b->diag,diag,bs=A->rmap->bs,bs2 = a->bs2,*pj,*v_pivots; 20 MatScalar *ba = b->a,*aa = a->a,*pv,*v,*rtmp,*multiplier,*v_work,*pc,*w; 21 22 PetscFunctionBegin; 23 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 24 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 25 ierr = PetscMalloc1(bs2*(n+1),&rtmp);CHKERRQ(ierr); 26 ierr = PetscMemzero(rtmp,(bs2*n+1)*sizeof(MatScalar));CHKERRQ(ierr); 27 /* generate work space needed by dense LU factorization */ 28 ierr = PetscMalloc3(bs,&v_work,bs2,&multiplier,bs,&v_pivots);CHKERRQ(ierr); 29 30 /* flops in while loop */ 31 bslog = 2*bs*bs2; 32 33 for (i=0; i<n; i++) { 34 nz = bi[i+1] - bi[i]; 35 ajtmp = bj + bi[i]; 36 for (j=0; j<nz; j++) { 37 ierr = PetscMemzero(rtmp+bs2*ajtmp[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 38 } 39 /* load in initial (unfactored row) */ 40 nz = ai[r[i]+1] - ai[r[i]]; 41 ajtmpold = aj + ai[r[i]]; 42 v = aa + bs2*ai[r[i]]; 43 for (j=0; j<nz; j++) { 44 ierr = PetscMemcpy(rtmp+bs2*ic[ajtmpold[j]],v+bs2*j,bs2*sizeof(MatScalar));CHKERRQ(ierr); 45 } 46 row = *ajtmp++; 47 while (row < i) { 48 pc = rtmp + bs2*row; 49 /* if (*pc) { */ 50 for (flg=0,k=0; k<bs2; k++) { 51 if (pc[k]!=0.0) { 52 flg = 1; 53 break; 54 } 55 } 56 if (flg) { 57 pv = ba + bs2*diag_offset[row]; 58 pj = bj + diag_offset[row] + 1; 59 PetscKernel_A_gets_A_times_B(bs,pc,pv,multiplier); 60 nz = bi[row+1] - diag_offset[row] - 1; 61 pv += bs2; 62 for (j=0; j<nz; j++) { 63 PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); 64 } 65 ierr = PetscLogFlops(bslog*(nz+1.0)-bs);CHKERRQ(ierr); 66 } 67 row = *ajtmp++; 68 } 69 /* finished row so stick it into b->a */ 70 pv = ba + bs2*bi[i]; 71 pj = bj + bi[i]; 72 nz = bi[i+1] - bi[i]; 73 for (j=0; j<nz; j++) { 74 ierr = PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 75 } 76 diag = diag_offset[i] - bi[i]; 77 /* invert diagonal block */ 78 w = pv + bs2*diag; 79 ierr = PetscKernel_A_gets_inverse_A(bs,w,v_pivots,v_work);CHKERRQ(ierr); 80 } 81 82 ierr = PetscFree(rtmp);CHKERRQ(ierr); 83 ierr = PetscFree3(v_work,multiplier,v_pivots);CHKERRQ(ierr); 84 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 85 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 86 87 C->ops->solve = MatSolve_SeqBAIJ_N_inplace; 88 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N_inplace; 89 C->assembled = PETSC_TRUE; 90 91 ierr = PetscLogFlops(1.333333333333*bs*bs2*b->mbs);CHKERRQ(ierr); /* from inverting diagonal blocks */ 92 PetscFunctionReturn(0); 93 } 94