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: }