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