xref: /petsc/src/ksp/pc/impls/deflation/deflationspace.c (revision 253c68fdaaed99ca2470d074a9941b951bbe7d35) 
1 
2 PetscScalar db2[] = {0.7071067811865476,0.7071067811865476};
3 
4 //dec low high
5 PetscScalar db4[] = {-0.12940952255092145,0.22414386804185735,0.836516303737469,0.48296291314469025};
6 //rec low high
7 //PetscScalar db4[] = {0.48296291314469025,0.836516303737469,0.22414386804185735,-0.12940952255092145};
8 
9 PetscScalar db8[] = {-0.010597401784997278,
10 0.032883011666982945,
11 0.030841381835986965,
12 -0.18703481171888114,
13 -0.02798376941698385,
14 0.6308807679295904,
15 0.7148465705525415,
16 0.23037781330885523};
17 
18 //PetscScalar db8[] = {0.23037781330885523,
19 //0.7148465705525415,
20 //0.6308807679295904,
21 //-0.02798376941698385,
22 //-0.18703481171888114,
23 //0.030841381835986965,
24 //0.032883011666982945,
25 //-0.010597401784997278};
26 
27 PetscScalar db16[] = {-0.00011747678400228192,
28 0.0006754494059985568,
29 -0.0003917403729959771,
30 -0.00487035299301066,
31 0.008746094047015655,
32 0.013981027917015516,
33 -0.04408825393106472,
34 -0.01736930100202211,
35 0.128747426620186,
36 0.00047248457399797254,
37 -0.2840155429624281,
38 -0.015829105256023893,
39 0.5853546836548691,
40 0.6756307362980128,
41 0.3128715909144659,
42 0.05441584224308161};
43 
44 //PetscScalar db16[] = {0.05441584224308161,
45 //0.3128715909144659,
46 //0.6756307362980128,
47 //0.5853546836548691,
48 //-0.015829105256023893,
49 //-0.2840155429624281,
50 //0.00047248457399797254,
51 //0.128747426620186,
52 //-0.01736930100202211,
53 //-0.04408825393106472,
54 //0.013981027917015516,
55 //0.008746094047015655,
56 //-0.00487035299301066,
57 //-0.0003917403729959771,
58 //0.0006754494059985568,
59 //-0.00011747678400228192};
60 
61 PetscScalar biorth22[] = {0.0,
62 -0.1767766952966369,
63 0.3535533905932738,
64 1.0606601717798214,
65 0.3535533905932738,
66 -0.1767766952966369};
67 
68 //PetscScalar biorth22[] = {0.0,
69 //0.3535533905932738,
70 //0.7071067811865476,
71 //0.3535533905932738,
72 //0.0,
73 //0.0};
74 
75 PetscScalar meyer[] = {0.0,-1.009999956941423e-12,8.519459636796214e-09,-1.111944952595278e-08,-1.0798819539621958e-08,6.066975741351135e-08,-1.0866516536735883e-07,8.200680650386481e-08,1.1783004497663934e-07,-5.506340565252278e-07,1.1307947017916706e-06,-1.489549216497156e-06,7.367572885903746e-07,3.20544191334478e-06,-1.6312699734552807e-05,6.554305930575149e-05,-0.0006011502343516092,-0.002704672124643725,0.002202534100911002,0.006045814097323304,-0.006387718318497156,-0.011061496392513451,0.015270015130934803,0.017423434103729693,-0.03213079399021176,-0.024348745906078023,0.0637390243228016,0.030655091960824263,-0.13284520043622938,-0.035087555656258346,0.44459300275757724,0.7445855923188063,0.44459300275757724,-0.035087555656258346,-0.13284520043622938,0.030655091960824263,0.0637390243228016,-0.024348745906078023,-0.03213079399021176,0.017423434103729693,0.015270015130934803,-0.011061496392513451,-0.006387718318497156,0.006045814097323304,0.002202534100911002,-0.002704672124643725,-0.0006011502343516092,6.554305930575149e-05,-1.6312699734552807e-05,3.20544191334478e-06,7.367572885903746e-07,-1.489549216497156e-06,1.1307947017916706e-06,-5.506340565252278e-07,1.1783004497663934e-07,8.200680650386481e-08,-1.0866516536735883e-07,6.066975741351135e-08,-1.0798819539621958e-08,-1.111944952595278e-08,8.519459636796214e-09,-1.009999956941423e-12};
76 
77 static PetscErrorCode PCDeflationCreateSpaceWave(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc,Mat *H)
78 {
79   PetscErrorCode ierr;
80   Mat defl;
81   PetscInt i,j,k,ilo,ihi,*Iidx;
82 
83   PetscFunctionBegin;
84   ierr = PetscMalloc1(ncoeffs,&Iidx);CHKERRQ(ierr);
85 
86   /* TODO pass A instead of PC? */
87   ierr = MatCreate(comm,&defl);CHKERRQ(ierr);
88   ierr = MatSetSizes(defl,m,n,M,N);CHKERRQ(ierr);
89   ierr = MatSetUp(defl);CHKERRQ(ierr);
90   ierr = MatSeqAIJSetPreallocation(defl,ncoeffs,NULL);CHKERRQ(ierr);
91   ierr = MatMPIAIJSetPreallocation(defl,ncoeffs,NULL,ncoeffs,NULL);CHKERRQ(ierr);
92   ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
93   ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
94 
95   /* Alg 735 Taswell: fvecmat */
96   k = ncoeffs -2;
97   if (trunc) k = k/2;
98 
99   ierr = MatGetOwnershipRange(defl,&ilo,&ihi);CHKERRQ(ierr);
100   for (i=0; i<ncoeffs; i++) {
101     Iidx[i] = i+ilo*2 -k;
102     if (Iidx[i] >= N) Iidx[i] = PETSC_MIN_INT;
103   }
104   for (i=ilo; i<ihi; i++) {
105     ierr = MatSetValues(defl,1,&i,ncoeffs,Iidx,coeffs,INSERT_VALUES);CHKERRQ(ierr);
106     for (j=0; j<ncoeffs; j++) {
107       Iidx[j] += 2;
108       if (Iidx[j] >= N) Iidx[j] = PETSC_MIN_INT;
109     }
110   }
111 
112   ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
113   ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
114 
115   ierr = PetscFree(Iidx);CHKERRQ(ierr);
116   *H = defl;
117   PetscFunctionReturn(0);
118 }
119 
120 PetscErrorCode PCDeflationGetSpaceHaar(PC pc,Mat *W,PetscInt size)
121 {
122   PetscErrorCode ierr;
123   Mat A,defl;
124   PetscInt i,j,len,ilo,ihi,*Iidx,m,M;
125   PetscScalar *col,val;
126 
127   PetscFunctionBegin;
128   /* Haar basis wavelet, level=size */
129   len = pow(2,size);
130   ierr = PetscMalloc1(len,&col);CHKERRQ(ierr);
131   ierr = PetscMalloc1(len,&Iidx);CHKERRQ(ierr);
132   val = 1./pow(2,size/2.);
133   for (i=0; i<len; i++) col[i] = val;
134 
135   /* TODO pass A instead of PC? */
136   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */
137   ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
138   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
139   ierr = MatCreate(PetscObjectComm((PetscObject)A),&defl);CHKERRQ(ierr);
140   ierr = MatSetSizes(defl,m,PETSC_DECIDE,M,ceil(M/(float)len));CHKERRQ(ierr);
141   ierr = MatSetUp(defl);CHKERRQ(ierr);
142   ierr = MatSeqAIJSetPreallocation(defl,size,NULL);CHKERRQ(ierr);
143   ierr = MatMPIAIJSetPreallocation(defl,size,NULL,size,NULL);CHKERRQ(ierr);
144   ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
145 
146   ierr = MatGetOwnershipRangeColumn(defl,&ilo,&ihi);CHKERRQ(ierr);
147   for (i=0; i<len; i++) Iidx[i] = i+ilo*len;
148   if (M%len && ihi == (int)ceil(M/(float)len)) ihi -= 1;
149   for (i=ilo; i<ihi; i++) {
150     ierr = MatSetValues(defl,len,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr);
151     for (j=0; j<len; j++) Iidx[j] += len;
152   }
153   if (M%len && ihi+1 == ceil(M/(float)len)) {
154     len = M%len;
155     val = 1./pow(pow(2,len),0.5);
156     for (i=0; i<len; i++) col[i] = val;
157     ierr = MatSetValues(defl,len,Iidx,1,&ihi,col,INSERT_VALUES);CHKERRQ(ierr);
158   }
159 
160   ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
161   ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
162 
163   ierr = PetscFree(col);CHKERRQ(ierr);
164   ierr = PetscFree(Iidx);CHKERRQ(ierr);
165   *W = defl;
166   PetscFunctionReturn(0);
167 }
168 
169 PetscErrorCode PCDeflationGetSpaceWave(PC pc,Mat *W,PetscInt size,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc)
170 {
171   PetscErrorCode ierr;
172   Mat A,*H,defl;
173   PetscInt i,m,M,Mdefl,Ndefl;
174   MPI_Comm comm;
175 
176   PetscFunctionBegin;
177   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
178   ierr = PetscMalloc1(size,&H);CHKERRQ(ierr);
179   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */
180   ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
181   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
182   Mdefl = M;
183   Ndefl = M;
184   for (i=0; i<size; i++) {
185     if (Mdefl%2)  {
186       if (trunc) {
187         Mdefl = (PetscInt)PetscCeilReal(Mdefl/2.);
188       } else {
189         Mdefl = (PetscInt)PetscFloorReal((ncoeffs+Mdefl-1)/2.);
190       }
191     } else {
192       Mdefl = Mdefl/2;
193     }
194     ierr = PCDeflationCreateSpaceWave(comm,PETSC_DECIDE,m,Mdefl,Ndefl,ncoeffs,coeffs,trunc,&H[i]);CHKERRQ(ierr);
195     ierr = MatGetLocalSize(H[i],&m,NULL);CHKERRQ(ierr);
196     Ndefl = Mdefl;
197   }
198   ierr = MatCreateComposite(comm,size,H,&defl);CHKERRQ(ierr);
199   ierr = MatCompositeSetType(defl,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
200   *W = defl;
201 
202   for (i=0; i<size; i++) {
203     ierr = MatDestroy(&H[i]);CHKERRQ(ierr);
204   }
205   ierr = PetscFree(H);CHKERRQ(ierr);
206   PetscFunctionReturn(0);
207 }
208 
209 PetscErrorCode PCDeflationGetSpaceAggregation(PC pc,Mat *W)
210 {
211   PetscErrorCode ierr;
212   Mat A,defl;
213   PetscInt i,ilo,ihi,*Iidx,m,M;
214   PetscScalar *col;
215   MPI_Comm comm;
216 
217   PetscFunctionBegin;
218   /* TODO pass A instead of PC? */
219   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */
220   ierr = MatGetOwnershipRangeColumn(A,&ilo,&ihi);CHKERRQ(ierr);
221   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
222   ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
223   ierr = MPI_Comm_size(comm,&m);CHKERRQ(ierr);
224   ierr = MatCreate(comm,&defl);CHKERRQ(ierr);
225   ierr = MatSetSizes(defl,ihi-ilo,1,M,m);CHKERRQ(ierr);
226   ierr = MatSetUp(defl);CHKERRQ(ierr);
227   ierr = MatSeqAIJSetPreallocation(defl,1,NULL);CHKERRQ(ierr);
228   ierr = MatMPIAIJSetPreallocation(defl,1,NULL,0,NULL);CHKERRQ(ierr);
229   ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
230   ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
231 
232 
233   ierr = PetscMalloc1(ihi-ilo,&col);CHKERRQ(ierr);
234   ierr = PetscMalloc1(ihi-ilo,&Iidx);CHKERRQ(ierr);
235   for (i=ilo; i<ihi; i++) {
236     Iidx[i-ilo] = i;
237     col[i-ilo] = 1;
238   }
239   ierr = MPI_Comm_rank(comm,&i);CHKERRQ(ierr);
240   ierr = MatSetValues(defl,ihi-ilo,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr);
241 
242   ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
243   ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
244 
245   ierr = PetscFree(col);CHKERRQ(ierr);
246   ierr = PetscFree(Iidx);CHKERRQ(ierr);
247   *W = defl;
248   PetscFunctionReturn(0);
249 }
250 
251 PetscErrorCode PCDeflationComputeSpace(PC pc)
252 {
253   PetscErrorCode ierr;
254   Mat defl;
255   PetscBool transp=PETSC_TRUE;
256   PC_Deflation *def = (PC_Deflation*)pc->data;
257 
258   /* TODO valid header */
259   PetscFunctionBegin;
260   if (def->spacesize < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PC_DEFLATION Space size specified: %d",def->spacesize);
261   switch (def->spacetype) {
262     case PC_DEFLATION_SPACE_HAAR:
263       transp = PETSC_FALSE;
264       ierr = PCDeflationGetSpaceHaar(pc,&defl,def->spacesize);CHKERRQ(ierr);break;
265     case PC_DEFLATION_SPACE_DB2:
266       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,2,db2,!def->extendsp);CHKERRQ(ierr);break;
267     case PC_DEFLATION_SPACE_DB4:
268       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,4,db4,!def->extendsp);CHKERRQ(ierr);break;
269     case PC_DEFLATION_SPACE_DB8:
270       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,8,db8,!def->extendsp);CHKERRQ(ierr);break;
271     case PC_DEFLATION_SPACE_DB16:
272       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,16,db16,!def->extendsp);CHKERRQ(ierr);break;
273     case PC_DEFLATION_SPACE_BIORTH22:
274       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,6,biorth22,!def->extendsp);CHKERRQ(ierr);break;
275     case PC_DEFLATION_SPACE_MEYER:
276       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,62,meyer,!def->extendsp);CHKERRQ(ierr);break;
277     case PC_DEFLATION_SPACE_AGGREGATION:
278       transp = PETSC_FALSE;
279       ierr = PCDeflationGetSpaceAggregation(pc,&defl);CHKERRQ(ierr);break;
280     default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PC_DEFLATION Space Type specified");
281   }
282 
283   ierr = PCDeflationSetSpace(pc,defl,transp);CHKERRQ(ierr);
284   ierr = MatDestroy(&defl);CHKERRQ(ierr);
285   PetscFunctionReturn(0);
286 }
287 
288