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