Actual source code: multigridmat.c

      1: #include "multigridmat.h"

      3: #include "../../indexset/taois_petsc.h"

      5: TaoMultiGridMatPetsc::TaoMultiGridMatPetsc(Mat MM):TaoMatPetsc(MM){
      6:   int i,nn=PETSCDAAPPMAXGRIDS;
      8:   for (i=0;i<nn;i++){
      9:     this->grid[i].da=0;
     10:     this->grid[i].H=0;
     11:     this->grid[i].R=0;
     12:     this->grid[i].RHS=0;
     13:     this->grid[i].Interpolate=0;
     14:     this->grid[i].CScale=0;
     15:     this->grid[i].W3=0;
     16:     this->grid[i].coloring=0;
     17:   }
     18:   this->nda=0;
     19:   return;
     20: }
     21:
     22: TaoMultiGridMatPetsc::~TaoMultiGridMatPetsc(){
     23:   this->TakeDown();
     24:   return;
     25: }
     26:

     30: int TaoMultiGridMatPetsc::TakeDown(){
     31:   int i,info,nn=this->nda;
     33:   for (i=0;i<nn;i++){
     34:     this->grid[i].da=0;
     35:     this->grid[i].H=0;
     36:     this->grid[i].R=0;
     37:     this->grid[i].RHS=0;
     38:     this->grid[i].Interpolate=0;
     39:     this->grid[i].CScale=0;
     40:     if (this->grid[i].W3){
     41:       info=VecDestroy(this->grid[i].W3);CHKERRQ(info);
     42:     }
     43:     this->grid[i].W3=0;
     44:   }
     45:   this->nda=0;
     46:   return(0);
     47: }

     52: int TaoMultiGridMatPetsc::SetUpMultiGrid(GridCtx*dagrid,int nn){
     53:   int i,info;
     55:   info=this->TakeDown(); CHKERRQ(info);
     56:   this->nda=nn;
     57:   this->ndamax=PETSCDAAPPMAXGRIDS;
     58:   for (i=0;i<nn;i++){
     59:     this->grid[i].H = dagrid[i].H;
     60:     this->grid[i].R = dagrid[i].R;
     61:     this->grid[i].RHS = dagrid[i].RHS;
     62:     this->grid[i].Interpolate = dagrid[i].Interpolate;
     63:     this->grid[i].CScale = dagrid[i].CScale;
     64:     this->grid[i].W3=0;
     65:     info=VecDuplicate(this->grid[i].R,&this->grid[i].W3); CHKERRQ(info);
     66:   }
     67:   info = this->SetMatrix(this->grid[nn-1].H,this->grid[nn-1].H,SAME_NONZERO_PATTERN); CHKERRQ(info);
     68:   info=TaoSelectSubset(TaoMaskFullSpace); CHKERRQ(info);
     69:   return(0);
     70: }

     75: int TaoMultiGridMatPetsc::SetDiagonal(TaoVec*tv){
     76:   int i,info;
     77:   Vec VCourse, VFine=(Vec)tv->VecObject;

     80:   info = MatDiagonalSet(this->grid[this->nda-1].H,VFine,INSERT_VALUES); CHKERRQ(info);
     81:   for (i=this->nda-1;i>0;i--){
     82:     VCourse=this->grid[i-1].R;
     83:     info=MatMultTranspose(this->grid[i].Interpolate,VFine,VCourse); CHKERRQ(info);
     84:     info=VecPointwiseMult(VCourse,this->grid[i].CScale,VCourse); CHKERRQ(info);
     85:     info=MatDiagonalSet(this->grid[i-1].H,VCourse,INSERT_VALUES); CHKERRQ(info);
     86:     VFine=VCourse;
     87:   }
     88:   return(0);
     89: }

     93: int TaoMultiGridMatPetsc::AddDiagonal(TaoVec*tv){
     94:   int i,info;
     95:   Vec VCourse,VFine = (Vec) tv->VecObject;

     98:   info = MatDiagonalSet(this->pm,VFine,ADD_VALUES); CHKERRQ(info);
     99:   for (i=this->nda-1;i>0;i--){
    100:     VCourse=this->grid[i-1].R;
    101:     info=MatMultTranspose(this->grid[i].Interpolate,VFine,VCourse); CHKERRQ(info);
    102:     info=VecPointwiseMult(VCourse,this->grid[i].CScale,VCourse); CHKERRQ(info);
    103:     info=MatDiagonalSet(this->grid[i-1].H,VCourse,ADD_VALUES); CHKERRQ(info);
    104:     VFine=VCourse;
    105:   }
    106:   return(0);
    107: }

    112: int TaoMultiGridMatPetsc::ShiftDiagonal(double c){
    113:   int i,info;
    114:   PetscScalar cc=c;
    116:   info = MatShift(pm,cc); CHKERRQ(info);
    117:   for (i=this->nda-1;i>0;i--){
    118:     info = MatShift(this->grid[i-1].H,cc); CHKERRQ(info);
    119:   }
    120:   return(0);
    121: }

    125: int TaoMultiGridMatPetsc::RowScale(TaoVec*tv){
    126:   int i,info;
    127:   Vec VCourse,VFine = (Vec) tv->VecObject;

    130:   info = MatDiagonalScale(this->pm,VFine,PETSC_NULL); CHKERRQ(info);
    131:   for (i=this->nda-1;i>0;i--){
    132:     VCourse=this->grid[i-1].R;
    133:     info=MatMultTranspose(this->grid[i].Interpolate,VFine,VCourse); CHKERRQ(info);
    134:     info=VecPointwiseMult(VCourse,this->grid[i].CScale,VCourse); CHKERRQ(info);
    135:     info=MatDiagonalScale(this->grid[i-1].H,VCourse,PETSC_NULL); CHKERRQ(info);
    136:     VFine=VCourse;
    137:   }
    138:   return(0);
    139: }

    143: int TaoMultiGridMatPetsc::ColScale(TaoVec*tv){
    144:   int i,info;
    145:   Vec VCourse,VFine = (Vec) tv->VecObject;

    148:   info = MatDiagonalScale(this->pm,PETSC_NULL,VFine); CHKERRQ(info);
    149:   for (i=this->nda-1;i>0;i--){
    150:     VCourse=this->grid[i-1].R;
    151:     info=MatMultTranspose(this->grid[i].Interpolate,VFine,VCourse); CHKERRQ(info);
    152:     info=VecPointwiseMult(VCourse,this->grid[i].CScale,VCourse); CHKERRQ(info);
    153:     info=MatDiagonalScale(this->grid[i-1].H,PETSC_NULL,VCourse); CHKERRQ(info);
    154:     VFine=VCourse;
    155:   }
    156:   return(0);
    157: }

    161: int TaoMultiGridMatPetsc::CreateReducedMatrix(TaoIndexSet*S1,TaoIndexSet*S2,TaoMat**MM){
    162:   int info;
    163:   Mat B,BPre;
    164:   MatStructure flag;
    165:   TaoMultiGridMatPetsc *MMM;

    168:   info = this->GetMatrix(&B,&BPre,&flag); CHKERRQ(info);
    169:   MMM = new TaoMultiGridMatPetsc(B);
    170:   info=MMM->SetUpMultiGrid(this->grid,this->nda); CHKERRQ(info);
    171:   info = MMM->SetMatrix(B,BPre,flag); CHKERRQ(info);
    172:   /*
    173:     info = MatDuplicate(B,BB); CHKERRQ(info);
    174:     info = MMM->SetMatrix(BB,BBPre,flag); CHKERRQ(info);
    175:   */
    176:   *MM=MMM;
    177:   return(0);
    178: }

    182: int TaoMultiGridMatPetsc::SetReducedMatrix(TaoMat*M,TaoIndexSet*S1,TaoIndexSet*S2){

    184:   int i,info;
    185:   TaoIndexSetPetsc *TRows=((TaoIndexSetPetsc *)S1);
    186:   TaoIndexSetPetsc *TCols=((TaoIndexSetPetsc *)S2);
    187:   TaoMatPetsc *MM=(TaoMatPetsc *)M;
    188:   Vec CFine,CCourse,RFine,RCourse;
    189:   Mat A,Apre,B,Bpre;

    192:   info=MM->GetMatrix(&A,&Apre,0); CHKERRQ(info);
    193:   info=this->GetMatrix(&B,&Bpre,0); CHKERRQ(info);
    194:   info=TRows->GetMask(&RFine); CHKERRQ(info);
    195:   info=TCols->GetMask(&CFine); CHKERRQ(info);
    196:   if (A!=B){
    197:     info=MatCopy(A,B,SAME_NONZERO_PATTERN); CHKERRQ(info);
    198:   }
    199:   info=MatGetDiagonal(B,this->grid[this->nda-1].W3); CHKERRQ(info);
    200:   info=MatDiagonalScale(B,RFine,CFine); CHKERRQ(info);
    201:   info=MatDiagonalSet(B,this->grid[this->nda-1].W3,INSERT_VALUES); CHKERRQ(info);
    202:   if (B!=Bpre){
    203:     info=MatCopy(Apre,Bpre,SAME_NONZERO_PATTERN); CHKERRQ(info);
    204:     info=MatGetDiagonal(Bpre,this->grid[this->nda-1].W3); CHKERRQ(info);
    205:     info=MatDiagonalScale(Bpre,RFine,CFine); CHKERRQ(info);
    206:     info=MatDiagonalSet(Bpre,this->grid[this->nda-1].W3,INSERT_VALUES); CHKERRQ(info);
    207:   }
    208:   for (i=this->nda-1;i>0;i--){
    209:     CCourse=this->grid[i-1].R;
    210:     RCourse=this->grid[i-1].RHS;
    211:     info=MatMultTranspose(this->grid[i].Interpolate,CFine,CCourse); CHKERRQ(info);
    212:     info=MatMultTranspose(this->grid[i].Interpolate,RFine,RCourse); CHKERRQ(info);
    213:     info=VecPointwiseMult(CCourse,this->grid[i].CScale,CCourse); CHKERRQ(info);
    214:     info=VecPointwiseMult(RCourse,this->grid[i].CScale,RCourse); CHKERRQ(info);
    215:     info=MatGetDiagonal(this->grid[i-1].H,this->grid[i-1].W3); CHKERRQ(info);
    216:     info=MatDiagonalScale(this->grid[i-1].H,CCourse,RCourse); CHKERRQ(info);
    217:     info=MatDiagonalSet(this->grid[i-1].H,this->grid[i-1].W3,INSERT_VALUES); CHKERRQ(info);
    218:     RFine=RCourse;
    219:     CFine=CCourse;
    220:   }
    221:   return(0);
    222: }