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