1 /*$Id: baijfact4.c,v 1.2 2001/01/15 21:45:50 bsmith Exp balay $*/ 2 /* 3 Factorization code for BAIJ format. 4 */ 5 #include "src/mat/impls/baij/seq/baij.h" 6 #include "src/vec/vecimpl.h" 7 #include "src/inline/ilu.h" 8 9 /* ----------------------------------------------------------- */ 10 #undef __FUNC__ 11 #define __FUNC__ "MatLUFactorNumeric_SeqBAIJ_N" 12 int MatLUFactorNumeric_SeqBAIJ_N(Mat A,Mat *B) 13 { 14 Mat C = *B; 15 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ *)C->data; 16 IS isrow = b->row,isicol = b->icol; 17 int *r,*ic,ierr,i,j,n = a->mbs,*bi = b->i,*bj = b->j; 18 int *ajtmpold,*ajtmp,nz,row,bslog,*ai=a->i,*aj=a->j,k,flg; 19 int *diag_offset=b->diag,diag,bs=a->bs,bs2 = a->bs2,*v_pivots,*pj; 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 = PetscMalloc(bs2*(n+1)*sizeof(MatScalar),&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 = PetscMalloc(bs*sizeof(int) + (bs+bs2)*sizeof(MatScalar),&v_work);CHKERRQ(ierr); 29 multiplier = v_work + bs; 30 v_pivots = (int*)(multiplier + bs2); 31 32 /* flops in while loop */ 33 bslog = 2*bs*bs2; 34 35 for (i=0; i<n; i++) { 36 nz = bi[i+1] - bi[i]; 37 ajtmp = bj + bi[i]; 38 for (j=0; j<nz; j++) { 39 ierr = PetscMemzero(rtmp+bs2*ajtmp[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 40 } 41 /* load in initial (unfactored row) */ 42 nz = ai[r[i]+1] - ai[r[i]]; 43 ajtmpold = aj + ai[r[i]]; 44 v = aa + bs2*ai[r[i]]; 45 for (j=0; j<nz; j++) { 46 ierr = PetscMemcpy(rtmp+bs2*ic[ajtmpold[j]],v+bs2*j,bs2*sizeof(MatScalar));CHKERRQ(ierr); 47 } 48 row = *ajtmp++; 49 while (row < i) { 50 pc = rtmp + bs2*row; 51 /* if (*pc) { */ 52 for (flg=0,k=0; k<bs2; k++) { if (pc[k]!=0.0) { flg = 1; break; }} 53 if (flg) { 54 pv = ba + bs2*diag_offset[row]; 55 pj = bj + diag_offset[row] + 1; 56 Kernel_A_gets_A_times_B(bs,pc,pv,multiplier); 57 nz = bi[row+1] - diag_offset[row] - 1; 58 pv += bs2; 59 for (j=0; j<nz; j++) { 60 Kernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); 61 } 62 PetscLogFlops(bslog*(nz+1)-bs); 63 } 64 row = *ajtmp++; 65 } 66 /* finished row so stick it into b->a */ 67 pv = ba + bs2*bi[i]; 68 pj = bj + bi[i]; 69 nz = bi[i+1] - bi[i]; 70 for (j=0; j<nz; j++) { 71 ierr = PetscMemcpy(pv+bs2*j,rtmp+bs2*pj[j],bs2*sizeof(MatScalar));CHKERRQ(ierr); 72 } 73 diag = diag_offset[i] - bi[i]; 74 /* invert diagonal block */ 75 w = pv + bs2*diag; 76 Kernel_A_gets_inverse_A(bs,w,v_pivots,v_work); 77 } 78 79 ierr = PetscFree(rtmp);CHKERRQ(ierr); 80 ierr = PetscFree(v_work);CHKERRQ(ierr); 81 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 82 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 83 C->factor = FACTOR_LU; 84 C->assembled = PETSC_TRUE; 85 PetscLogFlops(1.3333*bs*bs2*b->mbs); /* from inverting diagonal blocks */ 86 PetscFunctionReturn(0); 87 } 88