Actual source code: taomat_petsc.c
1: #include "tao_general.h"
3: #ifdef TAO_USE_PETSC
5: #include "taomat_petsc.h"
6: #include "../vector/taovec_petsc.h"
7: #include "../indexset/taois_petsc.h"
9: int MatD_Fischer(Mat, Vec, Vec, Vec, Vec, Vec, Vec, Vec, Vec);
10: int MatD_SFischer(Mat, Vec, Vec, Vec, Vec, double, Vec, Vec, Vec, Vec, Vec);
11: int MatSMFResetRowColumn(Mat,IS,IS);
15: /*@C
16: TaoWrapPetscMat - Creates a new TaoMat object using a PETSc matrix.
18: Input Parameter:
19: + M - a PETSc matrix
20: - MM - the address of a pointer to a TaoMatPetsc
22: Output Parameter:
23: . MM - the address of a pointer to new TaoMat
25: Note:
26: A TaoMatPetsc is an object with the methods of an abstract
27: TaoMat object. A TaoMatPetsc contains an implementation of the TaoMat
28: methods. Routines using these vectors should declare a pointer to
29: a TaoMat, assign this pointer to the address of a TaoMat object,
30: use the pointer to invoke methods on the object, and use this pointer
31: as an argument when calling other routines. This usage is different
32: from the usage of a PETSc Mat. In PETSc, applications will typically
33: declare a Mat, and pass it as an argument into routines. That is,
34: applications will typically declare a pointer to a TaoMat and use the
35: pointer, or declare a Mat and use it directly.
37: Note:
38: The user is repsonsible for destroying the Mat M, in addition to
39: to TaoMatPetsc vector MM. The Mat can be destroyed immediately
40: after this routine.
42: Level: developer
44: .seealso TaoMatGetPetscMat(), TaoMatDestroy()
45: @*/
46: int TaoWrapPetscMat( Mat M, TaoMatPetsc* *MM){
48: if (MM){ *MM = new TaoMatPetsc(M);}
49: return(0);
50: }
54: /*@C
55: TaoMatGetPetscMat - Sets the address of a Mat equal to the location
56: of the underlying Mat structure in this TaoMatPetsc object.
58: + MM - the TaoMatPetsc
59: - M - the address of Mat
61: Output Parameter:
62: . M - the address of a specific Mat
64: Note:
65: This routine does not create a Mat. It sets a pointer
66: to the location of an existing Mat.
68: Note:
69: The underlying PETSc Mat may also be obtained by
70: Mat M = (Mat) TMM->MatObject;
72: Note:
73: The function TaoMatPetsc::GetMatrix() will also return the Mat M, but this
74: routine is safer because it will check the object validity.
76: Level: advanced
78: .seealso TaoVecGetPetscVec()
79: @*/
80: int TaoMatGetPetscMat( TaoMat* MM, Mat *M){
82: if (MM && M){
83: *M=(Mat)MM->MatObject;
84: } else if (M){
85: *M=0;
86: }
87: return(0);
88: }
92: /* @C
93: TaoMatPetsc - Creates a new TaoMat object using a PETSc matrix.
95: Input Parameter:
96: + MM - a PETSc matrix
98: Level: advanced
100: Note:
101: The method TaoMatPetsc::SetMatrix(Mat M, Mat MPre, MatStructure) replaces the matrix MM with M.
102: The method TaoMatPetsc::GetMatrix() returns a pointer to the matrix MM, its preconditioner, and KSP flag
104: @ */
105: TaoMatPetsc::TaoMatPetsc(Mat MM)
106: :TaoMat(){
107: int info;
108: pm=0; pm_pre=0; MatObject=0;
109: preflag=DIFFERENT_NONZERO_PATTERN;
110: this->pmatviewer=PETSC_VIEWER_STDOUT_WORLD;
111: if (MM){
112: this->SetMatrix(MM,MM,this->preflag);
113: info = PetscLogInfo((MM,"Wrap a PETSc Mat to create a TaoMat .\n"));
114: }
115: return;
116: }
120: TaoMatPetsc::~TaoMatPetsc(){
121: int info;
122: if (pm){
123: info = PetscLogInfo((pm,"TAO: Destroy a PETSc Mat within a TaoMat .\n"));
124: PetscObjectDereference((PetscObject)pm);
125: PetscObjectDereference((PetscObject)pm_pre);
126: }
127: pm=0; pm_pre=0; MatObject=0;
128: return;
129: }
133: /* @C
134: GetMatrix - Gets the underlying PETSc matrix information
136: Output Parameters:
137: + M - a PETSc matrix
138: . MPre - a PETSc preconditioning matrix (use PETSC_NULL for no change)
139: - flag - flag associated with KSPSetOperators (PETSC_NULL for no change).
141: Level: advanced
143: .seealso TaoMatPetsc::SetMatrix()
144:
145: @ */
146: int TaoMatPetsc::GetMatrix(Mat *M, Mat *MPre, MatStructure *flag){
148: if (M) *M=this->pm;
149: if (flag) *flag=this->preflag;
150: if (MPre) *MPre=this->pm_pre;
151: return(0);
152: }
157: /* @C
158: SetMatrix - Changes the underlying PETSc matrix structrure
160: Input Parameters:
161: + M - a PETSc matrix
162: . MPre - a PETSc preconditioning matrix (use PETSC_NULL for no change)
163: - flag - flag associated with KSPSetOperators (PETSC_NULL for no change).
165: Level: advanced
167: Note:
168: The method TaoMatPetsc::GetMatrix() returns a pointer to the matrix M
170: @ */
171: int TaoMatPetsc::SetMatrix(Mat M, Mat MPre, MatStructure flag){
172: int info;
173: PetscMap cmap1,rmap1,cmap2,rmap2;
174: PetscTruth flg1,flg2;
177: if (M || MPre){
180: if (M!=MPre){
181: info=MatGetPetscMaps(this->pm,&rmap1,&cmap1);CHKERRQ(info);
182: info=MatGetPetscMaps(this->pm_pre,&rmap2,&cmap2);CHKERRQ(info);
183: info=PetscMapSame(rmap1,rmap1,&flg1);CHKERRQ(info);
184: info=PetscMapSame(rmap2,rmap2,&flg2);CHKERRQ(info);
185: if (flg1==PETSC_FALSE || flg2 == PETSC_FALSE){
186: SETERRQ(1,"TAO Error: Preconditioning Matrix does not match the Solution matrix");
187: }
188: }
189:
190: info = PetscObjectReference((PetscObject)M);CHKERRQ(info);
191: info = PetscObjectReference((PetscObject)MPre);CHKERRQ(info);
192: }
194: if (pm&&pm_pre){
195: info = PetscObjectDereference((PetscObject)pm);CHKERRQ(info);
196: info = PetscObjectDereference((PetscObject)pm_pre);CHKERRQ(info);
197: }
199: this->pm=M;
200: this->pm_pre=MPre;
201: this->preflag=flag;
202: this->MatObject=(void*)M;
204: return(0);
205: }
211: int TaoMatPetsc::Compatible(TaoVec* xx, TaoVec* yy, TaoTruth *flag){
212: int info;
213: PetscMap cmap,rmap,xmap,ymap;
214: PetscTruth flg;
215: Vec x,y;
218: *flag=TAO_TRUE;
219: if (xx==0 || yy==0){
220: *flag=TAO_FALSE;
221: return(0);
222: }
223: x=((TaoVecPetsc*)xx)->GetVec();
224: y=((TaoVecPetsc*)yy)->GetVec();
228: info=MatGetPetscMaps(this->pm,&rmap,&cmap);CHKERRQ(info);
229: info=VecGetPetscMap(x,&xmap);CHKERRQ(info);
230: info=VecGetPetscMap(y,&ymap);CHKERRQ(info);
231: info=PetscMapSame(rmap,xmap,&flg);CHKERRQ(info);
232: if (flg==PETSC_TRUE){
233: info=PetscMapSame(cmap,ymap,&flg);CHKERRQ(info);
234: }
235: if (flg==PETSC_FALSE) *flag=TAO_FALSE;
236: return(0);
237: }
241: int TaoMatPetsc::Clone(TaoMat* *ntm){
242: int info;
243: TaoMatPetsc *nptm;
244: Mat M,MPre,npm,nmpre;
245: MatStructure flag;
248: info = this->GetMatrix(&M,&MPre,&flag); CHKERRQ(info);
249: info = MatDuplicate(M,MAT_COPY_VALUES,&npm); CHKERRQ(info);
250: info = MatAssemblyBegin(npm,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
251: info = MatAssemblyEnd(npm,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
252: if (M!=MPre){
253: info = MatDuplicate(MPre,MAT_COPY_VALUES,&nmpre); CHKERRQ(info);
254: info = MatAssemblyBegin(nmpre,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
255: info = MatAssemblyEnd(nmpre,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
256: }
257: info = TaoWrapPetscMat(npm, &nptm); CHKERRQ(info);
258: *ntm = nptm;
259: info = nptm->SetMatrix(npm,nmpre,flag); CHKERRQ(info);
260: info = MatDestroy(npm); CHKERRQ(info);
261: if (M!=MPre){
262: info = MatDestroy(nmpre); CHKERRQ(info);
263: }
264: return(0);
265: }
269: int TaoMatPetsc::CopyFrom(TaoMat* tm){
270: int info;
271: TaoMatPetsc* tmm=(TaoMatPetsc*)tm;
272: Mat M,MPre;
273: MatStructure flag;
276: info = tmm->GetMatrix(&M,&MPre,&flag);CHKERRQ(info);
277: info = MatCopy(M,pm,SAME_NONZERO_PATTERN);
278: CHKERRQ(info);
279: if (pm!=pm_pre){
280: info = MatCopy(MPre,pm_pre,SAME_NONZERO_PATTERN); CHKERRQ(info);
281: }
282: return(0);
283: }
287: int TaoMatPetsc::GetDimensions( int *m, int *n ){
288: int info;
290: info=MatGetSize(pm,m,n); CHKERRQ(info);
291: return(0);
292: }
296: int TaoMatPetsc::Multiply(TaoVec* tv,TaoVec* tw){
297: Vec vv = (Vec) tv->VecObject;
298: Vec ww = (Vec) tw->VecObject;
299: int info;
301: info=MatMult(this->pm,vv,ww); CHKERRQ(info);
302: return(0);
303: }
307: int TaoMatPetsc::MultiplyAdd(TaoVec* tv,TaoVec* tw,TaoVec* ty){
308: Vec yy = (Vec) ty->VecObject;
309: Vec ww = (Vec) tw->VecObject;
310: Vec xx = (Vec) tv->VecObject;
311: int info;
313: info=MatMultAdd(this->pm,xx,ww,yy); CHKERRQ(info);
314: return(0);
315: }
319: int TaoMatPetsc::MultiplyTranspose(TaoVec* tv,TaoVec* tw){
320: Vec vv = (Vec) tv->VecObject;
321: Vec ww = (Vec) tw->VecObject;
322: int info;
324: info=MatMultTranspose(this->pm,vv,ww); CHKERRQ(info);
325: return(0);
326: }
330: int TaoMatPetsc::MultiplyTransposeAdd(TaoVec* tv,TaoVec* tw,TaoVec* ty){
331: Vec yy = (Vec) ty->VecObject;
332: Vec ww = (Vec) tw->VecObject;
333: Vec xx = (Vec) tv->VecObject;
334: int info;
336: info=MatMultTransposeAdd(this->pm,xx,ww,yy); CHKERRQ(info);
337: return(0);
338: }
342: int TaoMatPetsc::SetDiagonal(TaoVec* tv){
343: Vec vv = (Vec) tv->VecObject;
344: int info;
346: info = MatDiagonalSet(this->pm,vv,INSERT_VALUES); CHKERRQ(info);
347: if (this->pm!=this->pm_pre){
348: info = MatDiagonalSet(this->pm_pre,vv,INSERT_VALUES); CHKERRQ(info);
349: }
350: return(0);
351: }
355: int TaoMatPetsc::AddDiagonal(TaoVec* tv){
356: Vec vv = (Vec) tv->VecObject;
357: int info;
359: info = MatDiagonalSet(this->pm,vv,ADD_VALUES); CHKERRQ(info);
360: if (this->pm!=this->pm_pre){
361: info = MatDiagonalSet(this->pm_pre,vv,ADD_VALUES); CHKERRQ(info);
362: }
363: return(0);
364: }
368: int TaoMatPetsc::GetDiagonal(TaoVec* tv){
369: Vec vv = (Vec) tv->VecObject;
370: int info;
372: info = MatGetDiagonal(this->pm,vv); CHKERRQ(info);
373: return(0);
374: }
378: int TaoMatPetsc::ShiftDiagonal(double c){
379: int info;
380: PetscScalar cc=c;
382: info = MatShift(pm,cc); CHKERRQ(info);
383: if (this->pm!=this->pm_pre){
384: info = MatShift(this->pm_pre,cc); CHKERRQ(info);
385: }
386: return(0);
387: }
391: /*@C
392: SetPetscViewer
394: Input Parameter:
395: . viewer - a viewer object to be used with this->View() and MatView()
397: Level: advanced
399: @*/
400: int TaoMatPetsc::SetPetscViewer(PetscViewer viewer){
402: this->pmatviewer= viewer;
403: return(0);
404: }
409: int TaoMatPetsc::View(){
410: int info;
411: PetscTruth flg=PETSC_FALSE;
413: info = PetscOptionsHasName(PETSC_NULL,"-tao_mat_draw",&flg); CHKERRQ(info);
414: if (flg){
415: info = MatView(this->pm,PETSC_VIEWER_DRAW_WORLD); CHKERRQ(info);
416: } else {
417: info = MatView(this->pm,this->pmatviewer); CHKERRQ(info);
418: }
419: return(0);
420: }
425: int TaoMatPetsc::RowScale(TaoVec* tv){
426: Vec vv = (Vec) tv->VecObject;
427: int info;
430: info = MatDiagonalScale(pm, vv, PETSC_NULL); CHKERRQ(info);
431: if (this->pm!=this->pm_pre){
432: info = MatDiagonalScale(this->pm_pre,vv,PETSC_NULL); CHKERRQ(info);
433: }
434: return(0);
435: }
439: int TaoMatPetsc::ColScale(TaoVec* tv){
440: Vec vv = ((TaoVecPetsc *)tv)->GetVec();
441: int info;
444: info = MatDiagonalScale(pm, PETSC_NULL, vv); CHKERRQ(info);
445: if (this->pm!=this->pm_pre){
446: info = MatDiagonalScale(this->pm_pre,PETSC_NULL,vv); CHKERRQ(info);
447: }
448: return(0);
449: }
453: int TaoMatPetsc::CreateReducedMatrix(TaoIndexSet* S1,TaoIndexSet* S2,TaoMat* *MM){
454: int info;
455: TaoMatPetsc *M;
456: Mat B=0,BPre,A,APre;
457: TaoPetscISType type;
458: TaoIndexSetPetsc *TRows=((TaoIndexSetPetsc *)S1);
459: MatStructure flag;
460: PetscTruth assembled,flg;
463: info = TRows->GetReducedType(&type); CHKERRQ(info);
464: if (type==TaoMaskFullSpace){
465: info = this->GetMatrix(&B,&BPre,&flag); CHKERRQ(info);
466: PetscOptionsHasName(0,"-different_submatrix",&flg);
467: if (flg==PETSC_TRUE){
468: info = TaoWrapPetscMat(B, &M); CHKERRQ(info);
469: info = MatAssembled(B,&assembled); CHKERRQ(info);
470: if (assembled==PETSC_TRUE){
471: info = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&A); CHKERRQ(info);
472: info = M->SetMatrix(A,A,flag); CHKERRQ(info);
473: info = MatDestroy(A); CHKERRQ(info);
474: }
475: info = MatAssembled(BPre,&assembled); CHKERRQ(info);
476: if (B != BPre && assembled==PETSC_TRUE){
477: info = MatDuplicate(BPre,MAT_DO_NOT_COPY_VALUES,&APre); CHKERRQ(info);
478: info = M->SetMatrix(A,APre,flag); CHKERRQ(info);
479: info = MatDestroy(APre); CHKERRQ(info);
480: }
481: } else {
482: info = this->GetMatrix(&B,&BPre,&flag); CHKERRQ(info);
483: info = TaoWrapPetscMat(B, &M); CHKERRQ(info);
484: info = M->SetMatrix(B,BPre,flag); CHKERRQ(info);
485: }
486:
487: } else {
488: info = this->GetMatrix(&B,&BPre,&flag); CHKERRQ(info);
489: info = TaoWrapPetscMat(B, &M); CHKERRQ(info);
490: info = M->SetMatrix(B,BPre,flag); CHKERRQ(info);
491: }
492: *MM=M;
494: return(0);
495: }
499: int TaoMatPetsc::SetReducedMatrix(TaoMat* M,TaoIndexSet* S1,TaoIndexSet* S2){
501: int info;
502: TaoPetscISType type;
503: TaoIndexSetPetsc *TRows=((TaoIndexSetPetsc *)S1);
504: TaoIndexSetPetsc *TCols=((TaoIndexSetPetsc *)S2);
505: TaoMatPetsc *MM=((TaoMatPetsc *)M);
506: Mat A,Apre,B,BPre;
507: MatStructure flag;
510: info=this->GetMatrix(&B,&BPre,0); CHKERRQ(info);
511: info=MM->GetMatrix(&A,&Apre,0); CHKERRQ(info);
512: info = TRows->GetReducedType(&type); CHKERRQ(info);
513: if (type==TaoMaskFullSpace){
514: Vec DDiag,VRow,VCol;
515: info = this->GetMatrix(&B,&BPre,&flag); CHKERRQ(info);
516: info = MM->GetMatrix(&A,&Apre,&flag); CHKERRQ(info);
517: if (!A){
518: info = MatDuplicate(B,MAT_COPY_VALUES,&A); CHKERRQ(info);
519: if (B!=BPre){info = MatDuplicate(BPre,MAT_COPY_VALUES,&Apre); CHKERRQ(info);}
520: info = MM->SetMatrix(A,Apre,flag); CHKERRQ(info);
521: info=MatDestroy(A); CHKERRQ(info);
522: if (B!=BPre){info=MatDestroy(Apre); CHKERRQ(info);}
523: }
524: if (A!=B){ info=MatCopy(A,B,SAME_NONZERO_PATTERN); CHKERRQ(info); }
525: if (B!=BPre && Apre!=BPre){ info=MatCopy(Apre,BPre,SAME_NONZERO_PATTERN); CHKERRQ(info); }
526: info=TRows->GetMask(&VRow); CHKERRQ(info);
527: info=TCols->GetMask(&VCol); CHKERRQ(info);
528: info=VecDuplicate(VRow,&DDiag); CHKERRQ(info);
529: info=MatGetDiagonal(A,DDiag); CHKERRQ(info);
530: info=MatDiagonalScale(B,VRow,VCol); CHKERRQ(info);
531: info=MatDiagonalSet(B,DDiag,INSERT_VALUES); CHKERRQ(info);
532: if (B!=BPre){
533: info=MatGetDiagonal(Apre,DDiag); CHKERRQ(info);
534: info=MatDiagonalScale(BPre,VRow,VCol); CHKERRQ(info);
535: info=MatDiagonalSet(BPre,DDiag,INSERT_VALUES); CHKERRQ(info);
536: }
537: info=VecDestroy(DDiag); CHKERRQ(info);
538: } else {
539: IS LocalRows, LocalCols, AllCols;
540: int nn;
541: // info = this->GetMatrix(&B,&BPre,&flag); CHKERRQ(info);
542: info = MM->GetMatrix(&A,&Apre,&flag); CHKERRQ(info);
543: info = TRows->RedistributeIS(&LocalRows); CHKERRQ(info);
544: info = TCols->RedistributeIS(&LocalCols); CHKERRQ(info);
545: info = TCols->GetWholeIS(&AllCols); CHKERRQ(info);
546: info = ISGetLocalSize(LocalCols,&nn); CHKERRQ(info);
547: info = MatGetSubMatrix(A,LocalRows,AllCols,nn,MAT_INITIAL_MATRIX,&B);
548: CHKERRQ(info);
549: BPre=B;
551: if (A != Apre){
552: info = MatGetSubMatrix(Apre,LocalRows,AllCols,nn,MAT_INITIAL_MATRIX,&BPre);
553: CHKERRQ(info);
554: } else {
555: info = PetscObjectReference((PetscObject)BPre);CHKERRQ(info);
556: }
557:
558: info = this->SetMatrix(B,BPre,flag); CHKERRQ(info);
559: info = PetscObjectDereference((PetscObject)B);CHKERRQ(info);
560: info = PetscObjectDereference((PetscObject)BPre);CHKERRQ(info);
561: }
563: return(0);
564: }
568: int TaoMatPetsc::D_Fischer(TaoVec* tx, TaoVec* tf, TaoVec* tl, TaoVec* tu,
569: TaoVec* tt1, TaoVec* tt2, TaoVec* tda, TaoVec* tdb)
570: {
571: Vec x = ((TaoVecPetsc *)tx)->GetVec();
572: Vec f = ((TaoVecPetsc *)tf)->GetVec();
573: Vec l = ((TaoVecPetsc *)tl)->GetVec();
574: Vec u = ((TaoVecPetsc *)tu)->GetVec();
575: Vec da = ((TaoVecPetsc *)tda)->GetVec();
576: Vec db = ((TaoVecPetsc *)tdb)->GetVec();
577: Vec t1 = ((TaoVecPetsc *)tt1)->GetVec();
578: Vec t2 = ((TaoVecPetsc *)tt2)->GetVec();
579: int info;
582: info = MatD_Fischer(pm, x, f, l, u, t1, t2, da, db); CHKERRQ(info);
583: return(0);
584: }
588: int TaoMatPetsc::D_SFischer(TaoVec* tx, TaoVec* tf, TaoVec* tl, TaoVec* tu,
589: double mu,
590: TaoVec* tt1, TaoVec* tt2,
591: TaoVec* tda, TaoVec* tdb, TaoVec* tdm)
592: {
593: Vec x = ((TaoVecPetsc *)tx)->GetVec();
594: Vec f = ((TaoVecPetsc *)tf)->GetVec();
595: Vec l = ((TaoVecPetsc *)tl)->GetVec();
596: Vec u = ((TaoVecPetsc *)tu)->GetVec();
597: Vec t1 = ((TaoVecPetsc *)tt1)->GetVec();
598: Vec t2 = ((TaoVecPetsc *)tt2)->GetVec();
599: Vec da = ((TaoVecPetsc *)tda)->GetVec();
600: Vec db = ((TaoVecPetsc *)tdb)->GetVec();
601: Vec dm = ((TaoVecPetsc *)tdm)->GetVec();
602: int info;
605: if ((mu >= -TAO_EPSILON) && (mu <= TAO_EPSILON)) {
606: tdm->SetToZero();
607: D_Fischer(tx, tf, tl, tu, tt1, tt2, tda, tdb);
608: }
609: else {
610: info = MatD_SFischer(pm, x, f, l, u, mu, t1, t2, da, db, dm); CHKERRQ(info);
611: }
612: return(0);
613: }
617: int TaoMatPetsc::Norm1(double *nm){
618: int info;
619: PetscReal nnmm;
621: info = MatNorm(pm,NORM_1,&nnmm);CHKERRQ(info);
622: *nm=nnmm;
623: return(0);
624: }
626: #endif