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