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 unit_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 #undef __STRICT_ANSI__ // _controlfp is a non-standard function documented in MSDN
39 #include <float.h>
40 #include <stdio.h>
41
42 using namespace xlifepp;
43
44 namespace unit_HMatrix
45 {
46
unit_HMatrix(bool check)47 String unit_HMatrix(bool check)
48 {
49 String rootname = "unit_HMatrix";
50 trace_p->push(rootname);
51 std::stringstream out; // string stream receiving results
52 out.precision(testPrec);
53 verboseLevel(5);
54 //Trace::trackingMode=true;
55 //numberOfThreads(1);
56
57 //test on a sphere
58 Mesh meshd(Sphere(_center=Point(0.,0.,0.),_radius=1.,_nnodes=9,_domain_name="Omega"),_triangle,1,_subdiv);
59 Domain omega=meshd.domain("Omega");
60 Space W(omega,P0,"V",false);
61 Unknown u(W,"u"); TestFunction v(u,"v");
62 Vector<Real> x(W.dimSpace(),1.);
63 Kernel Gl=Laplace3dKernel();
64 IntegrationMethods ims(Sauter_Schwab,5,0.,_defaultRule,5,2.,_defaultRule,4,4.,_defaultRule,3);
65
66 Number nbp=50;
67
68 out<<"----------------- Laplace kernel single layer, dense matrix computation -------------------"<<eol;
69 std::cout<<"----------------- Laplace kernel single layer, dense matrix computation -------------------"<<eol;
70 BilinearForm blf=intg(omega,omega,u*Gl*v,ims,_noSymmetry);
71 elapsedTime();
72 TermMatrix A(blf,_denseRowStorage, "A");
73 elapsedTime("build dense matrix");
74 LargeMatrix<Real>& Alm=A.getLargeMatrix<Real>();
75 Vector<Real> Ax;
76 elapsedTime();
77 for(Number k=0; k<nbp; ++k) Ax=Alm*x;
78 out<<"dense matrix ("<<Alm.nbNonZero()<<") |A*x| = "<<norm(Ax)<<eol;
79 elapsedTime("dense matrix x vector");
80
81 Number nbox=10;
82 ClusterTree<FeDof> ct(W.feSpace()->dofs,_cardinalityBisection,nbox);
83 //ct.print(out);
84 elapsedTime("create Cluster tree");
85
86 HMatrix<Real,FeDof> hm(ct,ct,nbox,nbox);
87 elapsedTime("create HMatrix");
88 //ct.saveToFile("cluster_sphere.txt");
89 //hm.saveStructureToFile("hmatrix.dat");
90 hm.load(Alm); //load HMatrix from LargeMatrix
91 elapsedTime("load HMatrix from LargeMatrix");
92
93 //test matrix vector product
94 //==========================
95 elapsedTime();
96 Vector<Real> amx;
97 for(number_t k=0; k<50; k++) amx=Alm*x;
98 elapsedTime("A*x");
99 HMatrixNode<Real,FeDof>::initCounter();
100 Vector<Real> hmx;
101 for(number_t k=0; k<5; k++) hmx = hm*x;
102 elapsedTime("hm*x");
103 HMatrixNode<Real,FeDof>::stopCounter();
104 out<<"|amx-hmx| = "<<norm(amx-hmx)<<" matrixVectorProductCount = "<<HMatrixNode<Real,FeDof>::counter<<eol;
105 //test get leaves
106 std::list<HMatrixNode<Real,FeDof>* > leaves=hm.getLeaves();
107 std::list<HMatrixNode<Real,FeDof>* >::iterator itl=leaves.begin();
108 out<<" number of leaves = "<<leaves.size()<<eol;
109 //for(;itl!=leaves.end();++itl) (*itl)->printNode(out);
110
111 Number nt=20;
112 Number unsafe=0;
113 std::list<HMatrixNode<Real,FeDof>* > oleaves=hm.getLeaves(_row,nt,unsafe);
114 out<<"list of leaves with nt="<<nt<<" unsafe = "<<unsafe
115 <<" number of safe packets = "<<(oleaves.size()-unsafe)/nt<<eol;
116 itl=oleaves.begin();
117 //for(;itl!=oleaves.end();++itl) (*itl)->printNode(out);
118
119 out<<"norm2(am)="<<norm2(Alm)<<" norm2(hm)="<<hm.norm2()<<eol;
120
121 //test HMatrixNode stuff
122 //======================
123 // _clearfp();
124 // unsigned ucw = 0;
125 // _controlfp_s(&ucw, 0, _EM_INVALID|_EM_OVERFLOW); //activate Nan and overflow sigtrap
126 verboseLevel(20);
127
128 Vector<Real>::iterator ity;
129 HMatrixNode<Real,FeDof>* hno=hm.child(5);
130 out<<"----------------- HNode product -------------------"<<eol;
131 //hno->print(out);
132 Vector<Real> y0=Ax, y=Ax;
133 y0*=0; hno->multMatrixVector(x,y0);
134 // out<<"hno*x=";
135 // ity=y0.begin();
136 // for(number_t k=0;k<y0.size();++k, ++ity) if(*ity!=0) out<<" "<<k<<"->"<<*ity;
137 out<<"norm2(hno*x)="<<norm2(y0)<<eol;
138 out<<eol;
139 out<<"----------------- HNode hierarchical to LargeMatrix -------------------"<<eol;
140 hno->toLargeMatrix();
141 //hno->print(out);
142 y*=0; hno->multMatrixVector(x,y);
143 out<<"norm2(hno*x-Ax)="<<norm2(y-y0)<<eol;
144 out<<"----------------- HNode LargeMatrix to Hierarchical -------------------"<<eol;
145 hno->toHierarchical();
146 //hno->print(out);
147 y*=0; hno->multMatrixVector(x,y);
148 out<<"norm2(hno*x-Ax)="<<norm2(y-y0)<<eol;
149 out<<"----------------- HNode hierarchical to LowRankMatrix -------------------"<<eol;
150 hno->toApproximateMatrix(_lowRankApproximation,0,1.E-04);
151 //hno->print(out);
152 y*=0; hno->multMatrixVector(x,y);
153 out<<"norm2(hno*x-Ax)="<<norm2(y-y0)<<eol;
154 out<<"----------------- HNode LowRankMatrix to hierarchical -------------------"<<eol;
155 hno->toHierarchical();
156 //hno->print(out);
157 y*=0; hno->multMatrixVector(x,y);
158 out<<"norm2(hno*x-Ax)="<<norm2(y-y0)<<eol;
159
160 //test building Hmatrix from bilinear form
161 //========================================
162 Number sb=10;
163 HMatrixIM him0(_cardinalityBisection, sb,sb,_noHMApproximation,0.,ims);
164
165 // no approximation
166 out<<"----------------- Laplace kernel single layer, Hmatrix computation -------------------"<<eol;
167 theCout<<"----------------- Laplace kernel single layer, Hmatrix computation -------------------"<<eol;
168 BilinearForm blfhm0=intg(omega,omega,u*Gl*v,him0);
169 elapsedTime();
170 TermMatrix Ahm0(blfhm0,"Ahm0");
171 elapsedTime("build HMatrix exact");
172 HMatrix<real_t,FeDof>& hmA0=Ahm0.getHMatrix<Real>();
173 Vector<Real> Ahm0x;
174 elapsedTime();
175 for(Number k=0; k<nbp; ++k) Ahm0x=hmA0*x;
176 elapsedTime("exact Hmatrix x vector");
177 out<<"exact HM "; hmA0.printSummary(out); out<<" |A*x- hmA0*x| = "<<norm(Ax-Ahm0x)<<eol;
178
179 //test compression methods
180 DofCluster rowCluster=him0.rowCluster(), colCluster=him0.colCluster();
181 HMatrixIM him_svd(_svdCompression, 0.000001, ims, rowCluster, colCluster);
182 BilinearForm blfhm_svd=intg(omega,omega,u*Gl*v,him_svd);
183 elapsedTime();
184 TermMatrix Ahm_svd(blfhm_svd,"Ahm_svd");
185 elapsedTime("build HMatrix svd compression");
186 HMatrix<real_t,FeDof>& hmA_svd=Ahm_svd.getHMatrix<Real>();
187 Vector<Real> Ahm_svdx;
188 elapsedTime();
189 for(Number k=0; k<nbp; ++k) Ahm_svdx=hmA_svd*x;
190 elapsedTime("exact Hmatrix svd x vector");
191 out<<"svd HM "; hmA_svd.printSummary(out); out<<" |A*x- Ahm_svd*x| = "<<norm(Ax-Ahm_svdx)<<eol;
192
193 //test with a combination of a HMatrix and a LargeMatrix
194 out<<"----------------- Laplace kernel double layer - Id, dense matrix computation -------------------"<<eol;
195 theCout<<"----------------- Laplace kernel double layer - Id, dense matrix computation -------------------"<<eol;
196 BilinearForm blm=intg(omega,omega,u*ndotgrad_y(Gl)*v,ims)-0.5*intg(omega,u*v);;
197 TermMatrix B(blm,"B");
198 elapsedTime("build dense matrix");
199 Vector<Real> Bx;
200 LargeMatrix<Real>& Blm=B.getLargeMatrix<Real>();
201 elapsedTime();
202 for(Number k=0; k<nbp; ++k) Bx=Blm*x;
203 out<<"|B*x| = "<<norm(Bx)<<eol;
204 elapsedTime("dense matrix x vector");
205 out<<"----------------- Laplace kernel double layer - Id, Hmatrix computation -------------------"<<eol;
206 theCout<<"----------------- Laplace kernel double layer - Id, Hmatrix computation -------------------"<<eol;
207 BilinearForm blm0=intg(omega,omega,u*ndotgrad_y(Gl)*v,him0)-0.5*intg(omega,u*v);
208 TermMatrix B0(blm0, "B0");
209 elapsedTime("build exact HM/LM matrix");
210 Vector<Real> B0x;
211 HMatrix<real_t,FeDof>& B0hm=B0.getHMatrix<Real>();
212 elapsedTime();
213 for(Number k=0; k<nbp; ++k) B0x=B0hm*x;
214 out<<"exact HM/LM |B*x-B0*x| = "<<norm(Bx-B0x)<<eol;
215 elapsedTime("exact HM/LM matrix x vector");
216
217 //test for Helmholtz Kernel
218 out<<"----------------- Helmholtz kernel single layer, dense matrix computation -------------------"<<eol;
219 theCout<<"----------------- Helmholtz kernel single layer, dense matrix computation -------------------"<<eol;
220 Kernel Gh=Helmholtz3dKernel(2.);
221 Vector<Complex> z(W.dimSpace(),Complex(1.));
222 BilinearForm cblm=intg(omega,omega,u*Gh*v,ims);
223 TermMatrix C(cblm,"B");
224 elapsedTime("build full matrix");
225 Vector<Complex> Cx;
226 LargeMatrix<Complex>& Clm=C.getLargeMatrix<Complex>();
227 elapsedTime();
228 for(Number k=0; k<nbp; ++k) Cx=Clm*z;
229 out<<"|C*x| = "<<norm(Cx)<<eol;
230 elapsedTime("full matrix x vector");
231
232 // no approximation
233 out<<"----------------- Helmholtz kernel single layer, Hmatrix computation -------------------"<<eol;
234 theCout<<"----------------- Helmholtz kernel single layer, Hmatrix computation -------------------"<<eol;
235 BilinearForm cblfhm0=intg(omega,omega,u*Gh*v,him0);
236 elapsedTime();
237 TermMatrix Chm0(cblfhm0,"Chm0");
238 elapsedTime("build HMatrix exact");
239 HMatrix<Complex,FeDof>& hmC0=Chm0.getHMatrix<Complex>();
240 Vector<Complex> Chm0x;
241 elapsedTime();
242 for(Number k=0; k<nbp; ++k) Chm0x=hmC0*z;
243 elapsedTime("exact Hmatrix x vector");
244 out<<"exact HM/LM |C*x- hmC0*x| = "<<norm(Cx-Chm0x)<<eol;
245
246 //------------------------------------------------------------------------------------
247 // save results in a file or compare results with some references value in a file
248 //------------------------------------------------------------------------------------
249 trace_p->pop();
250 if(check) { return diffResFile(out, rootname); }
251 else { return saveResToFile(out, rootname); }
252 }
253
254 }
255