Actual source code: taolinearsolver_petsc.c

  1: #include "tao_general.h"

  3: #ifdef TAO_USE_PETSC

  5: #include "taolinearsolver_petsc.h"

  7: #include "../matrix/taomat_petsc.h"
  8: #include "../vector/taovec_petsc.h"

 10: TaoLinearSolverPetsc::TaoLinearSolverPetsc(KSP S):TaoLinearSolver(){
 11:   linear_its=0;
 12:   ksp=0;
 13:   this->pkspviewer=PETSC_VIEWER_STDOUT_WORLD;
 14:   this->SetKSP(S);
 15:   return;
 16: }

 18: TaoLinearSolverPetsc::~TaoLinearSolverPetsc(){
 19:   if (ksp){
 20:     KSPDestroy(ksp);
 21:   }
 22:   return;
 23: }


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

 31:    Input Parameter:
 32: .  S -  a KSP

 34:    Output Parameter:
 35: .  SS - new TaoMat

 37:    Note:  
 38:    A TaoLinearSolverPetsc is an object with the methods of an abstract
 39:    TaoLinearSolver object.  A TaoLinearSolverPetsc contains an implementation 
 40:    of the TaoLinearSolver methods.  Routines using these vectors should 
 41:    declare a pointer to a TaoLinearSolver, assign this pointer to the address 
 42:    of a TaoLinearSolver object, use the pointer to invoke methods on the 
 43:    object, and use this pointer as an argument when calling other routines.  
 44:    This usage is different from the usage of a PETSc KSP.  In PETSc, 
 45:    applications will typically declare a KSP, and pass it as an argument into 
 46:    routines.  That is, applications will typically declare a pointer to a 
 47:    TaoLinearSolver and use the pointer, or declare a KSP and use it directly.

 49:    Level: developer

 51: @*/
 52: int TaoWrapKSP( KSP S, TaoLinearSolverPetsc ** SS){
 55:   *SS = new TaoLinearSolverPetsc(S);
 56:   return(0);
 57: }

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

 64:    Input Parameter:
 65: +  TS - the TaoLinearSolver
 66: -  S -  the address of KSP

 68:    Output Parameter:
 69: .  S -  address of the underlying PETSc KSP object


 72:    Note:  
 73:    This routine does not create a KSP.  It sets a pointer
 74:    to the location of an existing KSP.

 76:    Level: advanced

 78: .seealso TaoVecGetPetscVec(), TaoMatGetPetscMat(), TaoAppGetKSP()
 79: @*/
 80: int TaoLinearSolverGetKSP( TaoLinearSolver *TS, KSP * S){
 82:   if (TS){
 83:     *S=((TaoLinearSolverPetsc *)TS)->GetKSP();
 85:   } else {
 86:     *S=0;
 87:   }
 88:   return(0);
 89: }


 94: int TaoLinearSolverPetsc::SetKSP(KSP ksp2){
 95:   int info;
 97:   if (ksp2){
 99:     PetscObjectReference((PetscObject)ksp2); 
100:   }
101:   if (ksp){
102:     info=KSPDestroy(ksp); CHKERRQ(info);
103:   }
104:   ksp=ksp2;
105:   return(0);
106: }

110: int TaoLinearSolverPetsc::GetKSP(KSP *ksp2){
112:   if (ksp2){
113:     *ksp2=this->ksp;
114:   }
115:   return(0);
116: }

120: int TaoLinearSolverPetsc::PreSolve(TaoMat* m1)
121: {
122:   TaoMatPetsc *MM = dynamic_cast <TaoMatPetsc *> (m1);
123:   Mat mm, mm_pre;
124:   MatStructure preflag;
125:   int info;

128:   info = MM->GetMatrix(&mm, &mm_pre, &preflag); CHKERRQ(info);
129:   info = KSPSetOperators(ksp, mm, mm_pre, preflag); CHKERRQ(info);
130:   return(0);
131: }

135: int TaoLinearSolverPetsc::SetTrustRadius(double rad)
136: {
137:   const KSPType ktype;
138:   int info;
139:   PetscTruth flg;


143:   info = KSPGetType(ksp, &ktype); CHKERRQ(info);

145:   info = PetscStrcmp((char *)ktype, KSPNASH, &flg); CHKERRQ(info);
146:   if (flg == PETSC_TRUE) {         
147:     info = KSPNASHSetRadius(ksp, rad); CHKERRQ(info);
148:   }

150:   info = PetscStrcmp((char *)ktype, KSPSTCG, &flg); CHKERRQ(info);
151:   if (flg == PETSC_TRUE) {         
152:     info = KSPSTCGSetRadius(ksp, rad); CHKERRQ(info);
153:   }

155:   info = PetscStrcmp((char *)ktype, KSPGLTR, &flg); CHKERRQ(info);
156:   if (flg == PETSC_TRUE) {         
157:     info = KSPGLTRSetRadius(ksp, rad); CHKERRQ(info);
158:   }

160:   return(0);
161: }

165: int TaoLinearSolverPetsc::GetNormDirection(double *norm_d)
166: {
167:   const KSPType ktype;
168:   int info;
169:   PetscTruth flg;


173:   info = KSPGetType(ksp, &ktype); CHKERRQ(info);

175:   info = PetscStrcmp((char *)ktype, KSPNASH, &flg); CHKERRQ(info);
176:   if (flg == PETSC_TRUE) {         
177:     info = KSPNASHGetNormD(ksp, norm_d); CHKERRQ(info);
178:   }

180:   info = PetscStrcmp((char *)ktype, KSPSTCG, &flg); CHKERRQ(info);
181:   if (flg == PETSC_TRUE) {         
182:     info = KSPSTCGGetNormD(ksp, norm_d); CHKERRQ(info);
183:   }

185:   info = PetscStrcmp((char *)ktype, KSPGLTR, &flg); CHKERRQ(info);
186:   if (flg == PETSC_TRUE) {         
187:     info = KSPGLTRGetNormD(ksp, norm_d); CHKERRQ(info);
188:   }

190:   return(0);
191: }

195: int TaoLinearSolverPetsc::GetObjFcn(double *o_fcn)
196: {
197:   const KSPType ktype;
198:   int info;
199:   PetscTruth flg;


203:   info = KSPGetType(ksp, &ktype); CHKERRQ(info);

205:   info = PetscStrcmp((char *)ktype, KSPNASH, &flg); CHKERRQ(info);
206:   if (flg == PETSC_TRUE) {         
207:     info = KSPNASHGetObjFcn(ksp, o_fcn); CHKERRQ(info);
208:   }

210:   info = PetscStrcmp((char *)ktype, KSPSTCG, &flg); CHKERRQ(info);
211:   if (flg == PETSC_TRUE) {         
212:     info = KSPSTCGGetObjFcn(ksp, o_fcn); CHKERRQ(info);
213:   }

215:   info = PetscStrcmp((char *)ktype, KSPGLTR, &flg); CHKERRQ(info);
216:   if (flg == PETSC_TRUE) {         
217:     info = KSPGLTRGetObjFcn(ksp, o_fcn); CHKERRQ(info);
218:   }

220:   return(0);
221: }

225: int TaoLinearSolverPetsc::GetMinEig(double *e_min)
226: {
227:   const KSPType ktype;
228:   int info;
229:   PetscTruth flg;


233:   *e_min = 0.0;

235:   info = KSPGetType(ksp, &ktype); CHKERRQ(info);
236:   info = PetscStrcmp((char *)ktype, KSPGLTR, &flg); CHKERRQ(info);
237:   if (flg == PETSC_TRUE) {         
238:     info = KSPGLTRGetMinEig(ksp, e_min); CHKERRQ(info);
239:   }

241:   return(0);
242: }

246: int TaoLinearSolverPetsc::GetLambda(double *lambda)
247: {
248:   const KSPType ktype;
249:   int info;
250:   PetscTruth flg;


254:   *lambda = 0.0;

256:   info = KSPGetType(ksp, &ktype); CHKERRQ(info);
257:   info = PetscStrcmp((char *)ktype, KSPGLTR, &flg); CHKERRQ(info);
258:   if (flg == PETSC_TRUE) {         
259:     info = KSPGLTRGetLambda(ksp, lambda); CHKERRQ(info);
260:   }

262:   return(0);
263: }

267: int TaoLinearSolverPetsc::Solve(TaoVec* tv, TaoVec* tw, TaoTruth *flag)
268: {
269:   TaoVecPetsc *pv = dynamic_cast <TaoVecPetsc *> (tv);
270:   TaoVecPetsc *pw = dynamic_cast <TaoVecPetsc *> (tw);
271:   int info;
272:   PetscInt its;


276:   info = KSPSolve(ksp,pv->GetVec(),pw->GetVec()); CHKERRQ(info);
277:   info = KSPGetIterationNumber(ksp, &its); CHKERRQ(info);

279:   this->linear_its=PetscMax(its,-its);
280:   if (its>0) *flag=TAO_TRUE; 
281:   else *flag=TAO_FALSE;

283:   return(0);
284: }

288: int TaoLinearSolverPetsc::SolveTrustRegion(TaoVec *b, TaoVec *x, 
289:                                            double r, TaoTruth *flg) 
290: {
291:   int info;

294:   info = this->SetTrustRadius(r); CHKERRQ(info);
295:   info = this->Solve(b, x, flg); CHKERRQ(info);
296:   return(0);
297: }

301: int TaoLinearSolverPetsc::GetNumberIterations(int * iters){
303:   *iters=linear_its;
304:   return(0);
305: }


310: int TaoLinearSolverPetsc::SetTolerances(double rtol, double atol,
311:                                         double dtol, int maxits) {
312:   int info;
314:   info = KSPSetTolerances(ksp, rtol, atol, dtol, maxits); CHKERRQ(info);
315:   return(0);
316: }

320: int TaoLinearSolverPetsc::View(){
321:   int info;
323:   info=KSPView(this->ksp,this->pkspviewer);CHKERRQ(info);
324:   return(0);
325: }

329: int TaoLinearSolverPetsc::SetOptions(){
330:   int info;
332:   info=KSPSetFromOptions(ksp); CHKERRQ(info);
333:   return(0);
334: }


337: #endif