Grid 0.7.0
Aggregates.h
Go to the documentation of this file.
1/*************************************************************************************
2
3 Grid physics library, www.github.com/paboyle/Grid
4
5 Source file: ./lib/algorithms/Aggregates.h
6
7 Copyright (C) 2015
8
9Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
10Author: Peter Boyle <paboyle@ph.ed.ac.uk>
11Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
12Author: paboyle <paboyle@ph.ed.ac.uk>
13
14 This program is free software; you can redistribute it and/or modify
15 it under the terms of the GNU General Public License as published by
16 the Free Software Foundation; either version 2 of the License, or
17 (at your option) any later version.
18
19 This program is distributed in the hope that it will be useful,
20 but WITHOUT ANY WARRANTY; without even the implied warranty of
21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 GNU General Public License for more details.
23
24 You should have received a copy of the GNU General Public License along
25 with this program; if not, write to the Free Software Foundation, Inc.,
26 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27
28 See the full license in the file "LICENSE" in the top level distribution directory
29*************************************************************************************/
30/* END LEGAL */
31#pragma once
32
33#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
34
35NAMESPACE_BEGIN(Grid);
36
37inline RealD AggregatePowerLaw(RealD x)
38{
39 // return std::pow(x,-4);
40 // return std::pow(x,-3);
41 return std::pow(x,-5);
42}
43
44template<class Fobj,class CComplex,int nbasis>
45class Aggregation {
46public:
47 constexpr int Nbasis(void) { return nbasis; };
48
49 typedef iVector<CComplex,nbasis > siteVector;
50 typedef Lattice<siteVector> CoarseVector;
51 typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
52
53 typedef Lattice< CComplex > CoarseScalar; // used for inner products on fine field
54 typedef Lattice<Fobj > FineField;
55
56 GridBase *CoarseGrid;
57 GridBase *FineGrid;
58 std::vector<Lattice<Fobj> > subspace;
59 int checkerboard;
60 int Checkerboard(void){return checkerboard;}
61 Aggregation(GridBase *_CoarseGrid,GridBase *_FineGrid,int _checkerboard) :
62 CoarseGrid(_CoarseGrid),
63 FineGrid(_FineGrid),
64 subspace(nbasis,_FineGrid),
65 checkerboard(_checkerboard)
66 {
67 };
68
69
70 void Orthogonalise(void){
71 CoarseScalar InnerProd(CoarseGrid);
72 // std::cout << GridLogMessage <<" Block Gramm-Schmidt pass 1"<<std::endl;
73 blockOrthogonalise(InnerProd,subspace);
74 }
75 void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
76 blockProject(CoarseVec,FineVec,subspace);
77 }
78 void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
79 FineVec.Checkerboard() = subspace[0].Checkerboard();
80 blockPromote(CoarseVec,FineVec,subspace);
81 }
82
83 virtual void CreateSubspaceRandom(GridParallelRNG &RNG) {
84 int nn=nbasis;
85 RealD scale;
86 FineField noise(FineGrid);
87 for(int b=0;b<nn;b++){
88 subspace[b] = Zero();
89 gaussian(RNG,noise);
90 scale = std::pow(norm2(noise),-0.5);
91 noise=noise*scale;
92 subspace[b] = noise;
93 }
94 }
95 virtual void CreateSubspace(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis)
96 {
97
98 RealD scale;
99
100 ConjugateGradient<FineField> CG(1.0e-3,400,false);
101 FineField noise(FineGrid);
102 FineField Mn(FineGrid);
103
104 for(int b=0;b<nn;b++){
105
106 subspace[b] = Zero();
107 gaussian(RNG,noise);
108 scale = std::pow(norm2(noise),-0.5);
109 noise=noise*scale;
110
111 hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise ["<<b<<"] <n|MdagM|n> "<<norm2(Mn)<<std::endl;
112
113 for(int i=0;i<4;i++){
114
115 CG(hermop,noise,subspace[b]);
116
117 noise = subspace[b];
118 scale = std::pow(norm2(noise),-0.5);
119 noise=noise*scale;
120
121 }
122
123 hermop.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(Mn)<<std::endl;
124 subspace[b] = noise;
125
126 }
127 }
128
129 virtual void CreateSubspaceGCR(GridParallelRNG &RNG,LinearOperatorBase<FineField> &DiracOp,int nn=nbasis)
130 {
131 RealD scale;
132
133 TrivialPrecon<FineField> simple_fine;
134 PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,30,DiracOp,simple_fine,12,12);
135 FineField noise(FineGrid);
136 FineField src(FineGrid);
137 FineField guess(FineGrid);
138 FineField Mn(FineGrid);
139
140 for(int b=0;b<nn;b++){
141
142 subspace[b] = Zero();
143 gaussian(RNG,noise);
144 scale = std::pow(norm2(noise),-0.5);
145 noise=noise*scale;
146
147 DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "noise ["<<b<<"] <n|Op|n> "<<innerProduct(noise,Mn)<<std::endl;
148
149 for(int i=0;i<2;i++){
150 // void operator() (const Field &src, Field &psi){
151#if 1
152 std::cout << GridLogMessage << " inverting on noise "<<std::endl;
153 src = noise;
154 guess=Zero();
155 GCR(src,guess);
156 subspace[b] = guess;
157#else
158 std::cout << GridLogMessage << " inverting on zero "<<std::endl;
159 src=Zero();
160 guess = noise;
161 GCR(src,guess);
162 subspace[b] = guess;
163#endif
164 noise = subspace[b];
165 scale = std::pow(norm2(noise),-0.5);
166 noise=noise*scale;
167
168 }
169
170 DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|Op|f> "<<innerProduct(noise,Mn)<<std::endl;
171 subspace[b] = noise;
172
173 }
174 }
175
177 // World of possibilities here. But have tried quite a lot of experiments (250+ jobs run on Summit)
178 // and this is the best I found
180
181 virtual void CreateSubspaceChebyshev(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,
182 int nn,
183 double hi,
184 double lo,
185 int orderfilter,
186 int ordermin,
187 int orderstep,
188 double filterlo
189 ) {
190
191 RealD scale;
192
193 FineField noise(FineGrid);
194 FineField Mn(FineGrid);
195 FineField tmp(FineGrid);
196
197 // New normalised noise
198 gaussian(RNG,noise);
199 scale = std::pow(norm2(noise),-0.5);
200 noise=noise*scale;
201
202 std::cout << GridLogMessage<<" Chebyshev subspace pass-1 : ord "<<orderfilter<<" ["<<lo<<","<<hi<<"]"<<std::endl;
203 std::cout << GridLogMessage<<" Chebyshev subspace pass-2 : nbasis"<<nn<<" min "
204 <<ordermin<<" step "<<orderstep
205 <<" lo"<<filterlo<<std::endl;
206
207 // Initial matrix element
208 hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
209
210 int b =0;
211 {
212 ComplexD ip;
213 // Filter
214 Chebyshev<FineField> Cheb(lo,hi,orderfilter);
215 Cheb(hermop,noise,Mn);
216 // normalise
217 scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
218 subspace[b] = Mn;
219
220 hermop.Op(Mn,tmp);
221 ip= innerProduct(Mn,tmp);
222 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|Op|n> "<<norm2(tmp)<<" "<<ip<<std::endl;
223
224 hermop.AdjOp(Mn,tmp);
225 ip = innerProduct(Mn,tmp);
226 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|AdjOp|n> "<<norm2(tmp)<<" "<<ip<<std::endl;
227 b++;
228 }
229
230 // Generate a full sequence of Chebyshevs
231 {
232 lo=filterlo;
233 noise=Mn;
234
235 FineField T0(FineGrid); T0 = noise;
236 FineField T1(FineGrid);
237 FineField T2(FineGrid);
238 FineField y(FineGrid);
239
240 FineField *Tnm = &T0;
241 FineField *Tn = &T1;
242 FineField *Tnp = &T2;
243
244 // Tn=T1 = (xscale M + mscale)in
245 RealD xscale = 2.0/(hi-lo);
246 RealD mscale = -(hi+lo)/(hi-lo);
247 hermop.HermOp(T0,y);
248 T1=y*xscale+noise*mscale;
249
250 for(int n=2;n<=ordermin+orderstep*(nn-2);n++){
251
252 hermop.HermOp(*Tn,y);
253
254 autoView( y_v , y, AcceleratorWrite);
255 autoView( Tn_v , (*Tn), AcceleratorWrite);
256 autoView( Tnp_v , (*Tnp), AcceleratorWrite);
257 autoView( Tnm_v , (*Tnm), AcceleratorWrite);
258 const int Nsimd = CComplex::Nsimd();
259 accelerator_for(ss, FineGrid->oSites(), Nsimd, {
260 coalescedWrite(y_v[ss],xscale*y_v(ss)+mscale*Tn_v(ss));
261 coalescedWrite(Tnp_v[ss],2.0*y_v(ss)-Tnm_v(ss));
262 });
263
264 // Possible more fine grained control is needed than a linear sweep,
265 // but huge productivity gain if this is simple algorithm and not a tunable
266 int m =1;
267 if ( n>=ordermin ) m=n-ordermin;
268 if ( (m%orderstep)==0 ) {
269 Mn=*Tnp;
270 scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
271 subspace[b] = Mn;
272
273
274 ComplexD ip;
275
276 hermop.Op(Mn,tmp);
277 ip= innerProduct(Mn,tmp);
278 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|Op|n> "<<norm2(tmp)<<" "<<ip<<std::endl;
279
280 hermop.AdjOp(Mn,tmp);
281 ip = innerProduct(Mn,tmp);
282 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|AdjOp|n> "<<norm2(tmp)<<" "<<ip<<std::endl;
283
284 b++;
285 }
286
287 // Cycle pointers to avoid copies
288 FineField *swizzle = Tnm;
289 Tnm =Tn;
290 Tn =Tnp;
291 Tnp =swizzle;
292
293 }
294 }
295 assert(b==nn);
296 }
297
298
299 virtual void CreateSubspacePolyCheby(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,
300 int nn,
301 double hi,
302 double lo1,
303 int orderfilter,
304 double lo2,
305 int orderstep)
306 {
307 RealD scale;
308
309 FineField noise(FineGrid);
310 FineField Mn(FineGrid);
311 FineField tmp(FineGrid);
312
313 // New normalised noise
314 gaussian(RNG,noise);
315 scale = std::pow(norm2(noise),-0.5);
316 noise=noise*scale;
317
318 std::cout << GridLogMessage<<" CreateSubspacePolyCheby "<<std::endl;
319 // Initial matrix element
320 hermop.Op(noise,Mn);
321 std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
322
323 int b =0;
324 {
325 // Filter
326 std::cout << GridLogMessage << "Cheby "<<lo1<<","<<hi<<" "<<orderstep<<std::endl;
327 Chebyshev<FineField> Cheb(lo1,hi,orderfilter);
328 Cheb(hermop,noise,Mn);
329 // normalise
330 scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
331 subspace[b] = Mn;
332 hermop.Op(Mn,tmp);
333 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
334 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|n> "<<norm2(Mn)<<std::endl;
335 }
336
337 // Generate a full sequence of Chebyshevs
338 for(int n=1;n<nn;n++){
339 std::cout << GridLogMessage << "Cheby "<<lo2<<","<<hi<<" "<<orderstep<<std::endl;
340 Chebyshev<FineField> Cheb(lo2,hi,orderstep);
341 Cheb(hermop,subspace[n-1],Mn);
342
343 for(int m=0;m<n;m++){
344 ComplexD c = innerProduct(subspace[m],Mn);
345 Mn = Mn - c*subspace[m];
346 }
347
348 // normalise
349 scale = std::pow(norm2(Mn),-0.5);
350 Mn=Mn*scale;
351
352 subspace[n]=Mn;
353
354 hermop.Op(Mn,tmp);
355 std::cout<<GridLogMessage << "filt ["<<n<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
356 std::cout<<GridLogMessage << "filt ["<<n<<"] <n|n> "<<norm2(Mn)<<std::endl;
357
358 }
359 }
360
361 virtual void CreateSubspaceChebyshev(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,
362 int nn,
363 double hi,
364 double lo,
365 int orderfilter
366 ) {
367
368 RealD scale;
369
370 FineField noise(FineGrid);
371 FineField Mn(FineGrid);
372 FineField tmp(FineGrid);
373
374 // New normalised noise
375 std::cout << GridLogMessage<<" Chebyshev subspace pure noise : ord "<<orderfilter<<" ["<<lo<<","<<hi<<"]"<<std::endl;
376 std::cout << GridLogMessage<<" Chebyshev subspace pure noise : nbasis "<<nn<<std::endl;
377
378
379 for(int b =0;b<nbasis;b++)
380 {
381 gaussian(RNG,noise);
382 scale = std::pow(norm2(noise),-0.5);
383 noise=noise*scale;
384
385 // Initial matrix element
386 hermop.Op(noise,Mn);
387 if(b==0) std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
388
389 // Filter
390 Chebyshev<FineField> Cheb(lo,hi,orderfilter);
391 Cheb(hermop,noise,Mn);
392 scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
393
394 // Refine
395 Chebyshev<FineField> PowerLaw(lo,hi,1000,AggregatePowerLaw);
396 noise = Mn;
397 PowerLaw(hermop,noise,Mn);
398 scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
399
400 // normalise
401 subspace[b] = Mn;
402 hermop.Op(Mn,tmp);
403 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
404 }
405
406 }
407
408 virtual void CreateSubspaceChebyshevPowerLaw(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,
409 int nn,
410 double hi,
411 int orderfilter
412 ) {
413
414 RealD scale;
415
416 FineField noise(FineGrid);
417 FineField Mn(FineGrid);
418 FineField tmp(FineGrid);
419
420 // New normalised noise
421 std::cout << GridLogMessage<<" Chebyshev subspace pure noise : ord "<<orderfilter<<" [0,"<<hi<<"]"<<std::endl;
422 std::cout << GridLogMessage<<" Chebyshev subspace pure noise : nbasis "<<nn<<std::endl;
423
424 for(int b =0;b<nbasis;b++)
425 {
426 gaussian(RNG,noise);
427 scale = std::pow(norm2(noise),-0.5);
428 noise=noise*scale;
429
430 // Initial matrix element
431 hermop.Op(noise,Mn);
432 if(b==0) std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
433 // Filter
434 Chebyshev<FineField> Cheb(0.0,hi,orderfilter,AggregatePowerLaw);
435 Cheb(hermop,noise,Mn);
436 // normalise
437 scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
438 subspace[b] = Mn;
439 hermop.Op(Mn,tmp);
440 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
441 }
442
443 }
444 virtual void CreateSubspaceChebyshevNew(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,
445 double hi
446 ) {
447
448 RealD scale;
449
450 FineField noise(FineGrid);
451 FineField Mn(FineGrid);
452 FineField tmp(FineGrid);
453
454 // New normalised noise
455 for(int b =0;b<nbasis;b++)
456 {
457 gaussian(RNG,noise);
458 scale = std::pow(norm2(noise),-0.5);
459 noise=noise*scale;
460
461 // Initial matrix element
462 hermop.Op(noise,Mn);
463 if(b==0) std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
464 // Filter
465 //#opt2(x) = acheb(x,3,90,300)* acheb(x,1,90,50) * acheb(x,0.5,90,200) * acheb(x,0.05,90,400) * acheb(x,0.01,90,1500)
466 /*266
467 Chebyshev<FineField> Cheb1(3.0,hi,300);
468 Chebyshev<FineField> Cheb2(1.0,hi,50);
469 Chebyshev<FineField> Cheb3(0.5,hi,300);
470 Chebyshev<FineField> Cheb4(0.05,hi,500);
471 Chebyshev<FineField> Cheb5(0.01,hi,2000);
472 */
473 /* 242 */
474 /*
475 Chebyshev<FineField> Cheb3(0.1,hi,300);
476 Chebyshev<FineField> Cheb2(0.02,hi,1000);
477 Chebyshev<FineField> Cheb1(0.003,hi,2000);
478 8?
479 */
480 /* How many??
481 */
482 Chebyshev<FineField> Cheb2(0.001,hi,2500); // 169 iters on HDCG after refine
483 Chebyshev<FineField> Cheb1(0.02,hi,600);
484
485 // Chebyshev<FineField> Cheb2(0.001,hi,1500);
486 // Chebyshev<FineField> Cheb1(0.02,hi,600);
487 Cheb1(hermop,noise,Mn); scale = std::pow(norm2(Mn),-0.5); noise=Mn*scale;
488 hermop.Op(noise,tmp); std::cout<<GridLogMessage << "Cheb1 <n|MdagM|n> "<<norm2(tmp)<<std::endl;
489 Cheb2(hermop,noise,Mn); scale = std::pow(norm2(Mn),-0.5); noise=Mn*scale;
490 hermop.Op(noise,tmp); std::cout<<GridLogMessage << "Cheb2 <n|MdagM|n> "<<norm2(tmp)<<std::endl;
491 // Cheb3(hermop,noise,Mn); scale = std::pow(norm2(Mn),-0.5); noise=Mn*scale;
492 // hermop.Op(noise,tmp); std::cout<<GridLogMessage << "Cheb3 <n|MdagM|n> "<<norm2(tmp)<<std::endl;
493 // Cheb4(hermop,noise,Mn); scale = std::pow(norm2(Mn),-0.5); noise=Mn*scale;
494 // hermop.Op(noise,tmp); std::cout<<GridLogMessage << "Cheb4 <n|MdagM|n> "<<norm2(tmp)<<std::endl;
495 // Cheb5(hermop,noise,Mn); scale = std::pow(norm2(Mn),-0.5); noise=Mn*scale;
496 // hermop.Op(noise,tmp); std::cout<<GridLogMessage << "Cheb5 <n|MdagM|n> "<<norm2(tmp)<<std::endl;
497 subspace[b] = noise;
498 hermop.Op(subspace[b],tmp);
499 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<< " norm " << norm2(noise)<<std::endl;
500 }
501
502 }
503
504 virtual void CreateSubspaceMultishift(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,
505 double Lo,double tol,int maxit)
506 {
507
508 RealD scale;
509
510 FineField noise(FineGrid);
511 FineField Mn(FineGrid);
512 FineField tmp(FineGrid);
513
514 // New normalised noise
515 std::cout << GridLogMessage<<" Multishift subspace : Lo "<<Lo<<std::endl;
516
517 // Filter
518 // [ 1/6(x+Lo) - 1/2(x+2Lo) + 1/2(x+3Lo) -1/6(x+4Lo) = Lo^3 /[ (x+1Lo)(x+2Lo)(x+3Lo)(x+4Lo) ]
519 //
520 // 1/(x+Lo) - 1/(x+2 Lo)
521 double epsilon = Lo/3;
522 std::vector<RealD> alpha({1.0/6.0,-1.0/2.0,1.0/2.0,-1.0/6.0});
523 std::vector<RealD> shifts({Lo,Lo+epsilon,Lo+2*epsilon,Lo+3*epsilon});
524 std::vector<RealD> tols({tol,tol,tol,tol});
525 std::cout << "sizes "<<alpha.size()<<" "<<shifts.size()<<" "<<tols.size()<<std::endl;
526
527 MultiShiftFunction msf(4,0.0,95.0);
528 std::cout << "msf constructed "<<std::endl;
529 msf.poles=shifts;
530 msf.residues=alpha;
531 msf.tolerances=tols;
532 msf.norm=0.0;
533 msf.order=alpha.size();
534 ConjugateGradientMultiShift<FineField> MSCG(maxit,msf);
535
536 for(int b =0;b<nbasis;b++)
537 {
538 gaussian(RNG,noise);
539 scale = std::pow(norm2(noise),-0.5);
540 noise=noise*scale;
541
542 // Initial matrix element
543 hermop.Op(noise,Mn);
544 if(b==0) std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
545
546 MSCG(hermop,noise,Mn);
547 scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
548 subspace[b] = Mn;
549 hermop.Op(Mn,tmp);
550 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
551
552 }
553
554 }
555 virtual void RefineSubspace(LinearOperatorBase<FineField> &hermop,
556 double Lo,double tol,int maxit)
557 {
558 FineField tmp(FineGrid);
559 for(int b =0;b<nbasis;b++)
560 {
561 ConjugateGradient<FineField> CGsloppy(tol,maxit,false);
562 ShiftedHermOpLinearOperator<FineField> ShiftedFineHermOp(hermop,Lo);
563 tmp=Zero();
564 CGsloppy(hermop,subspace[b],tmp);
565 RealD scale = std::pow(norm2(tmp),-0.5); tmp=tmp*scale;
566 subspace[b]=tmp;
567 hermop.Op(subspace[b],tmp);
568 std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
569 }
570 }
571 virtual void RefineSubspaceHDCG(LinearOperatorBase<FineField> &hermop,
572 TwoLevelADEF2mrhs<FineField,CoarseVector> & theHDCG,
573 int nrhs)
574 {
575 std::vector<FineField> src_mrhs(nrhs,FineGrid);
576 std::vector<FineField> res_mrhs(nrhs,FineGrid);
577 FineField tmp(FineGrid);
578 for(int b =0;b<nbasis;b+=nrhs)
579 {
580 tmp = subspace[b];
581 RealD scale = std::pow(norm2(tmp),-0.5); tmp=tmp*scale;
582 subspace[b] =tmp;
583 hermop.Op(subspace[b],tmp);
584 std::cout<<GridLogMessage << "before filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
585
586 for(int r=0;r<MIN(nbasis-b,nrhs);r++){
587 src_mrhs[r] = subspace[b+r];
588 }
589 for(int r=0;r<nrhs;r++){
590 res_mrhs[r] = Zero();
591 }
592 theHDCG(src_mrhs,res_mrhs);
593
594 for(int r=0;r<MIN(nbasis-b,nrhs);r++){
595 tmp = res_mrhs[r];
596 RealD scale = std::pow(norm2(tmp),-0.5); tmp=tmp*scale;
597 subspace[b+r]=tmp;
598 }
599 hermop.Op(subspace[b],tmp);
600 std::cout<<GridLogMessage << "after filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
601 }
602 }
603
604
605
606};
607NAMESPACE_END(Grid);
608
accelerator_for
#define accelerator_for(iterator, num, nsimd,...)
Definition Accelerator.h:609
AggregatePowerLaw
RealD AggregatePowerLaw(RealD x)
Definition Aggregates.h:37
innerProduct
ComplexD innerProduct(const Lattice< vobj > &left, const Lattice< vobj > &right)
Definition Lattice_reduction.h:290
norm2
RealD norm2(const Lattice< vobj > &arg)
Definition Lattice_reduction.h:202
gaussian
void gaussian(GridParallelRNG &rng, Lattice< vobj > &l)
Definition Lattice_rng.h:533
blockOrthogonalise
void blockOrthogonalise(Lattice< CComplex > &ip, std::vector< Lattice< vobj > > &Basis)
Definition Lattice_transfer.h:599
blockPromote
void blockPromote(const Lattice< iVector< CComplex, nbasis > > &coarseData, Lattice< vobj > &fineData, const VLattice &Basis)
Definition Lattice_transfer.h:662
blockProject
void blockProject(Lattice< iVector< CComplex, nbasis > > &coarseData, const Lattice< vobj > &fineData, const VLattice &Basis)
Definition Lattice_transfer.h:267
autoView
#define autoView(l_v, l, mode)
Definition Lattice_view.h:119
GridLogMessage
GridLogger GridLogMessage(1, "Message", GridLogColours, "NORMAL")
AcceleratorWrite
@ AcceleratorWrite
Definition MemoryManager.h:67
NAMESPACE_BEGIN
#define NAMESPACE_BEGIN(A)
Definition Namespace.h:35
NAMESPACE_END
#define NAMESPACE_END(A)
Definition Namespace.h:36
ComplexD
std::complex< RealD > ComplexD
Definition Simd.h:79
RealD
double RealD
Definition Simd.h:61
T2
#define T2
Definition StaggeredKernelsAsm.h:95
T1
#define T1
Definition StaggeredKernelsAsm.h:94
T0
#define T0
Definition StaggeredKernelsAsm.h:93
MIN
#define MIN(a, b)
Definition Zolotarev.cc:23
Aggregation::subspace
std::vector< Lattice< Fobj > > subspace
Definition Aggregates.h:58
Aggregation::Orthogonalise
void Orthogonalise(void)
Definition Aggregates.h:70
Aggregation::CoarseVector
Lattice< siteVector > CoarseVector
Definition Aggregates.h:50
Aggregation::FineGrid
GridBase * FineGrid
Definition Aggregates.h:57
Aggregation::CreateSubspaceChebyshev
virtual void CreateSubspaceChebyshev(GridParallelRNG &RNG, LinearOperatorBase< FineField > &hermop, int nn, double hi, double lo, int orderfilter, int ordermin, int orderstep, double filterlo)
Definition Aggregates.h:181
Aggregation::CreateSubspacePolyCheby
virtual void CreateSubspacePolyCheby(GridParallelRNG &RNG, LinearOperatorBase< FineField > &hermop, int nn, double hi, double lo1, int orderfilter, double lo2, int orderstep)
Definition Aggregates.h:299
Aggregation::RefineSubspace
virtual void RefineSubspace(LinearOperatorBase< FineField > &hermop, double Lo, double tol, int maxit)
Definition Aggregates.h:555
Aggregation::CreateSubspace
virtual void CreateSubspace(GridParallelRNG &RNG, LinearOperatorBase< FineField > &hermop, int nn=nbasis)
Definition Aggregates.h:95
Aggregation::checkerboard
int checkerboard
Definition Aggregates.h:59
Aggregation::CreateSubspaceChebyshevPowerLaw
virtual void CreateSubspaceChebyshevPowerLaw(GridParallelRNG &RNG, LinearOperatorBase< FineField > &hermop, int nn, double hi, int orderfilter)
Definition Aggregates.h:408
Aggregation::CreateSubspaceRandom
virtual void CreateSubspaceRandom(GridParallelRNG &RNG)
Definition Aggregates.h:83
Aggregation::PromoteFromSubspace
void PromoteFromSubspace(const CoarseVector &CoarseVec, FineField &FineVec)
Definition Aggregates.h:78
Aggregation::siteVector
iVector< CComplex, nbasis > siteVector
Definition Aggregates.h:49
Aggregation::CoarseScalar
Lattice< CComplex > CoarseScalar
Definition Aggregates.h:53
Aggregation::Checkerboard
int Checkerboard(void)
Definition Aggregates.h:60
Aggregation::CreateSubspaceGCR
virtual void CreateSubspaceGCR(GridParallelRNG &RNG, LinearOperatorBase< FineField > &DiracOp, int nn=nbasis)
Definition Aggregates.h:129
Aggregation::Aggregation
Aggregation(GridBase *_CoarseGrid, GridBase *_FineGrid, int _checkerboard)
Definition Aggregates.h:61
Aggregation::RefineSubspaceHDCG
virtual void RefineSubspaceHDCG(LinearOperatorBase< FineField > &hermop, TwoLevelADEF2mrhs< FineField, CoarseVector > &theHDCG, int nrhs)
Definition Aggregates.h:571
Aggregation::Nbasis
constexpr int Nbasis(void)
Definition Aggregates.h:47
Aggregation::CoarseGrid
GridBase * CoarseGrid
Definition Aggregates.h:56
Aggregation::ProjectToSubspace
void ProjectToSubspace(CoarseVector &CoarseVec, const FineField &FineVec)
Definition Aggregates.h:75
Aggregation::CreateSubspaceChebyshevNew
virtual void CreateSubspaceChebyshevNew(GridParallelRNG &RNG, LinearOperatorBase< FineField > &hermop, double hi)
Definition Aggregates.h:444
Aggregation::CoarseMatrix
Lattice< iMatrix< CComplex, nbasis > > CoarseMatrix
Definition Aggregates.h:51
Aggregation::CreateSubspaceChebyshev
virtual void CreateSubspaceChebyshev(GridParallelRNG &RNG, LinearOperatorBase< FineField > &hermop, int nn, double hi, double lo, int orderfilter)
Definition Aggregates.h:361
Aggregation::FineField
Lattice< Fobj > FineField
Definition Aggregates.h:54
Aggregation::CreateSubspaceMultishift
virtual void CreateSubspaceMultishift(GridParallelRNG &RNG, LinearOperatorBase< FineField > &hermop, double Lo, double tol, int maxit)
Definition Aggregates.h:504
LatticeAccelerator::Checkerboard
accelerator_inline int Checkerboard(void) const
Definition Lattice_view.h:37
LinearOperatorBase::Op
virtual void Op(const Field &in, Field &out)=0
LinearOperatorBase::AdjOp
virtual void AdjOp(const Field &in, Field &out)=0
LinearOperatorBase::HermOp
virtual void HermOp(const Field &in, Field &out)=0
MultiShiftFunction::poles
std::vector< RealD > poles
Definition MultiShiftFunction.h:37
MultiShiftFunction::tolerances
std::vector< RealD > tolerances
Definition MultiShiftFunction.h:39
MultiShiftFunction::residues
std::vector< RealD > residues
Definition MultiShiftFunction.h:38
Zero
Definition Simd.h:194