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