1 #ifndef lint 2 static char vcid[] = "$Id: mpiaijpc.c,v 1.21 1996/11/07 15:09:25 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 "src/mat/impls/aij/mpi/mpiaij.h" 11 #include "src/pc/pcimpl.h" 12 #include "src/vec/vecimpl.h" 13 #include "src/pc/impls/bjacobi/bjacobi.h" 14 #include "sles.h" 15 16 typedef struct { 17 Vec x, y; 18 } PC_BJacobi_MPIAIJ; 19 20 int PCDestroy_BJacobi_MPIAIJ(PetscObject obj) 21 { 22 PC pc = (PC) obj; 23 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 24 PC_BJacobi_MPIAIJ *bjac = (PC_BJacobi_MPIAIJ *) jac->data; 25 int ierr; 26 27 ierr = SLESDestroy(jac->sles[0]); CHKERRQ(ierr); 28 PetscFree(jac->sles); 29 ierr = VecDestroy(bjac->x); CHKERRQ(ierr); 30 ierr = VecDestroy(bjac->y); CHKERRQ(ierr); 31 if (jac->l_lens) PetscFree(jac->l_lens); 32 if (jac->g_lens) PetscFree(jac->g_lens); 33 PetscFree(bjac); PetscFree(jac); 34 return 0; 35 } 36 37 38 int PCSetUpOnBlocks_BJacobi_MPIAIJ(PC pc) 39 { 40 int ierr; 41 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 42 PC_BJacobi_MPIAIJ *bjac = (PC_BJacobi_MPIAIJ *) jac->data; 43 44 ierr = SLESSetUp(jac->sles[0],bjac->x,bjac->y); CHKERRQ(ierr); 45 return 0; 46 } 47 48 int PCApply_BJacobi_MPIAIJ(PC pc,Vec x, Vec y) 49 { 50 int ierr,its; 51 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 52 PC_BJacobi_MPIAIJ *bjac = (PC_BJacobi_MPIAIJ *) jac->data; 53 Scalar *x_array,*x_true_array, *y_array,*y_true_array; 54 55 /* 56 The VecPlaceArray() is to avoid having to copy the 57 y vector into the bjac->x vector. The reason for 58 the bjac->x vector is that we need a sequential vector 59 for the sequential solve. 60 */ 61 VecGetArray_Fast(x,x_array); 62 VecGetArray_Fast(y,y_array); 63 VecGetArray_Fast(bjac->x,x_true_array); 64 VecGetArray_Fast(bjac->y,y_true_array); 65 VecPlaceArray_Fast(bjac->x,x_array); 66 VecPlaceArray_Fast(bjac->y,y_array); 67 ierr = SLESSolve(jac->sles[0],bjac->x,bjac->y,&its); 68 VecPlaceArray_Fast(bjac->x,x_true_array); 69 VecPlaceArray_Fast(bjac->y,y_true_array); 70 return 0; 71 } 72 73 int PCSetUp_BJacobi_MPIAIJ(PC pc) 74 { 75 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 76 Mat mat = pc->mat, pmat = pc->pmat; 77 Mat_MPIAIJ *pmatin = (Mat_MPIAIJ *) pmat->data; 78 Mat_MPIAIJ *matin = 0; 79 int ierr, m; 80 SLES sles; 81 Vec x,y; 82 PC_BJacobi_MPIAIJ *bjac; 83 KSP subksp; 84 PC subpc; 85 MatType type; 86 87 if (jac->use_true_local) { 88 MatGetType(pc->mat,&type,PETSC_NULL); 89 if (type != MATMPIAIJ) SETERRQ(1,"PCSetUp_BJacobi_MPIAIJ:Incompatible matrix type."); 90 matin = (Mat_MPIAIJ *) mat->data; 91 } 92 93 /* set default direct solver with no Krylov method */ 94 if (!pc->setupcalled) { 95 char *prefix; 96 ierr = SLESCreate(MPI_COMM_SELF,&sles); CHKERRQ(ierr); 97 PLogObjectParent(pc,sles); 98 ierr = SLESGetKSP(sles,&subksp); CHKERRQ(ierr); 99 ierr = KSPSetType(subksp,KSPPREONLY); CHKERRQ(ierr); 100 ierr = SLESGetPC(sles,&subpc); CHKERRQ(ierr); 101 ierr = PCSetType(subpc,PCILU); CHKERRQ(ierr); 102 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 103 ierr = SLESSetOptionsPrefix(sles,prefix); CHKERRQ(ierr); 104 ierr = SLESAppendOptionsPrefix(sles,"sub_"); CHKERRQ(ierr); 105 ierr = SLESSetFromOptions(sles); CHKERRQ(ierr); 106 /* 107 This is not so good. The only reason we need to generate this vector 108 is so KSP may generate seq vectors for the local solves 109 */ 110 ierr = MatGetSize(pmatin->A,&m,&m); CHKERRQ(ierr); 111 ierr = VecCreateSeq(MPI_COMM_SELF,m,&x); CHKERRQ(ierr); 112 ierr = VecCreateSeq(MPI_COMM_SELF,m,&y); CHKERRQ(ierr); 113 PLogObjectParent(pmat,x); 114 PLogObjectParent(pmat,y); 115 116 pc->destroy = PCDestroy_BJacobi_MPIAIJ; 117 pc->apply = PCApply_BJacobi_MPIAIJ; 118 pc->setuponblocks = PCSetUpOnBlocks_BJacobi_MPIAIJ; 119 120 bjac = (PC_BJacobi_MPIAIJ *) PetscMalloc(sizeof(PC_BJacobi_MPIAIJ));CHKPTRQ(bjac); 121 PLogObjectMemory(pc,sizeof(PC_BJacobi_MPIAIJ)); 122 bjac->x = x; 123 bjac->y = y; 124 125 jac->sles = (SLES*) PetscMalloc( sizeof(SLES) ); CHKPTRQ(jac->sles); 126 jac->sles[0] = sles; 127 jac->data = (void *) bjac; 128 } 129 else { 130 sles = jac->sles[0]; 131 bjac = (PC_BJacobi_MPIAIJ *)jac->data; 132 } 133 /* if (jac->l_true[0] == USE_TRUE_MATRIX) { 134 ierr = SLESSetOperators(sles,matin->A,matin->A,pc->flag); 135 } 136 else */ if (jac->use_true_local) 137 ierr = SLESSetOperators(sles,matin->A,pmatin->A,pc->flag); 138 else 139 ierr = SLESSetOperators(sles,pmatin->A,pmatin->A,pc->flag); 140 CHKERRQ(ierr); 141 return 0; 142 } 143 144 145