1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file sys_HMatrix.cpp
19 \author E. Lunéville
20 \since 17 may 2016
21 \date 17 may 2016
22
23 Low level tests of Hmatrix class.
24 Almost functionalities are checked.
25 This function may either creates a reference file storing the results (check=false)
26 or compares results to those stored in the reference file (check=true)
27 It returns reporting informations in a string
28 */
29
30 #include "xlife++-libs.h"
31 #include "hierarchicalMatrix.h"
32 #include "testUtils.hpp"
33
34 #include <iostream>
35 #include <fstream>
36 #include <vector>
37
38 using namespace xlifepp;
39
40 namespace sys_HMatrix
41 {
42
sys_HMatrix(bool check)43 String sys_HMatrix(bool check)
44 {
45 String rootname = "sys_HMatrix";
46 trace_p->push(rootname);
47 std::stringstream out; // string stream receiving results
48 out.precision(testPrec);
49 verboseLevel(20);
50 //Trace::trackingMode=true;
51 numberOfThreads(4);
52
53 Real eps=1.E-08;
54 Real eps2=1.E-04;
55 Number nr=5;
56 Number order=5;
57 IntegrationMethods ims(_SauterSchwabIM,order,0.,_defaultRule,4,2.,_defaultRule,3);
58 bool laplace_sl=true, laplace_dl=true, helmholtz_sl=false;
59 bool hm_exact = false, hm_exact_no_sym=false,
60 hm_svd=false, hm_svd_rank=false, hm_svd_eps=false,
61 hm_aca = false, hm_aca_part=false, hm_aca_plus=true,
62 hm_rsvd_rank = false, hm_rsvd_eps=false, hm_r3svd=false;
63 MeshGenerator mg=_gmsh; //_subdiv
64
65 for(number_t k=5; k<30; k+=5)
66 {
67 //test on a sphere
68 Mesh meshd(Sphere(_center=Point(0.,0.,0.),_radius=1.,_nnodes=k,_domain_name="Omega"),_triangle,1,mg);
69 Domain omega=meshd.domain("Omega");
70 Space W(omega,P0,"V",false);
71 Unknown u(W,"u"); TestFunction v(u,"v");
72 Vector<Real> x(W.dimSpace(),1.);
73 Number sb=std::max(number_t(10),W.dimSpace()/100);
74 Number nbp=50;
75 Kernel Gl=Laplace3dKernel();
76 HMatrixIM him0(_cardinalityBisection, sb,sb,_noHMApproximation,0.,ims);
77 HMatrixIM him(_cardinalityBisection, sb,sb,_svdCompression,0.,ims);
78 HMatrixIM him1(_cardinalityBisection, sb,sb,_svdCompression,eps,ims);
79 HMatrixIM him2(_cardinalityBisection, sb,sb,_svdCompression,nr,ims);
80 HMatrixIM him3(_cardinalityBisection, sb,sb,_acaFull,eps2, ims);
81 HMatrixIM him4(_cardinalityBisection, sb,sb,_acaPartial,eps2,ims);
82 HMatrixIM him5(_cardinalityBisection, sb,sb,_acaPlus,eps2,ims);
83 HMatrixIM him6(_cardinalityBisection, sb,sb,_rsvdCompression,nr,ims);
84 HMatrixIM him8(_cardinalityBisection, sb,sb,_rsvdCompression,eps2,ims);
85 HMatrixIM him7(_cardinalityBisection, sb,sb,_r3svdCompression,eps,ims);
86 thePrintStream<<"============================================================="<<eol;
87 thePrintStream<<" sphere mesh k = "<<k<<", nbdofs = "<<W.dimSpace()<<", block size = "<<sb<<eol;
88 thePrintStream<<"============================================================="<<eol;
89 std::cout<<"---------- sphere mesh k = "<<k<<", nbdofs = "<<W.dimSpace()<<", block size = "<<sb<<" ----------"<<eol;
90 std::cout<<"-----------------------------------------------------------"<<eol;
91
92 if(laplace_sl)
93 {
94 thePrintStream<<"============================================================="<<eol;
95 thePrintStream<<" Test for Laplace kernel, single layer"<<eol;
96 thePrintStream<<"============================================================="<<eol;
97 std::cout<<"------------ Test for Laplace kernel, single layer ----------"<<eol;
98 std::cout<<"----------------- Dense matrix computation -------------------"<<eol;
99
100 BilinearForm blf=intg(omega,omega,u*Gl*v,ims);
101 elapsedTime();
102 TermMatrix A(blf, "A");
103 elapsedTime("build full matrix");
104 LargeMatrix<Real>& Alm=A.getLargeMatrix<Real>();
105 Vector<Real> Ax;
106 elapsedTime();
107 for(Number k=0; k<nbp; ++k) Ax=Alm*x;
108 thePrintStream<<"|A*x| = "<<norm(Ax)<<eol;
109 elapsedTime("full matrix x vector");
110 clear(A);
111
112 if(hm_exact) //no approximation sym
113 {
114 std::cout<<"------------------- Hmatrix no approximation -------------------"<<eol;
115 BilinearForm blfhm0=intg(omega,omega,u*Gl*v,him0);
116 elapsedTime();
117 TermMatrix Ahm0(blfhm0,"Ahm0");
118 elapsedTime("build HMatrix exact ");
119 HMatrix<Real,FeDof>& hmA0=Ahm0.getHMatrix<Real>();
120 Vector<Real> Ahm0x;
121 elapsedTime();
122 for(Number k=0; k<nbp; ++k) Ahm0x=hmA0*x;
123 elapsedTime("exact Hmatrix x vector");
124 thePrintStream<<"exact |A*x- hmA0*x| = "<<norm(Ax-Ahm0x)<<eol;
125 }
126
127 if(hm_exact_no_sym) // no approximation no symmetry
128 {
129 std::cout<<"------------------- Hmatrix no approximation -------------------"<<eol;
130 BilinearForm blfhm00=intg(omega,omega,u*Gl*v,him0,_noSymmetry);
131 elapsedTime();
132 TermMatrix Ahm00(blfhm00,"Ahm00");
133 elapsedTime("build HMatrix exact non sym");
134 HMatrix<Real,FeDof>& hmA00=Ahm00.getHMatrix<Real>();
135 Vector<Real> Ahm00x;
136 elapsedTime();
137 for(Number k=0; k<nbp; ++k) Ahm00x=hmA00*x;
138 elapsedTime("exact Hmatrix non sym x vector");
139 thePrintStream<<"exact non sym |A*x- hmA00*x| = "<<norm(Ax-Ahm00x)<<eol;
140 }
141
142 if(hm_svd)//full svd
143 {
144 std::cout<<"------------------- Hmatrix full svd -------------------"<<eol;
145 BilinearForm blfhm=intg(omega,omega,u*Gl*v,him);
146 elapsedTime();
147 TermMatrix Ahm(blfhm,"Ahm");
148 elapsedTime("build HMatrix full svd");
149 HMatrix<real_t,FeDof>& hmA=Ahm.getHMatrix<Real>();
150 Vector<Real> Ahmx;
151 elapsedTime();
152 for(Number k=0; k<nbp; ++k) Ahmx=hmA*x;
153 elapsedTime("full svd Hmatrix x vector");
154 thePrintStream<<"full svd |A*x-hmA*x| = "<<norm(Ax-Ahmx)<<eol;
155 }
156
157 if(hm_svd_eps) //svd with truncation of singular values sigma_i> eps
158 {
159 std::cout<<"------------------- Hmatrix truncated svd, sigma_n > "<<eps<<eol;
160 BilinearForm blfhm1=intg(omega,omega,u*Gl*v,him1);
161 elapsedTime();
162 TermMatrix Ahm1(blfhm1,"Ahm1");
163 elapsedTime("build HMatrix eps truncated svd");
164 HMatrix<real_t,FeDof>& hmA1=Ahm1.getHMatrix<Real>();
165 Vector<Real> Ahm1x;
166 elapsedTime();
167 for(Number k=0; k<nbp; ++k) Ahm1x=hmA1*x;
168 elapsedTime("eps truncated svd Hmatrix x vector");
169 thePrintStream<<"truncated svd (eps="<<eps<<") |A*x-hmA1*x| = "<<norm(Ax-Ahm1x)<<" average rank = "<<hmA1.averageRank()<<eol;
170 }
171
172 if(hm_svd_rank) //svd with truncation of n first singular values
173 {
174 std::cout<<"------------------- Hmatrix truncated svd, n < "<<nr<<eol;
175 BilinearForm blfhm2=intg(omega,omega,u*Gl*v,him2);
176 elapsedTime();
177 TermMatrix Ahm2(blfhm2,"Ahm2");
178 elapsedTime("build HMatrix n truncated svd");
179 HMatrix<real_t,FeDof>& hmA2=Ahm2.getHMatrix<Real>();
180 Vector<Real> Ahm2x;
181 elapsedTime();
182 for(Number k=0; k<nbp; ++k) Ahm2x=hmA2*x;
183 elapsedTime("n truncated svd Hmatrix x vector");
184 thePrintStream<<"truncated svd (n="<<nr<<") |A*x-hmA2*x| = "<<norm(Ax-Ahm2x)<<" average rank = "<<hmA2.averageRank()<<eol;
185 }
186
187 if(hm_aca) //compression ACA full
188 {
189 std::cout<<"------------------- Hmatrix ACA full -------------------"<<eol;
190 BilinearForm blfhm3=intg(omega,omega,u*Gl*v,him3);
191 elapsedTime();
192 TermMatrix Ahm3(blfhm3,"Ahm3");
193 elapsedTime("build HMatrix ACA full");
194 HMatrix<real_t,FeDof>& hmA3=Ahm3.getHMatrix<Real>();
195 Vector<Real> Ahm3x;
196 elapsedTime();
197 for(Number k=0; k<nbp; ++k) Ahm3x=hmA3*x;
198 elapsedTime("ACA full Hmatrix x vector");
199 thePrintStream<<"ACA full (eps="<<eps<<") |A*x-hmA3*x| = "<<norm(Ax-Ahm3x)<<" average rank = "<<hmA3.averageRank()<<eol;
200 }
201
202 if(hm_aca_part)//compression ACA partial
203 {
204 std::cout<<"------------------- Hmatrix ACA partial -------------------"<<eol;
205 BilinearForm blfhm4=intg(omega,omega,u*Gl*v,him4);
206 elapsedTime();
207 TermMatrix Ahm4(blfhm4,"Ahm4");
208 elapsedTime("build HMatrix ACA partial");
209 HMatrix<real_t,FeDof>& hmA4=Ahm4.getHMatrix<Real>();
210 Vector<Real> Ahm4x;
211 elapsedTime();
212 for(Number k=0; k<nbp; ++k) Ahm4x=hmA4*x;
213 elapsedTime("ACA partial Hmatrix x vector");
214 thePrintStream<<"ACA partial (eps="<<eps<<") |A*x-hmA4*x| = "<<norm(Ax-Ahm4x)<<" average rank = "<<hmA4.averageRank()<<eol;
215 }
216
217 if(hm_aca_plus)//compression ACA plus
218 {
219 std::cout<<"------------------- Hmatrix ACA + -------------------"<<eol;
220 BilinearForm blfhm5=intg(omega,omega,u*Gl*v,him5);
221 elapsedTime();
222 TermMatrix Ahm5(blfhm5,"Ahm5");
223 elapsedTime("build HMatrix ACA+");
224 HMatrix<real_t,FeDof>& hmA5=Ahm5.getHMatrix<Real>();
225 Vector<Real> Ahm5x;
226 elapsedTime();
227 for(Number k=0; k<nbp; ++k) Ahm5x=hmA5*x;
228 elapsedTime("ACA+ Hmatrix x vector");
229 thePrintStream<<"ACA+ (eps="<<eps<<") |A*x-hmA5*x| = "<<norm(Ax-Ahm5x)<<" average rank = "<<hmA5.averageRank()<<eol;
230 }
231
232 if(hm_rsvd_rank) //rsvd compression rank
233 {
234 std::cout<<"------------------- Hmatrix rsvd rank -------------------"<<eol;
235 BilinearForm blfhm6=intg(omega,omega,u*Gl*v,him6);
236 elapsedTime();
237 TermMatrix Ahm6(blfhm6,"Ahm6");
238 elapsedTime("build HMatrix rsvd_rank");
239 HMatrix<real_t,FeDof>& hmA6=Ahm6.getHMatrix<Real>();
240 Vector<Real> Ahm6x;
241 elapsedTime();
242 for(Number k=0; k<nbp; ++k) Ahm6x=hmA6*x;
243 elapsedTime("rsvd_rank Hmatrix x vector");
244 thePrintStream<<"rsvd (nr="<<nr<<") |A*x-hmA6*x| = "<<norm(Ax-Ahm6x)<<" average rank = "<<hmA6.averageRank()<<eol;
245 }
246
247
248 if(hm_rsvd_eps) //rsvd compression eps
249 {
250 std::cout<<"------------------- Hmatrix rsvd eps -------------------"<<eol;
251 BilinearForm blfhm8=intg(omega,omega,u*Gl*v,him8);
252 elapsedTime();
253 TermMatrix Ahm8(blfhm8,"Ahm8");
254 elapsedTime("build HMatrix rsvd_eps");
255 HMatrix<real_t,FeDof>& hmA8=Ahm8.getHMatrix<Real>();
256 Vector<Real> Ahm8x;
257 elapsedTime();
258 for(Number k=0; k<nbp; ++k) Ahm8x=hmA8*x;
259 elapsedTime("rsvd_eps Hmatrix x vector");
260 thePrintStream<<"rsvd (eps="<<eps2<<") |A*x-hmA8*x| = "<<norm(Ax-Ahm8x)<<" average rank = "<<hmA8.averageRank()<<eol;
261 }
262
263 if(hm_r3svd) //r3svd compression
264 {
265 std::cout<<"------------------- Hmatrix r3svd -------------------"<<eol;
266 BilinearForm blfhm7=intg(omega,omega,u*Gl*v,him7);
267 elapsedTime();
268 TermMatrix Ahm7(blfhm7,"Ahm7");
269 elapsedTime("build HMatrix r3svd");
270 HMatrix<real_t,FeDof>& hmA7=Ahm7.getHMatrix<Real>();
271 Vector<Real> Ahm7x;
272 elapsedTime();
273 for(Number k=0; k<nbp; ++k) Ahm7x=hmA7*x;
274 elapsedTime("r3svd Hmatrix x vector");
275 thePrintStream<<"r3svd (eps="<<eps<<") |A*x-hmA7*x| = "<<norm(Ax-Ahm7x)<<" average rank = "<<hmA7.averageRank()<<eol;
276 }
277 }
278
279 //test with a combination of a HMatrix and a LargeMatrix
280 //======================================================
281 if(laplace_dl)
282 {
283 thePrintStream<<"============================================================="<<eol;
284 thePrintStream<<" Test for Laplace kernel, double layer, hybrid HM/LM"<<eol;
285 thePrintStream<<"============================================================="<<eol;
286 std::cout<<"----- Test for Laplace kernel, double layer, hybrid HM/LM ------"<<eol;
287 std::cout<<"------------------ Dense matrix computation --------------------"<<eol;
288 BilinearForm blm=intg(omega,omega,u*ndotgrad_y(Gl)*v,ims)-0.5*intg(omega,u*v);
289 elapsedTime();
290 TermMatrix B(blm,"B");
291 elapsedTime("build full matrix");
292 Vector<Real> Bx;
293 LargeMatrix<Real>& Blm=B.getLargeMatrix<Real>();
294 elapsedTime();
295 for(Number k=0; k<nbp; ++k) Bx=Blm*x;
296 thePrintStream<<"|B*x| = "<<norm(Bx)<<eol;
297 elapsedTime("full matrix x vector");
298 clear(B);
299
300 if(hm_exact) //no approximation
301 {
302 std::cout<<"------------------- Hmatrix no approximation -------------------"<<eol;
303 BilinearForm blm0=intg(omega,omega,u*ndotgrad_y(Gl)*v,him0)-0.5*intg(omega,u*v);
304 TermMatrix B0(blm0, "B0");
305 elapsedTime("build exact HM/LM matrix");
306 Vector<Real> B0x;
307 HMatrix<real_t,FeDof>& B0hm=B0.getHMatrix<Real>();
308 elapsedTime();
309 for(Number k=0; k<nbp; ++k) B0x=B0hm*x;
310 thePrintStream<<"exact HM/LM |B*x-B0*x| = "<<norm(Bx-B0x)<<eol;
311 elapsedTime("exact HM/LM matrix x vector");
312 }
313
314 if(hm_svd_eps) //svd with truncation of singular values sigma_i> eps
315 {
316 std::cout<<"------------------- Hmatrix truncated svd, sigma_n > "<<eps<<eol;
317 BilinearForm blm1=intg(omega,omega,u*ndotgrad_y(Gl)*v,him1)-0.5*intg(omega,u*v);
318 TermMatrix B1(blm1, "B1");
319 elapsedTime("build SVD eps HM/LM matrix");
320 Vector<Real> B1x;
321 HMatrix<real_t,FeDof>& B1hm=B1.getHMatrix<Real>();
322 elapsedTime();
323 for(Number k=0; k<nbp; ++k) B1x=B1hm*x;
324 thePrintStream<<"truncated svd (eps="<<eps<<") |B*x-B1*x| = "<<norm(Bx-B1x)<<" average rank = "<<B1hm.averageRank()<<eol;
325 elapsedTime("SVD eps HM/LM matrix x vector");
326 }
327
328 if(hm_svd_rank) //svd with truncation of n first singular values
329 {
330 std::cout<<"------------------- Hmatrix truncated svd, n < "<<nr<<eol;
331 BilinearForm blm2=intg(omega,omega,u*ndotgrad_y(Gl)*v,him2)-0.5*intg(omega,u*v);
332 TermMatrix B2(blm2, "B2");
333 elapsedTime("build SVD n HM/LM matrix");
334 Vector<Real> B2x;
335 HMatrix<real_t,FeDof>& B2hm=B2.getHMatrix<Real>();
336 elapsedTime();
337 for(Number k=0; k<nbp; ++k) B2x=B2hm*x;
338 thePrintStream<<"truncated svd (n="<<nr<<") |B*x-B2*x| = "<<norm(Bx-B2x)<<" average rank = "<<B2hm.averageRank()<<eol;
339 elapsedTime("SVD eps HM/LM matrix x vector");
340 }
341
342 if(hm_aca) //ACA compression
343 {
344 std::cout<<"------------------- Hmatrix ACA full -------------------"<<eol;
345 BilinearForm blm3=intg(omega,omega,u*ndotgrad_y(Gl)*v,him3)-0.5*intg(omega,u*v);
346 TermMatrix B3(blm3, "B3");
347 elapsedTime("build ACA matrix");
348 Vector<Real> B3x;
349 HMatrix<real_t,FeDof>& B3hm=B3.getHMatrix<Real>();
350 elapsedTime();
351 for(Number k=0; k<nbp; ++k) B3x=B3hm*x;
352 thePrintStream<<"ACA full (eps="<<eps<<") |B*x-B3*x| = "<<norm(Bx-B3x)<<" average rank = "<<B3hm.averageRank()<<eol;
353 elapsedTime("ACA full HM/LM matrix x vector");
354 }
355
356
357 if(hm_aca_part) //partial ACA compression
358 {
359 std::cout<<"------------------- Hmatrix ACA partial -------------------"<<eol;
360 BilinearForm blm4=intg(omega,omega,u*ndotgrad_y(Gl)*v,him4)-0.5*intg(omega,u*v);
361 TermMatrix B4(blm4, "B4");
362 elapsedTime("build ACA Part matrix");
363 Vector<Real> B4x;
364 HMatrix<real_t,FeDof>& B4hm=B4.getHMatrix<Real>();
365 elapsedTime();
366 for(Number k=0; k<nbp; ++k) B4x=B4hm*x;
367 thePrintStream<<"ACA partial (eps="<<eps<<") |B*x-B4*x| = "<<norm(Bx-B4x)<<" average rank = "<<B4hm.averageRank()<<eol;
368 elapsedTime("ACA Part HM/LM matrix x vector");
369 }
370
371 if(hm_aca_plus)//ACA+ compression
372 {
373 std::cout<<"------------------- Hmatrix ACA + -------------------"<<eol;
374 BilinearForm blm5=intg(omega,omega,u*ndotgrad_y(Gl)*v,him5)-0.5*intg(omega,u*v);
375 TermMatrix B5(blm5, "B5");
376 elapsedTime("build ACA Part matrix");
377 Vector<Real> B5x;
378 HMatrix<real_t,FeDof>& B5hm=B5.getHMatrix<Real>();
379 elapsedTime();
380 for(Number k=0; k<nbp; ++k) B5x=B5hm*x;
381 thePrintStream<<"ACA plus (eps="<<eps<<")|B*x-B5*x| = "<<norm(Bx-B5x)<<" average rank = "<<B5hm.averageRank()<<eol;
382 elapsedTime("ACA+ HM/LM matrix x vector");
383 }
384
385 if(hm_rsvd_rank)//rsvd compression
386 {
387 std::cout<<"------------------- Hmatrix rsvd_rank -------------------"<<eol;
388 BilinearForm blm6=intg(omega,omega,u*ndotgrad_y(Gl)*v,him6)-0.5*intg(omega,u*v);
389 TermMatrix B6(blm6, "B6");
390 elapsedTime("build RSVD_rank matrix");
391 Vector<Real> B6x;
392 HMatrix<real_t,FeDof>& B6hm=B6.getHMatrix<Real>();
393 elapsedTime();
394 for(Number k=0; k<nbp; ++k) B6x=B6hm*x;
395 thePrintStream<<"rsvd_rank (nr="<<nr<<") |B*x-B6*x| = "<<norm(Bx-B6x)<<" average rank = "<<B6hm.averageRank()<<eol;
396 elapsedTime("rsvd_rank matrix x vector");
397 }
398
399 if(hm_rsvd_eps)//rsvd_eps compression
400 {
401 std::cout<<"------------------- Hmatrix rsvd_eps-------------------"<<eol;
402 BilinearForm blm8=intg(omega,omega,u*ndotgrad_y(Gl)*v,him8)-0.5*intg(omega,u*v);
403 TermMatrix B8(blm8, "B8");
404 elapsedTime("build R3SVD matrix");
405 Vector<Real> B8x;
406 HMatrix<real_t,FeDof>& B8hm=B8.getHMatrix<Real>();
407 elapsedTime();
408 for(Number k=0; k<nbp; ++k) B8x=B8hm*x;
409 thePrintStream<<"rsvd (eps="<<eps2<<") |B*x-B8*x| = "<<norm(Bx-B8x)<<" average rank = "<<B8hm.averageRank()<<eol;
410 elapsedTime("rsvd_eps matrix x vector");
411 }
412
413 if(hm_r3svd) //r3svd compression
414 {
415 std::cout<<"------------------- Hmatrix r3svd -------------------"<<eol;
416 BilinearForm blm7=intg(omega,omega,u*ndotgrad_y(Gl)*v,him7)-0.5*intg(omega,u*v);
417 TermMatrix B7(blm7, "B7");
418 elapsedTime("build R3SVD matrix");
419 Vector<Real> B7x;
420 HMatrix<real_t,FeDof>& B7hm=B7.getHMatrix<Real>();
421 elapsedTime();
422 for(Number k=0; k<nbp; ++k) B7x=B7hm*x;
423 thePrintStream<<"r3svd (eps="<<eps<<") |B*x-B7*x| = "<<norm(Bx-B7x)<<" average rank = "<<B7hm.averageRank()<<eol;
424 elapsedTime("R3SVD HM/LM matrix x vector");
425 }
426 }
427
428 // test for Helmholtz Kernel
429 // ============================================================================
430 if(helmholtz_sl)
431 {
432 thePrintStream<<"============================================================="<<eol;
433 thePrintStream<<" Test for Helmhotz kernel, single layer"<<eol;
434 thePrintStream<<"============================================================="<<eol;
435 std::cout<<"---------- Test for Helmhotz kernel, single layer ------------"<<eol;
436 std::cout<<"----------------- Dense matrix computation -------------------"<<eol;
437 Kernel Gh=Helmholtz3dKernel(2.);
438 Vector<Complex> z(W.dimSpace(),Complex(1.));
439 BilinearForm cblm=intg(omega,omega,u*Gh*v,ims);
440 elapsedTime();
441 TermMatrix C(cblm,"C");
442 elapsedTime("build full matrix");
443 Vector<Complex> Cx;
444 LargeMatrix<Complex>& Clm=C.getLargeMatrix<Complex>();
445 elapsedTime();
446 for(Number k=0; k<nbp; ++k) Cx=Clm*z;
447 thePrintStream<<"|C*x| = "<<norm(Cx)<<eol;
448 elapsedTime("full matrix x vector");
449 clear(C);
450
451 if(hm_exact) //no approximation sym
452 {
453 std::cout<<"------------------- Hmatrix no approximation -------------------"<<eol;
454 BilinearForm cblfhm0=intg(omega,omega,u*Gh*v,him0);
455 elapsedTime();
456 TermMatrix Chm0(cblfhm0,"Chm0");
457 elapsedTime("build HMatrix exact");
458 HMatrix<Complex,FeDof>& hmC0=Chm0.getHMatrix<Complex>();
459 Vector<Complex> Chm0x;
460 elapsedTime();
461 for(Number k=0; k<nbp; ++k) Chm0x=hmC0*z;
462 elapsedTime("exact Hmatrix x vector");
463 thePrintStream<<"exact |C*x- hmC0*x| = "<<norm(Cx-Chm0x)<<eol;
464 }
465
466 if(hm_exact_no_sym) //no approximation no sym
467 {
468 std::cout<<"------------------- Hmatrix no approximation no symmetry -------------------"<<eol;
469 BilinearForm cblfhm00=intg(omega,omega,u*Gh*v,him0,_noSymmetry);
470 elapsedTime();
471 TermMatrix Chm00(cblfhm00,"Chm00");
472 elapsedTime("build HMatrix exact non sym");
473 HMatrix<Complex,FeDof>& hmC00=Chm00.getHMatrix<Complex>();
474 Vector<Complex> Chm00x;
475 elapsedTime();
476 for(Number k=00; k<nbp; ++k) Chm00x=hmC00*z;
477 elapsedTime("exact Hmatrix non sym x vector");
478 thePrintStream<<"exact non sym |C*x- hmC00*x| = "<<norm(Cx-Chm00x)<<eol;
479 }
480
481 if(hm_svd) //full svd
482 {
483 std::cout<<"------------------- Hmatrix full svd -------------------"<<eol;
484 BilinearForm cblfhm=intg(omega,omega,u*Gh*v,him);
485 elapsedTime();
486 TermMatrix Chm(cblfhm,"Chm");
487 elapsedTime("build HMatrix full svd");
488 HMatrix<Complex,FeDof>& hmC=Chm.getHMatrix<Complex>();
489 Vector<Complex> Chmx;
490 elapsedTime();
491 for(Number k=0; k<nbp; ++k) Chmx=hmC*z;
492 elapsedTime("full svd Hmatrix x vector");
493 thePrintStream<<"full svd |C*x-hmC*x| = "<<norm(Cx-Chmx)<<eol;
494 }
495
496 if(hm_svd_eps) //svd with truncation of singular values sigma_i> eps
497 {
498 std::cout<<"------------------- Hmatrix trucated svd, sigma_n > "<<eps<<" -------------------"<<eol;
499 BilinearForm cblfhm1=intg(omega,omega,u*Gh*v,him1);
500 elapsedTime();
501 TermMatrix Chm1(cblfhm1,"Chm1");
502 elapsedTime("build HMatrix eps truncated svd");
503 HMatrix<Complex,FeDof>& hmC1=Chm1.getHMatrix<Complex>();
504 Vector<Complex> Chm1x;
505 elapsedTime();
506 for(Number k=0; k<nbp; ++k) Chm1x=hmC1*z;
507 elapsedTime("eps truncated svd Hmatrix x vector");
508 thePrintStream<<"truncated svd (eps="<<eps<<") |C*x-hmC1*x| = "<<norm(Cx-Chm1x)<<" average rank = "<<hmC1.averageRank()<<eol;
509 }
510
511
512 if(hm_svd_rank) // svd with truncation of singular values (n<=nr)
513 {
514 std::cout<<"------------------- Hmatrix truncated svd, n < "<<nr<<" -------------------"<<eol;
515 BilinearForm cblfhm2=intg(omega,omega,u*Gh*v,him2);
516 elapsedTime();
517 TermMatrix Chm2(cblfhm2,"Chm2");
518 elapsedTime("build HMatrix n truncated svd");
519 HMatrix<Complex,FeDof>& hmC2=Chm2.getHMatrix<Complex>();
520 Vector<Complex> Chm2x;
521 elapsedTime();
522 for(Number k=0; k<nbp; ++k) Chm2x=hmC2*z;
523 elapsedTime("n truncated svd Hmatrix x vector");
524 thePrintStream<<"truncated svd (n="<<nr<<") |C*x-hmC2*x| = "<<norm(Cx-Chm2x)<<" average rank = "<<hmC2.averageRank()<<eol;
525 }
526
527 if(hm_aca) //compression ACA full
528 {
529 std::cout<<"------------------- Hmatrix ACA full -------------------"<<eol;
530 BilinearForm cblfhm3=intg(omega,omega,u*Gh*v,him3);
531 elapsedTime();
532 TermMatrix Chm3(cblfhm3,"Chm3");
533 elapsedTime("build HMatrix ACA full");
534 HMatrix<Complex,FeDof>& hmC3=Chm3.getHMatrix<Complex>();
535 Vector<Complex> Chm3x;
536 elapsedTime();
537 for(Number k=0; k<nbp; ++k) Chm3x=hmC3*z;
538 elapsedTime("ACA full Hmatrix x vector");
539 thePrintStream<<"ACA full (eps="<<eps<<") |C*x-hmC3*x| = "<<norm(Cx-Chm3x)<<" average rank = "<<hmC3.averageRank()<<eol;
540 }
541
542 if(hm_aca_part)//compression ACA partial
543 {
544 std::cout<<"------------------- Hmatrix ACA partial -------------------"<<eol;
545 BilinearForm cblfhm4=intg(omega,omega,u*Gh*v,him4);
546 elapsedTime();
547 TermMatrix Chm4(cblfhm4,"Chm4");
548 elapsedTime("build HMatrix ACA partial");
549 HMatrix<Complex,FeDof>& hmC4=Chm4.getHMatrix<Complex>();
550 Vector<Complex> Chm4x;
551 elapsedTime();
552 for(Number k=0; k<nbp; ++k) Chm4x=hmC4*z;
553 elapsedTime("ACA partial Hmatrix x vector");
554 thePrintStream<<"ACA partial (eps="<<eps<<") |C*x-hmC4*x| = "<<norm(Cx-Chm4x)<<" average rank = "<<hmC4.averageRank()<<eol;
555 }
556
557 if(hm_aca_plus) //compression ACA plus
558 {
559 std::cout<<"------------------- Hmatrix ACA + -------------------"<<eol;
560 BilinearForm cblfhm5=intg(omega,omega,u*Gh*v,him5);
561 elapsedTime();
562 TermMatrix Chm5(cblfhm5,"Chm5");
563 elapsedTime("build HMatrix ACA plus");
564 HMatrix<Complex,FeDof>& hmC5=Chm5.getHMatrix<Complex>();
565 Vector<Complex> Chm5x;
566 elapsedTime();
567 for(Number k=0; k<nbp; ++k) Chm5x=hmC5*z;
568 elapsedTime("ACA plus Hmatrix x vector");
569 thePrintStream<<"ACA plus (eps="<<eps<<") |C*x-hmC5*x| = "<<norm(Cx-Chm5x)<<" average rank = "<<hmC5.averageRank()<<eol;
570 }
571
572 if(hm_rsvd_rank)//rsvd rank compression
573 {
574 std::cout<<"------------------- Hmatrix rsvd -------------------"<<eol;
575 BilinearForm cblfhm6=intg(omega,omega,u*Gh*v,him6);
576 elapsedTime();
577 TermMatrix Chm6(cblfhm6,"Chm6");
578 elapsedTime("build HMatrix rsvd_rank");
579 HMatrix<Complex,FeDof>& hmC6=Chm6.getHMatrix<Complex>();
580 Vector<Complex> Chm6x;
581 elapsedTime();
582 for(Number k=0; k<nbp; ++k) Chm6x=hmC6*z;
583 elapsedTime("rsvd_rank Hmatrix x vector");
584 thePrintStream<<"rsvd (nr="<<nr<<") |C*x-hmC6*x| = "<<norm(Cx-Chm6x)<<" average rank = "<<hmC6.averageRank()<<eol;
585 }
586
587 if(hm_rsvd_eps)//rsvd eps compression
588 {
589 std::cout<<"------------------- Hmatrix rsvd eps -------------------"<<eol;
590 BilinearForm cblfhm8=intg(omega,omega,u*Gh*v,him8);
591 elapsedTime();
592 TermMatrix Chm8(cblfhm8,"Chm8");
593 elapsedTime("build HMatrix rsvd_eps");
594 HMatrix<Complex,FeDof>& hmC8=Chm8.getHMatrix<Complex>();
595 Vector<Complex> Chm8x;
596 elapsedTime();
597 for(Number k=0; k<nbp; ++k) Chm8x=hmC8*z;
598 elapsedTime("rsvd_eps Hmatrix x vector");
599 thePrintStream<<"rsvd (eps="<<eps2<<") |C*x-hmC8*x| = "<<norm(Cx-Chm8x)<<" average rank = "<<hmC8.averageRank()<<eol;
600 }
601
602 if(hm_r3svd) //r3svd compression
603 {
604 std::cout<<"------------------- Hmatrix r3svd -------------------"<<eol;
605 BilinearForm cblfhm7=intg(omega,omega,u*Gh*v,him7);
606 elapsedTime();
607 TermMatrix Chm7(cblfhm7,"Chm7");
608 elapsedTime("build HMatrix r3svd");
609 HMatrix<Complex,FeDof>& hmC7=Chm7.getHMatrix<Complex>();
610 Vector<Complex> Chm7x;
611 elapsedTime();
612 for(Number k=0; k<nbp; ++k) Chm7x=hmC7*z;
613 elapsedTime("r3svd Hmatrix x vector");
614 thePrintStream<<"r3svd (eps="<<eps<<") |C*x-hmC7*x| = "<<norm(Cx-Chm7x)<<" average rank = "<<hmC7.averageRank()<<eol;
615 }
616 }
617 }
618
619 //------------------------------------------------------------------------------------
620 // save results in a file or compare results with some references value in a file
621 //------------------------------------------------------------------------------------
622 trace_p->pop();
623 if(check) { return diffResFile(out, rootname); }
624 else { return saveResToFile(out, rootname); }
625 }
626
627 }
628