1*6699817fSBarry Smith #ifndef lint 2*6699817fSBarry Smith static char vcid[] = "$Id: bjmpiaij.c,v 1.20 1995/10/01 21:52:09 bsmith Exp $"; 3*6699817fSBarry Smith #endif 4*6699817fSBarry Smith /* 5*6699817fSBarry Smith Defines a block Jacobi preconditioner for the MPIAIJ format. 6*6699817fSBarry Smith At the moment only supports a single block per processor. 7*6699817fSBarry Smith This file knows about storage formats for MPIAIJ matrices. 8*6699817fSBarry Smith This code is nearly identical to that for the MPIROW format; 9*6699817fSBarry Smith */ 10*6699817fSBarry Smith #include "src/mat/impls/aij/mpi/mpiaij.h" 11*6699817fSBarry Smith #include "pcimpl.h" 12*6699817fSBarry Smith #include "bjacobi.h" 13*6699817fSBarry Smith #include "sles.h" 14*6699817fSBarry Smith 15*6699817fSBarry Smith typedef struct { 16*6699817fSBarry Smith Vec x; 17*6699817fSBarry Smith } PC_BJacobiMPIAIJ; 18*6699817fSBarry Smith 19*6699817fSBarry Smith int PCDestroy_BJacobiMPIAIJ(PetscObject obj) 20*6699817fSBarry Smith { 21*6699817fSBarry Smith PC pc = (PC) obj; 22*6699817fSBarry Smith PC_BJacobi *jac = (PC_BJacobi *) pc->data; 23*6699817fSBarry Smith PC_BJacobiMPIAIJ *bjac = (PC_BJacobiMPIAIJ *) jac->data; 24*6699817fSBarry Smith int i,ierr; 25*6699817fSBarry Smith 26*6699817fSBarry Smith for ( i=0; i<jac->n_local; i++ ) { 27*6699817fSBarry Smith ierr = SLESDestroy(jac->sles[i]); CHKERRQ(ierr); 28*6699817fSBarry Smith } 29*6699817fSBarry Smith PETSCFREE(jac->sles); 30*6699817fSBarry Smith ierr = VecDestroy(bjac->x); CHKERRQ(ierr); 31*6699817fSBarry Smith PETSCFREE(bjac); PETSCFREE(jac); 32*6699817fSBarry Smith return 0; 33*6699817fSBarry Smith } 34*6699817fSBarry Smith 35*6699817fSBarry Smith int PCApply_BJacobiMPIAIJ(PC pc,Vec x, Vec y) 36*6699817fSBarry Smith { 37*6699817fSBarry Smith int ierr,its; 38*6699817fSBarry Smith PC_BJacobi *jac = (PC_BJacobi *) pc->data; 39*6699817fSBarry Smith PC_BJacobiMPIAIJ *bjac = (PC_BJacobiMPIAIJ *) jac->data; 40*6699817fSBarry Smith 41*6699817fSBarry Smith ierr = SLESSolve(jac->sles[0],x,bjac->x,&its); CHKERRQ(ierr); 42*6699817fSBarry Smith ierr = VecCopy(bjac->x,y); CHKERRQ(ierr); 43*6699817fSBarry Smith return 0; 44*6699817fSBarry Smith } 45*6699817fSBarry Smith 46*6699817fSBarry Smith int PCSetUp_BJacobiMPIAIJ(PC pc) 47*6699817fSBarry Smith { 48*6699817fSBarry Smith PC_BJacobi *jac = (PC_BJacobi *) pc->data; 49*6699817fSBarry Smith Mat mat = pc->mat, pmat = pc->pmat; 50*6699817fSBarry Smith Mat_MPIAIJ *pmatin = (Mat_MPIAIJ *) pmat->data; 51*6699817fSBarry Smith Mat_MPIAIJ *matin = 0; 52*6699817fSBarry Smith int ierr, numtids = pmatin->numtids,m; 53*6699817fSBarry Smith SLES sles; 54*6699817fSBarry Smith Vec x; 55*6699817fSBarry Smith PC_BJacobiMPIAIJ *bjac; 56*6699817fSBarry Smith KSP subksp; 57*6699817fSBarry Smith PC subpc; 58*6699817fSBarry Smith MatType type; 59*6699817fSBarry Smith 60*6699817fSBarry Smith if (jac->use_true_local) { 61*6699817fSBarry Smith MatGetType(pc->mat,&type); 62*6699817fSBarry Smith if (type != MATMPIAIJ) SETERRQ(1,"PCSetUp_BJacobiMPIAIJ:Incompatible matrix type."); 63*6699817fSBarry Smith matin = (Mat_MPIAIJ *) mat->data; 64*6699817fSBarry Smith } 65*6699817fSBarry Smith 66*6699817fSBarry Smith if (numtids != jac->n) { 67*6699817fSBarry Smith SETERRQ(1,"PCSetUp_BJacobiMPIAIJ:Can do only 1 block per processor\n"); 68*6699817fSBarry Smith } 69*6699817fSBarry Smith 70*6699817fSBarry Smith /* set default direct solver with no Krylov method */ 71*6699817fSBarry Smith if (!pc->setupcalled) { 72*6699817fSBarry Smith ierr = SLESCreate(MPI_COMM_SELF,&sles); CHKERRQ(ierr); 73*6699817fSBarry Smith PLogObjectParent(pc,sles); 74*6699817fSBarry Smith ierr = SLESGetKSP(sles,&subksp); CHKERRQ(ierr); 75*6699817fSBarry Smith ierr = KSPSetMethod(subksp,KSPPREONLY); CHKERRQ(ierr); 76*6699817fSBarry Smith ierr = SLESGetPC(sles,&subpc); CHKERRQ(ierr); 77*6699817fSBarry Smith ierr = PCSetMethod(subpc,PCLU); CHKERRQ(ierr); 78*6699817fSBarry Smith ierr = SLESSetOptionsPrefix(sles,"-sub_"); CHKERRQ(ierr); 79*6699817fSBarry Smith ierr = SLESSetFromOptions(sles); CHKERRQ(ierr); 80*6699817fSBarry Smith /* 81*6699817fSBarry Smith This is not so good. The only reason we need to generate this vector 82*6699817fSBarry Smith is so KSP may generate seq vectors for the local solves 83*6699817fSBarry Smith */ 84*6699817fSBarry Smith ierr = MatGetSize(pmatin->A,&m,&m); CHKERRQ(ierr); 85*6699817fSBarry Smith ierr = VecCreateSeq(MPI_COMM_SELF,m,&x); CHKERRQ(ierr); 86*6699817fSBarry Smith PLogObjectParent(pmat,x); 87*6699817fSBarry Smith 88*6699817fSBarry Smith pc->destroy = PCDestroy_BJacobiMPIAIJ; 89*6699817fSBarry Smith pc->apply = PCApply_BJacobiMPIAIJ; 90*6699817fSBarry Smith 91*6699817fSBarry Smith bjac = (PC_BJacobiMPIAIJ *) PETSCMALLOC(sizeof(PC_BJacobiMPIAIJ)); 92*6699817fSBarry Smith CHKPTRQ(bjac); 93*6699817fSBarry Smith PLogObjectMemory(pc,sizeof(PC_BJacobiMPIAIJ)); 94*6699817fSBarry Smith bjac->x = x; 95*6699817fSBarry Smith 96*6699817fSBarry Smith jac->sles = (SLES*) PETSCMALLOC( sizeof(SLES) ); CHKPTRQ(jac->sles); 97*6699817fSBarry Smith jac->sles[0] = sles; 98*6699817fSBarry Smith jac->data = (void *) bjac; 99*6699817fSBarry Smith } 100*6699817fSBarry Smith else { 101*6699817fSBarry Smith sles = jac->sles[0]; 102*6699817fSBarry Smith } 103*6699817fSBarry Smith if (jac->use_true_local) 104*6699817fSBarry Smith ierr = SLESSetOperators(sles,matin->A,pmatin->A,pc->flag); 105*6699817fSBarry Smith else 106*6699817fSBarry Smith ierr = SLESSetOperators(sles,pmatin->A,pmatin->A,pc->flag); 107*6699817fSBarry Smith CHKERRQ(ierr); 108*6699817fSBarry Smith 109*6699817fSBarry Smith return 0; 110*6699817fSBarry Smith } 111*6699817fSBarry Smith 112