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