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