xref: /petsc/src/ksp/pc/impls/deflation/deflationspace.c (revision e53e0a0d0486bbcc0147c336cd92e1a6a50ba191)
1*e53e0a0dSJakub Kruzik #include <../src/ksp/pc/impls/deflation/deflation.h>
2*e53e0a0dSJakub Kruzik 
3*e53e0a0dSJakub Kruzik #if defined(HAVE_SLEPC)
4*e53e0a0dSJakub Kruzik #include <slepceps.h>
5*e53e0a0dSJakub Kruzik #include <slepcbv.h>
6*e53e0a0dSJakub Kruzik #endif
7*e53e0a0dSJakub Kruzik 
8*e53e0a0dSJakub Kruzik const char *const PCDeflationSpaceTypes[] = {
9*e53e0a0dSJakub Kruzik   "haar",
10*e53e0a0dSJakub Kruzik   "jacket-haar",
11*e53e0a0dSJakub Kruzik   "db2",
12*e53e0a0dSJakub Kruzik   "db4",
13*e53e0a0dSJakub Kruzik   "db8",
14*e53e0a0dSJakub Kruzik   "db16",
15*e53e0a0dSJakub Kruzik   "biorth22",
16*e53e0a0dSJakub Kruzik   "meyer",
17*e53e0a0dSJakub Kruzik   "aggregation",
18*e53e0a0dSJakub Kruzik   "slepc",
19*e53e0a0dSJakub Kruzik   "slepc-cheap",
20*e53e0a0dSJakub Kruzik   "user",
21*e53e0a0dSJakub Kruzik   "DdefSpaceType",
22*e53e0a0dSJakub Kruzik   "Ddef_SPACE_",
23*e53e0a0dSJakub Kruzik   0
24*e53e0a0dSJakub Kruzik };
25*e53e0a0dSJakub Kruzik 
26*e53e0a0dSJakub Kruzik PetscScalar db2[] = {0.7071067811865476,0.7071067811865476};
27*e53e0a0dSJakub Kruzik 
28*e53e0a0dSJakub Kruzik //dec low high
29*e53e0a0dSJakub Kruzik PetscScalar db4[] = {-0.12940952255092145,0.22414386804185735,0.836516303737469,0.48296291314469025};
30*e53e0a0dSJakub Kruzik //rec low high
31*e53e0a0dSJakub Kruzik //PetscScalar db4[] = {0.48296291314469025,0.836516303737469,0.22414386804185735,-0.12940952255092145};
32*e53e0a0dSJakub Kruzik 
33*e53e0a0dSJakub Kruzik PetscScalar db8[] = {-0.010597401784997278,
34*e53e0a0dSJakub Kruzik 0.032883011666982945,
35*e53e0a0dSJakub Kruzik 0.030841381835986965,
36*e53e0a0dSJakub Kruzik -0.18703481171888114,
37*e53e0a0dSJakub Kruzik -0.02798376941698385,
38*e53e0a0dSJakub Kruzik 0.6308807679295904,
39*e53e0a0dSJakub Kruzik 0.7148465705525415,
40*e53e0a0dSJakub Kruzik 0.23037781330885523};
41*e53e0a0dSJakub Kruzik 
42*e53e0a0dSJakub Kruzik //PetscScalar db8[] = {0.23037781330885523,
43*e53e0a0dSJakub Kruzik //0.7148465705525415,
44*e53e0a0dSJakub Kruzik //0.6308807679295904,
45*e53e0a0dSJakub Kruzik //-0.02798376941698385,
46*e53e0a0dSJakub Kruzik //-0.18703481171888114,
47*e53e0a0dSJakub Kruzik //0.030841381835986965,
48*e53e0a0dSJakub Kruzik //0.032883011666982945,
49*e53e0a0dSJakub Kruzik //-0.010597401784997278};
50*e53e0a0dSJakub Kruzik 
51*e53e0a0dSJakub Kruzik PetscScalar db16[] = {-0.00011747678400228192,
52*e53e0a0dSJakub Kruzik 0.0006754494059985568,
53*e53e0a0dSJakub Kruzik -0.0003917403729959771,
54*e53e0a0dSJakub Kruzik -0.00487035299301066,
55*e53e0a0dSJakub Kruzik 0.008746094047015655,
56*e53e0a0dSJakub Kruzik 0.013981027917015516,
57*e53e0a0dSJakub Kruzik -0.04408825393106472,
58*e53e0a0dSJakub Kruzik -0.01736930100202211,
59*e53e0a0dSJakub Kruzik 0.128747426620186,
60*e53e0a0dSJakub Kruzik 0.00047248457399797254,
61*e53e0a0dSJakub Kruzik -0.2840155429624281,
62*e53e0a0dSJakub Kruzik -0.015829105256023893,
63*e53e0a0dSJakub Kruzik 0.5853546836548691,
64*e53e0a0dSJakub Kruzik 0.6756307362980128,
65*e53e0a0dSJakub Kruzik 0.3128715909144659,
66*e53e0a0dSJakub Kruzik 0.05441584224308161};
67*e53e0a0dSJakub Kruzik 
68*e53e0a0dSJakub Kruzik //PetscScalar db16[] = {0.05441584224308161,
69*e53e0a0dSJakub Kruzik //0.3128715909144659,
70*e53e0a0dSJakub Kruzik //0.6756307362980128,
71*e53e0a0dSJakub Kruzik //0.5853546836548691,
72*e53e0a0dSJakub Kruzik //-0.015829105256023893,
73*e53e0a0dSJakub Kruzik //-0.2840155429624281,
74*e53e0a0dSJakub Kruzik //0.00047248457399797254,
75*e53e0a0dSJakub Kruzik //0.128747426620186,
76*e53e0a0dSJakub Kruzik //-0.01736930100202211,
77*e53e0a0dSJakub Kruzik //-0.04408825393106472,
78*e53e0a0dSJakub Kruzik //0.013981027917015516,
79*e53e0a0dSJakub Kruzik //0.008746094047015655,
80*e53e0a0dSJakub Kruzik //-0.00487035299301066,
81*e53e0a0dSJakub Kruzik //-0.0003917403729959771,
82*e53e0a0dSJakub Kruzik //0.0006754494059985568,
83*e53e0a0dSJakub Kruzik //-0.00011747678400228192};
84*e53e0a0dSJakub Kruzik 
85*e53e0a0dSJakub Kruzik PetscScalar biorth22[] = {0.0,
86*e53e0a0dSJakub Kruzik -0.1767766952966369,
87*e53e0a0dSJakub Kruzik 0.3535533905932738,
88*e53e0a0dSJakub Kruzik 1.0606601717798214,
89*e53e0a0dSJakub Kruzik 0.3535533905932738,
90*e53e0a0dSJakub Kruzik -0.1767766952966369};
91*e53e0a0dSJakub Kruzik 
92*e53e0a0dSJakub Kruzik //PetscScalar biorth22[] = {0.0,
93*e53e0a0dSJakub Kruzik //0.3535533905932738,
94*e53e0a0dSJakub Kruzik //0.7071067811865476,
95*e53e0a0dSJakub Kruzik //0.3535533905932738,
96*e53e0a0dSJakub Kruzik //0.0,
97*e53e0a0dSJakub Kruzik //0.0};
98*e53e0a0dSJakub Kruzik 
99*e53e0a0dSJakub Kruzik PetscReal 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};
100*e53e0a0dSJakub Kruzik 
101*e53e0a0dSJakub Kruzik static PetscErrorCode PCDeflationCreateSpaceJacketHaar(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscBool jacket,Mat *H)
102*e53e0a0dSJakub Kruzik {
103*e53e0a0dSJakub Kruzik   PetscErrorCode ierr;
104*e53e0a0dSJakub Kruzik   Mat defl;
105*e53e0a0dSJakub Kruzik   PetscInt i,j,ilo,ihi,alloc=2,*Iidx;
106*e53e0a0dSJakub Kruzik   PetscReal val,*row;
107*e53e0a0dSJakub Kruzik 
108*e53e0a0dSJakub Kruzik   PetscFunctionBegin;
109*e53e0a0dSJakub Kruzik   if (jacket) alloc = 3;
110*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(alloc,&row);CHKERRQ(ierr);
111*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(alloc,&Iidx);CHKERRQ(ierr);
112*e53e0a0dSJakub Kruzik 
113*e53e0a0dSJakub Kruzik   val = 1./pow(2,0.5);
114*e53e0a0dSJakub Kruzik   row[0] = val;
115*e53e0a0dSJakub Kruzik   row[1] = val;
116*e53e0a0dSJakub Kruzik 
117*e53e0a0dSJakub Kruzik   /* TODO pass A instead of PC? */
118*e53e0a0dSJakub Kruzik   ierr = MatCreate(comm,&defl);CHKERRQ(ierr);
119*e53e0a0dSJakub Kruzik   ierr = MatSetSizes(defl,m,n,M,N);CHKERRQ(ierr);
120*e53e0a0dSJakub Kruzik   ierr = MatSetUp(defl);CHKERRQ(ierr);
121*e53e0a0dSJakub Kruzik   ierr = MatSeqAIJSetPreallocation(defl,alloc,NULL);CHKERRQ(ierr);
122*e53e0a0dSJakub Kruzik   ierr = MatMPIAIJSetPreallocation(defl,alloc,NULL,alloc,NULL);CHKERRQ(ierr);
123*e53e0a0dSJakub Kruzik   ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
124*e53e0a0dSJakub Kruzik   ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
125*e53e0a0dSJakub Kruzik 
126*e53e0a0dSJakub Kruzik   ierr = MatGetOwnershipRange(defl,&ilo,&ihi);CHKERRQ(ierr);
127*e53e0a0dSJakub Kruzik   for (i=0; i<2; i++) Iidx[i] = i+ilo*2;
128*e53e0a0dSJakub Kruzik   if (jacket && ihi==M) ihi -=2;
129*e53e0a0dSJakub Kruzik   if (ihi<ilo) SETERRQ1(comm,PETSC_ERR_ARG_WRONG,"To many cores to assemble Jacket Haar matrix with %d rows",M);
130*e53e0a0dSJakub Kruzik   for (i=ilo; i<ihi; i++) {
131*e53e0a0dSJakub Kruzik     ierr = MatSetValues(defl,1,&i,2,Iidx,row,INSERT_VALUES);CHKERRQ(ierr);
132*e53e0a0dSJakub Kruzik     for (j=0; j<2; j++) Iidx[j] += 2;
133*e53e0a0dSJakub Kruzik   }
134*e53e0a0dSJakub Kruzik   if (jacket && ihi == M-2) {
135*e53e0a0dSJakub Kruzik     for (i=0; i<3; i++) Iidx[i] = i+ilo*2;
136*e53e0a0dSJakub Kruzik     row[0] = 0.5; row[1] = 0.5; row[2] = val;
137*e53e0a0dSJakub Kruzik     ierr = MatSetValues(defl,1,&ihi,3,Iidx,row,INSERT_VALUES);CHKERRQ(ierr);
138*e53e0a0dSJakub Kruzik     ihi += 1;
139*e53e0a0dSJakub Kruzik     row[2] = -row[2];
140*e53e0a0dSJakub Kruzik     ierr = MatSetValues(defl,1,&ihi,3,Iidx,row,INSERT_VALUES);CHKERRQ(ierr);
141*e53e0a0dSJakub Kruzik   }
142*e53e0a0dSJakub Kruzik 
143*e53e0a0dSJakub Kruzik   ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
144*e53e0a0dSJakub Kruzik   ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
145*e53e0a0dSJakub Kruzik 
146*e53e0a0dSJakub Kruzik   ierr = PetscFree(row);CHKERRQ(ierr);
147*e53e0a0dSJakub Kruzik   ierr = PetscFree(Iidx);CHKERRQ(ierr);
148*e53e0a0dSJakub Kruzik   *H = defl;
149*e53e0a0dSJakub Kruzik   PetscFunctionReturn(0);
150*e53e0a0dSJakub Kruzik }
151*e53e0a0dSJakub Kruzik 
152*e53e0a0dSJakub Kruzik static PetscErrorCode PCDeflationCreateSpaceWave(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc,Mat *H)
153*e53e0a0dSJakub Kruzik {
154*e53e0a0dSJakub Kruzik   PetscErrorCode ierr;
155*e53e0a0dSJakub Kruzik   Mat defl;
156*e53e0a0dSJakub Kruzik   PetscInt i,j,k,ilo,ihi,*Iidx;
157*e53e0a0dSJakub Kruzik 
158*e53e0a0dSJakub Kruzik   PetscFunctionBegin;
159*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(ncoeffs,&Iidx);CHKERRQ(ierr);
160*e53e0a0dSJakub Kruzik 
161*e53e0a0dSJakub Kruzik   /* TODO pass A instead of PC? */
162*e53e0a0dSJakub Kruzik   ierr = MatCreate(comm,&defl);CHKERRQ(ierr);
163*e53e0a0dSJakub Kruzik   ierr = MatSetSizes(defl,m,n,M,N);CHKERRQ(ierr);
164*e53e0a0dSJakub Kruzik   ierr = MatSetUp(defl);CHKERRQ(ierr);
165*e53e0a0dSJakub Kruzik   ierr = MatSeqAIJSetPreallocation(defl,ncoeffs,NULL);CHKERRQ(ierr);
166*e53e0a0dSJakub Kruzik   ierr = MatMPIAIJSetPreallocation(defl,ncoeffs,NULL,ncoeffs,NULL);CHKERRQ(ierr);
167*e53e0a0dSJakub Kruzik   ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
168*e53e0a0dSJakub Kruzik   ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
169*e53e0a0dSJakub Kruzik 
170*e53e0a0dSJakub Kruzik   /* Alg 735 Taswell: fvecmat */
171*e53e0a0dSJakub Kruzik   k = ncoeffs -2;
172*e53e0a0dSJakub Kruzik   if (trunc) k = k/2;
173*e53e0a0dSJakub Kruzik 
174*e53e0a0dSJakub Kruzik   ierr = MatGetOwnershipRange(defl,&ilo,&ihi);CHKERRQ(ierr);
175*e53e0a0dSJakub Kruzik   for (i=0; i<ncoeffs; i++) {
176*e53e0a0dSJakub Kruzik     Iidx[i] = i+ilo*2 -k;
177*e53e0a0dSJakub Kruzik     if (Iidx[i] >= N) Iidx[i] = PETSC_MIN_INT;
178*e53e0a0dSJakub Kruzik   }
179*e53e0a0dSJakub Kruzik   for (i=ilo; i<ihi; i++) {
180*e53e0a0dSJakub Kruzik     ierr = MatSetValues(defl,1,&i,ncoeffs,Iidx,coeffs,INSERT_VALUES);CHKERRQ(ierr);
181*e53e0a0dSJakub Kruzik     for (j=0; j<ncoeffs; j++) {
182*e53e0a0dSJakub Kruzik       Iidx[j] += 2;
183*e53e0a0dSJakub Kruzik       if (Iidx[j] >= N) Iidx[j] = PETSC_MIN_INT;
184*e53e0a0dSJakub Kruzik     }
185*e53e0a0dSJakub Kruzik   }
186*e53e0a0dSJakub Kruzik 
187*e53e0a0dSJakub Kruzik   ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
188*e53e0a0dSJakub Kruzik   ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
189*e53e0a0dSJakub Kruzik 
190*e53e0a0dSJakub Kruzik   ierr = PetscFree(Iidx);CHKERRQ(ierr);
191*e53e0a0dSJakub Kruzik   *H = defl;
192*e53e0a0dSJakub Kruzik   PetscFunctionReturn(0);
193*e53e0a0dSJakub Kruzik }
194*e53e0a0dSJakub Kruzik 
195*e53e0a0dSJakub Kruzik PetscErrorCode PCDeflationGetSpaceHaar(PC pc,Mat *W,PetscInt size)
196*e53e0a0dSJakub Kruzik {
197*e53e0a0dSJakub Kruzik   PetscErrorCode ierr;
198*e53e0a0dSJakub Kruzik   Mat A,defl;
199*e53e0a0dSJakub Kruzik   PetscInt i,j,len,ilo,ihi,*Iidx,m,M;
200*e53e0a0dSJakub Kruzik   PetscReal *col,val;
201*e53e0a0dSJakub Kruzik 
202*e53e0a0dSJakub Kruzik   PetscFunctionBegin;
203*e53e0a0dSJakub Kruzik   /* Haar basis wavelet, level=size */
204*e53e0a0dSJakub Kruzik   len = pow(2,size);
205*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(len,&col);CHKERRQ(ierr);
206*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(len,&Iidx);CHKERRQ(ierr);
207*e53e0a0dSJakub Kruzik   val = 1./pow(2,size/2.);
208*e53e0a0dSJakub Kruzik   for (i=0; i<len; i++) col[i] = val;
209*e53e0a0dSJakub Kruzik 
210*e53e0a0dSJakub Kruzik   /* TODO pass A instead of PC? */
211*e53e0a0dSJakub Kruzik   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */
212*e53e0a0dSJakub Kruzik   ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
213*e53e0a0dSJakub Kruzik   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
214*e53e0a0dSJakub Kruzik   ierr = MatCreate(PetscObjectComm((PetscObject)A),&defl);CHKERRQ(ierr);
215*e53e0a0dSJakub Kruzik   ierr = MatSetSizes(defl,m,PETSC_DECIDE,M,ceil(M/(float)len));CHKERRQ(ierr);
216*e53e0a0dSJakub Kruzik   ierr = MatSetUp(defl);CHKERRQ(ierr);
217*e53e0a0dSJakub Kruzik   ierr = MatSeqAIJSetPreallocation(defl,size,NULL);CHKERRQ(ierr);
218*e53e0a0dSJakub Kruzik   ierr = MatMPIAIJSetPreallocation(defl,size,NULL,size,NULL);CHKERRQ(ierr);
219*e53e0a0dSJakub Kruzik   ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
220*e53e0a0dSJakub Kruzik 
221*e53e0a0dSJakub Kruzik   ierr = MatGetOwnershipRangeColumn(defl,&ilo,&ihi);CHKERRQ(ierr);
222*e53e0a0dSJakub Kruzik   for (i=0; i<len; i++) Iidx[i] = i+ilo*len;
223*e53e0a0dSJakub Kruzik   if (M%len && ihi == (int)ceil(M/(float)len)) ihi -= 1;
224*e53e0a0dSJakub Kruzik   for (i=ilo; i<ihi; i++) {
225*e53e0a0dSJakub Kruzik     ierr = MatSetValues(defl,len,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr);
226*e53e0a0dSJakub Kruzik     for (j=0; j<len; j++) Iidx[j] += len;
227*e53e0a0dSJakub Kruzik   }
228*e53e0a0dSJakub Kruzik   if (M%len && ihi+1 == ceil(M/(float)len)) {
229*e53e0a0dSJakub Kruzik     len = M%len;
230*e53e0a0dSJakub Kruzik     val = 1./pow(pow(2,len),0.5);
231*e53e0a0dSJakub Kruzik     for (i=0; i<len; i++) col[i] = val;
232*e53e0a0dSJakub Kruzik     ierr = MatSetValues(defl,len,Iidx,1,&ihi,col,INSERT_VALUES);CHKERRQ(ierr);
233*e53e0a0dSJakub Kruzik   }
234*e53e0a0dSJakub Kruzik 
235*e53e0a0dSJakub Kruzik   ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
236*e53e0a0dSJakub Kruzik   ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
237*e53e0a0dSJakub Kruzik 
238*e53e0a0dSJakub Kruzik   ierr = PetscFree(col);CHKERRQ(ierr);
239*e53e0a0dSJakub Kruzik   ierr = PetscFree(Iidx);CHKERRQ(ierr);
240*e53e0a0dSJakub Kruzik   *W = defl;
241*e53e0a0dSJakub Kruzik   PetscFunctionReturn(0);
242*e53e0a0dSJakub Kruzik }
243*e53e0a0dSJakub Kruzik 
244*e53e0a0dSJakub Kruzik PetscErrorCode PCDeflationGetSpaceJacketHaar(PC pc,Mat *W,PetscInt size)
245*e53e0a0dSJakub Kruzik {
246*e53e0a0dSJakub Kruzik   PetscErrorCode ierr;
247*e53e0a0dSJakub Kruzik   Mat A,*H,defl;
248*e53e0a0dSJakub Kruzik   PetscInt i,m,M,Mdefl,Ndefl;
249*e53e0a0dSJakub Kruzik   PetscBool jh;
250*e53e0a0dSJakub Kruzik   MPI_Comm comm;
251*e53e0a0dSJakub Kruzik 
252*e53e0a0dSJakub Kruzik   PetscFunctionBegin;
253*e53e0a0dSJakub Kruzik   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
254*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(size,&H);CHKERRQ(ierr);
255*e53e0a0dSJakub Kruzik   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */
256*e53e0a0dSJakub Kruzik   ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
257*e53e0a0dSJakub Kruzik   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
258*e53e0a0dSJakub Kruzik   Mdefl = M;
259*e53e0a0dSJakub Kruzik   Ndefl = M;
260*e53e0a0dSJakub Kruzik   for (i=0; i<size; i++) {
261*e53e0a0dSJakub Kruzik     if (Mdefl%2)  {
262*e53e0a0dSJakub Kruzik       jh=PETSC_TRUE;
263*e53e0a0dSJakub Kruzik       Mdefl = Mdefl/2 +1;
264*e53e0a0dSJakub Kruzik     } else {
265*e53e0a0dSJakub Kruzik       jh=PETSC_FALSE;
266*e53e0a0dSJakub Kruzik       Mdefl = Mdefl/2;
267*e53e0a0dSJakub Kruzik     }
268*e53e0a0dSJakub Kruzik     ierr = PCDeflationCreateSpaceJacketHaar(comm,PETSC_DECIDE,m,Mdefl,Ndefl,jh,&H[i]);CHKERRQ(ierr);
269*e53e0a0dSJakub Kruzik     ierr = MatGetLocalSize(H[i],&m,NULL);CHKERRQ(ierr);
270*e53e0a0dSJakub Kruzik     Ndefl = Mdefl;
271*e53e0a0dSJakub Kruzik   }
272*e53e0a0dSJakub Kruzik   //ierr = MatCreateProd(comm,size,H,&defl);CHKERRQ(ierr);
273*e53e0a0dSJakub Kruzik   //ierr = MatCreateComposite(comm,size,H,&defl);CHKERRQ(ierr);
274*e53e0a0dSJakub Kruzik   //ierr = MatCompositeSetType(defl,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
275*e53e0a0dSJakub Kruzik   /* TODO allow implicit */
276*e53e0a0dSJakub Kruzik   //ierr = MatCompositeMerge(defl);CHKERRQ(ierr);
277*e53e0a0dSJakub Kruzik   Mat newmat;
278*e53e0a0dSJakub Kruzik   defl = H[0];
279*e53e0a0dSJakub Kruzik   for (i=0; i<size-1; i++) {
280*e53e0a0dSJakub Kruzik     ierr = MatMatMult(H[i+1],defl,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr);
281*e53e0a0dSJakub Kruzik     ierr = MatDestroy(&defl);CHKERRQ(ierr);
282*e53e0a0dSJakub Kruzik     defl = newmat ;
283*e53e0a0dSJakub Kruzik   }
284*e53e0a0dSJakub Kruzik 
285*e53e0a0dSJakub Kruzik   ierr = MatTranspose(defl,MAT_INITIAL_MATRIX,W);CHKERRQ(ierr);
286*e53e0a0dSJakub Kruzik 
287*e53e0a0dSJakub Kruzik   ierr = MatDestroy(&defl);CHKERRQ(ierr);
288*e53e0a0dSJakub Kruzik   for (i=1; i<size; i++) {
289*e53e0a0dSJakub Kruzik     ierr = MatDestroy(&H[i]);CHKERRQ(ierr);
290*e53e0a0dSJakub Kruzik   }
291*e53e0a0dSJakub Kruzik   ierr = PetscFree(H);CHKERRQ(ierr);
292*e53e0a0dSJakub Kruzik   PetscFunctionReturn(0);
293*e53e0a0dSJakub Kruzik }
294*e53e0a0dSJakub Kruzik 
295*e53e0a0dSJakub Kruzik PetscErrorCode PCDeflationGetSpaceWave(PC pc,Mat *W,PetscInt size,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc)
296*e53e0a0dSJakub Kruzik {
297*e53e0a0dSJakub Kruzik   PetscErrorCode ierr;
298*e53e0a0dSJakub Kruzik   Mat A,*H,defl;
299*e53e0a0dSJakub Kruzik   PetscInt i,m,M,Mdefl,Ndefl;
300*e53e0a0dSJakub Kruzik   MPI_Comm comm;
301*e53e0a0dSJakub Kruzik 
302*e53e0a0dSJakub Kruzik   PetscFunctionBegin;
303*e53e0a0dSJakub Kruzik   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
304*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(size,&H);CHKERRQ(ierr);
305*e53e0a0dSJakub Kruzik   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */
306*e53e0a0dSJakub Kruzik   ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
307*e53e0a0dSJakub Kruzik   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
308*e53e0a0dSJakub Kruzik   Mdefl = M;
309*e53e0a0dSJakub Kruzik   Ndefl = M;
310*e53e0a0dSJakub Kruzik   for (i=0; i<size; i++) {
311*e53e0a0dSJakub Kruzik     if (Mdefl%2)  {
312*e53e0a0dSJakub Kruzik       if (trunc) {
313*e53e0a0dSJakub Kruzik         Mdefl = (PetscInt)PetscCeilReal(Mdefl/2.);
314*e53e0a0dSJakub Kruzik       } else {
315*e53e0a0dSJakub Kruzik         Mdefl = (PetscInt)PetscFloorReal((ncoeffs+Mdefl-1)/2.);
316*e53e0a0dSJakub Kruzik       }
317*e53e0a0dSJakub Kruzik     } else {
318*e53e0a0dSJakub Kruzik       Mdefl = Mdefl/2;
319*e53e0a0dSJakub Kruzik     }
320*e53e0a0dSJakub Kruzik     ierr = PCDeflationCreateSpaceWave(comm,PETSC_DECIDE,m,Mdefl,Ndefl,ncoeffs,coeffs,trunc,&H[i]);CHKERRQ(ierr);
321*e53e0a0dSJakub Kruzik     ierr = MatGetLocalSize(H[i],&m,NULL);CHKERRQ(ierr);
322*e53e0a0dSJakub Kruzik     Ndefl = Mdefl;
323*e53e0a0dSJakub Kruzik   }
324*e53e0a0dSJakub Kruzik   ierr = MatCreateComposite(comm,size,H,&defl);CHKERRQ(ierr);
325*e53e0a0dSJakub Kruzik   ierr = MatCompositeSetType(defl,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
326*e53e0a0dSJakub Kruzik   *W = defl;
327*e53e0a0dSJakub Kruzik 
328*e53e0a0dSJakub Kruzik   for (i=0; i<size; i++) {
329*e53e0a0dSJakub Kruzik     ierr = MatDestroy(&H[i]);CHKERRQ(ierr);
330*e53e0a0dSJakub Kruzik   }
331*e53e0a0dSJakub Kruzik   ierr = PetscFree(H);CHKERRQ(ierr);
332*e53e0a0dSJakub Kruzik   PetscFunctionReturn(0);
333*e53e0a0dSJakub Kruzik }
334*e53e0a0dSJakub Kruzik 
335*e53e0a0dSJakub Kruzik PetscErrorCode PCDeflationGetSpaceAggregation(PC pc,Mat *W)
336*e53e0a0dSJakub Kruzik {
337*e53e0a0dSJakub Kruzik   PetscErrorCode ierr;
338*e53e0a0dSJakub Kruzik   Mat A,defl;
339*e53e0a0dSJakub Kruzik   PetscInt i,ilo,ihi,*Iidx,m,M;
340*e53e0a0dSJakub Kruzik   PetscReal *col;
341*e53e0a0dSJakub Kruzik   MPI_Comm comm;
342*e53e0a0dSJakub Kruzik 
343*e53e0a0dSJakub Kruzik   PetscFunctionBegin;
344*e53e0a0dSJakub Kruzik   /* TODO pass A instead of PC? */
345*e53e0a0dSJakub Kruzik   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */
346*e53e0a0dSJakub Kruzik   ierr = MatGetOwnershipRangeColumn(A,&ilo,&ihi);CHKERRQ(ierr);
347*e53e0a0dSJakub Kruzik   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
348*e53e0a0dSJakub Kruzik   ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
349*e53e0a0dSJakub Kruzik   ierr = MPI_Comm_size(comm,&m);CHKERRQ(ierr);
350*e53e0a0dSJakub Kruzik   ierr = MatCreate(comm,&defl);CHKERRQ(ierr);
351*e53e0a0dSJakub Kruzik   ierr = MatSetSizes(defl,ihi-ilo,1,M,m);CHKERRQ(ierr);
352*e53e0a0dSJakub Kruzik   ierr = MatSetUp(defl);CHKERRQ(ierr);
353*e53e0a0dSJakub Kruzik   ierr = MatSeqAIJSetPreallocation(defl,1,NULL);CHKERRQ(ierr);
354*e53e0a0dSJakub Kruzik   ierr = MatMPIAIJSetPreallocation(defl,1,NULL,0,NULL);CHKERRQ(ierr);
355*e53e0a0dSJakub Kruzik   ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
356*e53e0a0dSJakub Kruzik   ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
357*e53e0a0dSJakub Kruzik 
358*e53e0a0dSJakub Kruzik 
359*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(ihi-ilo,&col);CHKERRQ(ierr);
360*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(ihi-ilo,&Iidx);CHKERRQ(ierr);
361*e53e0a0dSJakub Kruzik   for (i=ilo; i<ihi; i++) {
362*e53e0a0dSJakub Kruzik     Iidx[i-ilo] = i;
363*e53e0a0dSJakub Kruzik     col[i-ilo] = 1;
364*e53e0a0dSJakub Kruzik   }
365*e53e0a0dSJakub Kruzik   ierr = MPI_Comm_rank(comm,&i);CHKERRQ(ierr);
366*e53e0a0dSJakub Kruzik   ierr = MatSetValues(defl,ihi-ilo,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr);
367*e53e0a0dSJakub Kruzik 
368*e53e0a0dSJakub Kruzik   ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
369*e53e0a0dSJakub Kruzik   ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
370*e53e0a0dSJakub Kruzik 
371*e53e0a0dSJakub Kruzik   ierr = PetscFree(col);CHKERRQ(ierr);
372*e53e0a0dSJakub Kruzik   ierr = PetscFree(Iidx);CHKERRQ(ierr);
373*e53e0a0dSJakub Kruzik   *W = defl;
374*e53e0a0dSJakub Kruzik   PetscFunctionReturn(0);
375*e53e0a0dSJakub Kruzik }
376*e53e0a0dSJakub Kruzik 
377*e53e0a0dSJakub Kruzik PetscErrorCode PCDeflationGetSpaceSLEPc(PC pc,Mat *W,PetscInt size,PetscBool cheapCP)
378*e53e0a0dSJakub Kruzik {
379*e53e0a0dSJakub Kruzik #if defined(HAVE_SLEPC)
380*e53e0a0dSJakub Kruzik   PetscErrorCode ierr;
381*e53e0a0dSJakub Kruzik   PC_Deflation *def = (PC_Deflation*)pc->data;
382*e53e0a0dSJakub Kruzik   Mat A,defl;
383*e53e0a0dSJakub Kruzik   Vec vec;
384*e53e0a0dSJakub Kruzik   EPS eps;
385*e53e0a0dSJakub Kruzik   PetscScalar *data,*dataScaled,eigval;
386*e53e0a0dSJakub Kruzik   PetscInt i,nconv,m,M,n=PETSC_DECIDE;
387*e53e0a0dSJakub Kruzik   PetscBool slepcinit;
388*e53e0a0dSJakub Kruzik   MPI_Comm comm;
389*e53e0a0dSJakub Kruzik 
390*e53e0a0dSJakub Kruzik   PetscFunctionBegin;
391*e53e0a0dSJakub Kruzik   ierr = SlepcInitialized(&slepcinit);CHKERRQ(ierr);
392*e53e0a0dSJakub Kruzik   if (!slepcinit) {
393*e53e0a0dSJakub Kruzik     ierr = SlepcInitialize(NULL,NULL,(char*)0,(char*)0);CHKERRQ(ierr);
394*e53e0a0dSJakub Kruzik     slepcinit = PETSC_TRUE;
395*e53e0a0dSJakub Kruzik   }
396*e53e0a0dSJakub Kruzik   ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */
397*e53e0a0dSJakub Kruzik   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
398*e53e0a0dSJakub Kruzik   ierr = EPSCreate(comm,&eps);CHKERRQ(ierr);
399*e53e0a0dSJakub Kruzik   ierr = EPSSetOperators(eps,A,NULL);CHKERRQ(ierr);
400*e53e0a0dSJakub Kruzik   ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr); /* Implemented only for def */
401*e53e0a0dSJakub Kruzik   ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr);
402*e53e0a0dSJakub Kruzik   ierr = EPSSetDimensions(eps,size,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
403*e53e0a0dSJakub Kruzik   ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
404*e53e0a0dSJakub Kruzik 
405*e53e0a0dSJakub Kruzik   ierr = EPSSolve(eps);CHKERRQ(ierr);
406*e53e0a0dSJakub Kruzik   ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr);
407*e53e0a0dSJakub Kruzik   if (nconv < size) SETERRQ2(comm,PETSC_ERR_CONV_FAILED,"SLEPc: Number of converged eigenpairs (%d) is less than requested (%d)",nconv,size);
408*e53e0a0dSJakub Kruzik   ierr = MatCreateVecs(A,NULL,&vec);CHKERRQ(ierr);
409*e53e0a0dSJakub Kruzik   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
410*e53e0a0dSJakub Kruzik   ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
411*e53e0a0dSJakub Kruzik   ierr = PetscSplitOwnership(comm,&n,&size);CHKERRQ(ierr);
412*e53e0a0dSJakub Kruzik   ierr = PetscMalloc1(m*size,&data);CHKERRQ(ierr);
413*e53e0a0dSJakub Kruzik   /* TODO check that eigenvalue is not 0 -> vec is not in Ker A */
414*e53e0a0dSJakub Kruzik   for (i=0; i<size; i++) {
415*e53e0a0dSJakub Kruzik     ierr = VecPlaceArray(vec,&data[i*m]);CHKERRQ(ierr);
416*e53e0a0dSJakub Kruzik     ierr = EPSGetEigenvector(eps,i,vec,NULL);CHKERRQ(ierr);
417*e53e0a0dSJakub Kruzik     ierr = VecResetArray(vec);CHKERRQ(ierr);
418*e53e0a0dSJakub Kruzik   }
419*e53e0a0dSJakub Kruzik   ierr = MatCreateDense(comm,m,n,M,size,data,&defl);CHKERRQ(ierr);
420*e53e0a0dSJakub Kruzik 
421*e53e0a0dSJakub Kruzik   if (cheapCP) {
422*e53e0a0dSJakub Kruzik     ierr = PetscMalloc1(m*size,&dataScaled);CHKERRQ(ierr);
423*e53e0a0dSJakub Kruzik     for (i=0; i<size; i++) {
424*e53e0a0dSJakub Kruzik         ierr = VecPlaceArray(vec,&dataScaled[i*m]);CHKERRQ(ierr);
425*e53e0a0dSJakub Kruzik         ierr = EPSGetEigenpair(eps,i,&eigval,NULL,vec,NULL);CHKERRQ(ierr);
426*e53e0a0dSJakub Kruzik         ierr = VecScale(vec,eigval);CHKERRQ(ierr);
427*e53e0a0dSJakub Kruzik         ierr = VecResetArray(vec);CHKERRQ(ierr);
428*e53e0a0dSJakub Kruzik     }
429*e53e0a0dSJakub Kruzik     ierr = MatCreateDense(comm,m,n,M,size,dataScaled,&def->AW);CHKERRQ(ierr);
430*e53e0a0dSJakub Kruzik   }
431*e53e0a0dSJakub Kruzik 
432*e53e0a0dSJakub Kruzik   //ierr = EPSGetBV(eps,&bv);CHKERRQ(ierr);
433*e53e0a0dSJakub Kruzik   //ierr = BVCreateMat(bv,&defl);CHKERRQ(ierr);
434*e53e0a0dSJakub Kruzik   *W = defl;
435*e53e0a0dSJakub Kruzik 
436*e53e0a0dSJakub Kruzik   ierr = EPSDestroy(&eps);CHKERRQ(ierr);
437*e53e0a0dSJakub Kruzik   if (slepcinit) ierr = SlepcFinalize();CHKERRQ(ierr);
438*e53e0a0dSJakub Kruzik   PetscFunctionReturn(0);
439*e53e0a0dSJakub Kruzik #else
440*e53e0a0dSJakub Kruzik   SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_CONV_FAILED,"Not compiled with SLEPc support (call make HAVE_SLEPC)");
441*e53e0a0dSJakub Kruzik #endif
442*e53e0a0dSJakub Kruzik }
443*e53e0a0dSJakub Kruzik 
444*e53e0a0dSJakub Kruzik PetscErrorCode PCDeflationComputeSpace(PC pc)
445*e53e0a0dSJakub Kruzik {
446*e53e0a0dSJakub Kruzik   PetscErrorCode ierr;
447*e53e0a0dSJakub Kruzik   Mat defl;
448*e53e0a0dSJakub Kruzik   PetscBool transp=PETSC_TRUE;
449*e53e0a0dSJakub Kruzik   PC_Deflation *def = (PC_Deflation*)pc->data;
450*e53e0a0dSJakub Kruzik 
451*e53e0a0dSJakub Kruzik   /* TODO valid header */
452*e53e0a0dSJakub Kruzik   PetscFunctionBegin;
453*e53e0a0dSJakub Kruzik   if (def->spacesize < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PC_DEFLATION Space size specified: %d",def->spacesize);
454*e53e0a0dSJakub Kruzik   switch (def->spacetype) {
455*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_HAAR:
456*e53e0a0dSJakub Kruzik       transp = PETSC_FALSE;
457*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceHaar(pc,&defl,def->spacesize);CHKERRQ(ierr);break;
458*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_JACKET_HAAR:
459*e53e0a0dSJakub Kruzik       transp = PETSC_FALSE;
460*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceJacketHaar(pc,&defl,def->spacesize);CHKERRQ(ierr);break;
461*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_DB2:
462*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,2,db2,!def->extendsp);CHKERRQ(ierr);break;
463*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_DB4:
464*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,4,db4,!def->extendsp);CHKERRQ(ierr);break;
465*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_DB8:
466*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,8,db8,!def->extendsp);CHKERRQ(ierr);break;
467*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_DB16:
468*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,16,db16,!def->extendsp);CHKERRQ(ierr);break;
469*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_BIORTH22:
470*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,6,biorth22,!def->extendsp);CHKERRQ(ierr);break;
471*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_MEYER:
472*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,62,meyer,!def->extendsp);CHKERRQ(ierr);break;
473*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_AGGREGATION:
474*e53e0a0dSJakub Kruzik       transp = PETSC_FALSE;
475*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceAggregation(pc,&defl);CHKERRQ(ierr);break;
476*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_SLEPC:
477*e53e0a0dSJakub Kruzik       transp = PETSC_FALSE;
478*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceSLEPc(pc,&defl,def->spacesize,PETSC_FALSE);CHKERRQ(ierr);break;
479*e53e0a0dSJakub Kruzik     case PC_DEFLATION_SPACE_SLEPC_CHEAP:
480*e53e0a0dSJakub Kruzik       transp = PETSC_FALSE;
481*e53e0a0dSJakub Kruzik       ierr = PCDeflationGetSpaceSLEPc(pc,&defl,def->spacesize,PETSC_TRUE);CHKERRQ(ierr);break;
482*e53e0a0dSJakub Kruzik     default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PC_DEFLATION Space Type specified");
483*e53e0a0dSJakub Kruzik   }
484*e53e0a0dSJakub Kruzik 
485*e53e0a0dSJakub Kruzik   ierr = PCDeflationSetSpace(pc,defl,transp);CHKERRQ(ierr);
486*e53e0a0dSJakub Kruzik   PetscFunctionReturn(0);
487*e53e0a0dSJakub Kruzik }
488*e53e0a0dSJakub Kruzik 
489