Actual source code: taolinearsolver_petsc.c

  1: #include "tao_general.h"

  3: #ifdef TAO_USE_PETSC

  5: #include "taolinearsolver_petsc.h"

  7: #include "src/ksp/ksp/kspimpl.h"
  8: #include "src/ksp/ksp/impls/qcg/qcg.h"

 10: #include "../matrix/taomat_petsc.h"
 11:  #include taovec.h
 12: // #include "../vector/taovec_petsc.h"

 14: int KSPDuplicate(KSP ksp1, KSP* kspout);

 17: TaoLinearSolverPetsc::TaoLinearSolverPetsc(KSP S):TaoLinearSolver(){
 18:   linear_its=0;
 19:   rrmap=0;rcmap=0;
 20:   ksp=0;
 21:   this->pkspviewer=PETSC_VIEWER_STDOUT_WORLD;
 22:   this->SetKSP(S);
 23:   return;
 24: }

 26: TaoLinearSolverPetsc::~TaoLinearSolverPetsc(){
 27:   if (ksp){
 28:     KSPDestroy(ksp);
 29:   }
 30:   if (this->rrmap && this->rcmap){
 31:     PetscObjectDereference((PetscObject)this->rrmap);
 32:     PetscObjectDereference((PetscObject)this->rcmap);
 33:   };
 34:   return;
 35: }

 39: /*@C
 40:    TaoWrapKSP - Create a new TaoLinearSolver object using PETSc KSP.

 42:    Input Parameter:
 43: .  S -  a KSP

 45:    Output Parameter:
 46: .  SS - new TaoMat

 48:    Note:  
 49:    A TaoLinearSolverPetsc is an object with the methods of an abstract
 50:    TaoLinearSolver object.  A TaoLinearSolverPetsc contains an implementation 
 51:    of the TaoLinearSolver methods.  Routines using these vectors should 
 52:    declare a pointer to a TaoLinearSolver, assign this pointer to the address 
 53:    of a TaoLinearSolver object, use the pointer to invoke methods on the 
 54:    object, and use this pointer as an argument when calling other routines.  
 55:    This usage is different from the usage of a PETSc KSP.  In PETSc, 
 56:    applications will typically declare a KSP, and pass it as an argument into 
 57:    routines.  That is, applications will typically declare a pointer to a 
 58:    TaoLinearSolver and use the pointer, or declare a KSP and use it directly.

 60:    Level: developer

 62: @*/
 63: int TaoWrapKSP( KSP S, TaoLinearSolverPetsc ** SS){
 66:   *SS = new  TaoLinearSolverPetsc(S);
 67:   return(0);
 68: }

 72: /*@C
 73:    TaoLinearSolverGetKSP - If the TaoLinearSolver is of the TaoLinearSolverPetsc type, it gets the underlying PETSc KSP 

 75:    Input Parameter:
 76: +  TS - the TaoLinearSolver
 77: -  S -  the address of KSP

 79:    Output Parameter:
 80: .  S -  address of the underlying PETSc KSP object


 83:    Note:  
 84:    This routine does not create a KSP.  It sets a pointer
 85:    to the location of an existing KSP.

 87:    Level: advanced

 89: .seealso TaoVecGetPetscVec(), TaoMatGetPetscMat(), TaoGetKSP()
 90: @*/
 91: int TaoLinearSolverGetKSP( TaoLinearSolver *TS, KSP * S){
 93:   if (TS){
 94:     *S=((TaoLinearSolverPetsc *)TS)->GetKSP();
 96:   } else {
 97:     *S=0;
 98:   }
 99:   return(0);
100: }

104: int TaoLinearSolverPetsc::SetKSP(KSP ksp2){
105:   int info;
107:   if (ksp2){
109:     PetscObjectReference((PetscObject)ksp2);
110:   }
111:   if (this->ksp){
112:     info=KSPDestroy(ksp); CHKERRQ(info);
113:   }
114:   this->ksp=ksp2;
115:   this->LinearSolverObject=(void*)ksp2;
116:   return(0);
117: }
120: int TaoLinearSolverPetsc::GetKSP(KSP *ksp2){
122:   if (ksp2){
123:     *ksp2=this->ksp;
124:   }
125:   return(0);
126: }

130: int TaoLinearSolverPetsc::PreSolve(TaoMat* m1){
131:   int info;
132:   TaoMatPetsc *MM=(TaoMatPetsc*)m1;
133:   Mat mm, mm_pre;
134:   MatStructure preflag;
135:   PetscMap rmap1,cmap1;
136:   PetscTruth samemap=PETSC_TRUE;
137:   KSP tksp, ksp1=this->ksp;

140:   info = MM->GetMatrix(&mm,&mm_pre,&preflag);CHKERRQ(info);
141:   info = MatGetPetscMaps(mm,&rmap1,&cmap1);CHKERRQ(info);

143:   if (this->rrmap && this->rcmap){
144:     info=PetscMapSame(rmap1,this->rrmap, &samemap);CHKERRQ(info);
145:     if (samemap){
146:       info=PetscMapSame(cmap1,this->rcmap, &samemap);CHKERRQ(info);
147:     }
148:     PetscObjectDereference((PetscObject)this->rrmap);
149:     PetscObjectDereference((PetscObject)this->rcmap);
150:   }

152:   if (samemap==PETSC_FALSE){
153:     info = KSPDuplicate(ksp1,&tksp); CHKERRQ(info);
154:     info = this->SetKSP(tksp); CHKERRQ(info);
155:     info = KSPDestroy(tksp);  CHKERRQ(info);
156:   }

158:   this->rrmap=rmap1; this->rcmap=cmap1;

160:   PetscObjectReference((PetscObject)rmap1);
161:   PetscObjectReference((PetscObject)cmap1);

163:   info = KSPSetOperators(this->ksp,mm,mm_pre,preflag);CHKERRQ(info);
164:   return(0);
165: }

169: int TaoLinearSolverPetsc::SetTrustRadius(double rad){

171:   KSPType ktype;
172:   KSP_QCG *qcgP;
173:   int info;
174:   PetscTruth flg;


178:   info=KSPGetType(ksp,&ktype); CHKERRQ(info);
179:   info = PetscStrcmp((char *)ktype,KSPQCG,&flg); CHKERRQ(info);
180:   if (flg==PETSC_TRUE){
181:     qcgP = (KSP_QCG*) ksp->data;
182:     qcgP->delta=rad;
183:   }
184:   return(0);
185: }

189: int TaoLinearSolverPetsc::Solve(TaoVec* v1, TaoVec* v2, TaoTruth *flag){
190:   Vec vv=(Vec) v1->VecObject;
191:   Vec ww=(Vec) v2->VecObject;
192:   int info,its;

195:   info = KSPSolve(ksp,vv,ww); CHKERRQ(info);
196:   info = KSPGetIterationNumber(ksp, &its); CHKERRQ(info);

198:   this->linear_its=PetscMax(its,-its);
199:   if (its>0) *flag=TAO_TRUE;
200:   else *flag=TAO_FALSE;
201:   return(0);
202: }

206: int  TaoLinearSolverPetsc::MinQuadraticTrustRegion(TaoVec*tb,TaoVec*tx,double r, TaoTruth*flg){
207:   int info;

210:   info=KSPSetType(ksp,KSPQCG); CHKERRQ(info);
211:   info = this->SetTrustRadius(r);CHKERRQ(info);
212:   info = this->Solve(tb,tx,flg);CHKERRQ(info);
213:   info = tx->Negate(); CHKERRQ(info);
214:   return(0);
215: }

219: int TaoLinearSolverPetsc::GetNumberIterations(int * iters){
221:   *iters=linear_its;
222:   return(0);
223: }

227: int TaoLinearSolverPetsc::SetTolerances(double rtol, double atol, double dtol, int maxits){
228:   int info;

231:   info = KSPSetTolerances(ksp,rtol,atol,dtol,maxits);CHKERRQ(info);
232:   return(0);
233: }

237: int TaoLinearSolverPetsc::View(){
238:   int info;
240:   info=KSPView(this->ksp,this->pkspviewer);CHKERRQ(info);
241:   return(0);
242: }

246: int TaoLinearSolverPetsc::SetOptions(){
247:   int info;
249:   info=KSPSetFromOptions(ksp); CHKERRQ(info);
250:   return(0);
251: }

255: int KSPDuplicate(KSP ksp1, KSP* kspout){
256:   const char *prefix;
257:   PetscReal rtol,atol,dtol;
258:   int maxits,info;
259:   KSPType ktype;
260:   KSP ksp2;
261:   PCType ptype;
262:   PC  pc1,pc2;
263:   PetscTruth flg;
264:   MPI_Comm comm;
266: 
267:   info = PetscObjectGetComm((PetscObject)ksp1,&comm); CHKERRQ(info);
268:   info=KSPCreate(comm,&ksp2); CHKERRQ(info);

270:   info=KSPGetOptionsPrefix( ksp1,&prefix); CHKERRQ(info);
271:   info=KSPSetOptionsPrefix( ksp2, prefix); CHKERRQ(info);

273:   info=KSPGetType(ksp1,&ktype); CHKERRQ(info);
274:   if (ktype){
275:     info=KSPSetType(ksp2, ktype); CHKERRQ(info);
276:   }
277:   info=KSPGetTolerances(ksp1,&rtol,&atol,&dtol,&maxits); CHKERRQ(info);
278:   info=KSPSetTolerances(ksp2, rtol, atol, dtol, PETSC_DEFAULT); CHKERRQ(info);

280:   info=KSPGetPC( ksp1,&pc1); CHKERRQ(info);
281:   info=KSPGetPC( ksp2,&pc2); CHKERRQ(info);

283:   info=KSPSetFromOptions(ksp2); CHKERRQ(info);

285:   info=PCGetType(pc1,&ptype); CHKERRQ(info);
286:   if (ptype){
287:     info=PetscStrncmp(ptype,PCMG,3,&flg); CHKERRQ(info);
288:     if (flg){
289:       info=PCSetType(pc2,PCBJACOBI); CHKERRQ(info);
290:     } else {
291:       info=PCSetType(pc2, ptype); CHKERRQ(info);
292:     }
293:   }

295:   *kspout=ksp2;
296:   return(0);
297: }

301: int TaoLinearSolverPetsc::Duplicate(TaoLinearSolver**M){
302:   int info;
303:   KSP ksp1=this->GetKSP();
304:   KSP ksp2;

307:   info = KSPDuplicate(ksp1,&ksp2); CHKERRQ(info);
308:   *M = new TaoLinearSolverPetsc(ksp2);
309:   return(0);
310: }

312: #endif