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