1 #ifndef lint 2 static char vcid[] = "$Id: mpiaijpc.c,v 1.7 1995/12/03 02:42:38 bsmith Exp bsmith $"; 3 #endif 4 /* 5 Defines a block Jacobi preconditioner for the MPIAIJ format. 6 Handles special case of single block per processor. 7 This file knows about storage formats for MPIAIJ matrices. 8 The general case is handled in aijpc.c 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 ierr; 25 26 ierr = SLESDestroy(jac->sles[0]); CHKERRQ(ierr); 27 PetscFree(jac->sles); 28 ierr = VecDestroy(bjac->x); CHKERRQ(ierr); 29 PetscFree(bjac); PetscFree(jac); 30 return 0; 31 } 32 33 int PCApply_BJacobiMPIAIJ(PC pc,Vec x, Vec y) 34 { 35 int ierr,its; 36 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 37 PC_BJacobiMPIAIJ *bjac = (PC_BJacobiMPIAIJ *) jac->data; 38 Scalar *array,*true_array; 39 40 /* 41 The VecPlaceArray() is to avoid having to copy the 42 y vector into the bjac->x vector. The reason for 43 the bjac->x vector is that we need a sequential vector 44 for the sequential solve. 45 */ 46 ierr = VecGetArray(y,&array); CHKERRQ(ierr); 47 ierr = VecGetArray(bjac->x,&true_array); CHKERRQ(ierr); 48 ierr = VecPlaceArray(bjac->x,array); CHKERRQ(ierr); 49 ierr = SLESSolve(jac->sles[0],x,bjac->x,&its); CHKERRQ(ierr); 50 ierr = VecPlaceArray(bjac->x,true_array); CHKERRQ(ierr); 51 52 return 0; 53 } 54 55 int PCSetUp_BJacobiMPIAIJ(PC pc) 56 { 57 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 58 Mat mat = pc->mat, pmat = pc->pmat; 59 Mat_MPIAIJ *pmatin = (Mat_MPIAIJ *) pmat->data; 60 Mat_MPIAIJ *matin = 0; 61 int ierr, m; 62 SLES sles; 63 Vec x; 64 PC_BJacobiMPIAIJ *bjac; 65 KSP subksp; 66 PC subpc; 67 MatType type; 68 69 if (jac->use_true_local) { 70 MatGetType(pc->mat,&type); 71 if (type != MATMPIAIJ) SETERRQ(1,"PCSetUp_BJacobiMPIAIJ:Incompatible matrix type."); 72 matin = (Mat_MPIAIJ *) mat->data; 73 } 74 75 /* set default direct solver with no Krylov method */ 76 if (!pc->setupcalled) { 77 ierr = SLESCreate(MPI_COMM_SELF,&sles); CHKERRQ(ierr); 78 PLogObjectParent(pc,sles); 79 ierr = SLESGetKSP(sles,&subksp); CHKERRQ(ierr); 80 ierr = KSPSetMethod(subksp,KSPPREONLY); CHKERRQ(ierr); 81 ierr = SLESGetPC(sles,&subpc); CHKERRQ(ierr); 82 ierr = PCSetMethod(subpc,PCLU); CHKERRQ(ierr); 83 ierr = SLESSetOptionsPrefix(sles,"-sub_"); CHKERRQ(ierr); 84 ierr = SLESSetFromOptions(sles); CHKERRQ(ierr); 85 /* 86 This is not so good. The only reason we need to generate this vector 87 is so KSP may generate seq vectors for the local solves 88 */ 89 ierr = MatGetSize(pmatin->A,&m,&m); CHKERRQ(ierr); 90 ierr = VecCreateSeq(MPI_COMM_SELF,m,&x); CHKERRQ(ierr); 91 PLogObjectParent(pmat,x); 92 93 pc->destroy = PCDestroy_BJacobiMPIAIJ; 94 pc->apply = PCApply_BJacobiMPIAIJ; 95 96 bjac = (PC_BJacobiMPIAIJ *) PetscMalloc(sizeof(PC_BJacobiMPIAIJ));CHKPTRQ(bjac); 97 PLogObjectMemory(pc,sizeof(PC_BJacobiMPIAIJ)); 98 bjac->x = x; 99 100 jac->sles = (SLES*) PetscMalloc( sizeof(SLES) ); CHKPTRQ(jac->sles); 101 jac->sles[0] = sles; 102 jac->data = (void *) bjac; 103 } 104 else { 105 sles = jac->sles[0]; 106 } 107 if (jac->l_true[0] == USE_TRUE_MATRIX) { 108 ierr = SLESSetOperators(sles,matin->A,matin->A,pc->flag); 109 } 110 else if (jac->use_true_local) 111 ierr = SLESSetOperators(sles,matin->A,pmatin->A,pc->flag); 112 else 113 ierr = SLESSetOperators(sles,pmatin->A,pmatin->A,pc->flag); 114 CHKERRQ(ierr); 115 116 return 0; 117 } 118 119 120