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