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