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