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_CgSolver.cpp
19 \author ManhHa NGUYEN
20 \since 13 Sep 2012
21 \date 29 Oct 2012
22
23 Low level tests of CgSolver.
24 Only tested with symmetric matrices. 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 information in a string
28 Only symmetric positive-definite matrix can be used as a preconditioner
29 http://en.wikipedia.org/wiki/Conjugate_gradient_method
30 Test order:
31 1. Real matrix and real vector
32 1.1. Without a preconditioner
33 a. DENSE STORAGE
34 b. CS STORAGE
35 c. SKYLINE STORAGE
36 1.2 With a preconditioner
37 a. DENSE STORAGE
38 b. CS STORAGE
39 c. SKYLINE STORAGE
40 2. Complex matrix and complex vector
41 The same as Real
42 */
43
44 #include "xlife++-libs.h"
45 #include "testUtils.hpp"
46
47 #include <iostream>
48 #include <fstream>
49 #include <vector>
50
51 using namespace xlifepp;
52
53 namespace unit_CgSolver
54 {
55
unit_CgSolver(bool check)56 String unit_CgSolver(bool check)
57 {
58 String rootname = "unit_CgSolver";
59 trace_p->push(rootname);
60 std::stringstream out; // string stream receiving results
61 out.precision(testPrec);
62
63 verboseLevel(3);
64
65 const int rowNum = 3;
66 const int colNum = 3;
67
68 const std::string rMatrixDataSym("matrix3x3Sym.data");
69 const std::string rMatrixDataSkewSym("matrix3x3SkewSym.data");
70 const std::string rMatrixDataNoSym("matrix3x3NoSym.data");
71 const std::string rMatrixDataSymPosDef("matrix3x3SymPosDef.data");
72
73 const std::string cMatrixDataSym("cmatrix3x3Sym.data");
74 const std::string cMatrixDataNoSym("cmatrix3x3NoSym.data");
75 const std::string cMatrixDataSymSelfAjoint("cmatrix3x3SymSelfAjoint.data");
76 const std::string cMatrixDataSymSkewAjoint("cmatrix3x3SymSkewAjoint.data");
77 const std::string cMatrixDataSymPosDef("cmatrix3x3SymPosDef.data");
78
79 CgSolver cgSolverRow, cgSolverCol, cgSolverDual, cgSolverSym;
80 CgSolver cgSolverRowCplx, cgSolverColCplx, cgSolverDualCplx, cgSolverSymCplx;
81
82 //////////////////////////////////////////////////////////////////////////
83 //
84 // Vectors
85 //
86 /////////////////////////////////////////////////////////////////////////
87 // Real Vector B
88 Vector<Real> rvB(rowNum, 1.);
89 //out << "Real vector B is: " << rvB << std::endl;
90
91 // Real initial value
92 Vector<Real> rvX0(rowNum, 0.);
93 out << "Real vector x0 is: " << rvX0 << std::endl;
94
95 // Vector real unknown
96 Vector<Real> rvX(rowNum);
97 Vector<Real> rvXe(rowNum);
98 for(number_t i=0;i<rowNum; i++) rvXe(i+1)=Real(i+1);
99 out << "Real vector Xe is: " << rvXe << std::endl;
100
101 //-------------------------
102 // Vector complex B
103 Vector<Complex> cvB(rowNum, 1.);
104 //out << "The Complex vector B is: " << cvB << std::endl;
105
106 // Vector initial value
107 Vector<Complex> cvX0(rowNum, 0.);
108 out << "The Complex vector x0 is: " << cvX0 << std::endl;
109
110 // Vector complex unknown
111 Vector<Complex> cvX(rowNum);
112 Vector<Complex> cvXe(rowNum);
113 for(number_t i=0;i<rowNum; i++) cvXe(i+1)=Complex(1.,1.)*Real(i+1);
114 out << "Complex vector Xe is: " << cvXe << std::endl;
115
116 //////////////////////////////////////////////////////////////////////////
117 //
118 // Real Large Matrix (Dense Type)
119 //
120 /////////////////////////////////////////////////////////////////////////
121 // Symmetric
122 LargeMatrix<Real> rMatDenseRowSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _dense, _row);
123 LargeMatrix<Real> rMatDenseColSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _dense, _col);
124 LargeMatrix<Real> rMatDenseDualSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _dense, _dual);
125 LargeMatrix<Real> rMatDenseSymSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, _dense, _symmetric);
126
127 //////////////////////////////////////////////////////////////////////////
128 //
129 // Complex Large Matrix (Dense Type)
130 //
131 /////////////////////////////////////////////////////////////////////////
132 // Symmetric
133 LargeMatrix<Complex> cMatDenseRowSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _dense, _row);
134 LargeMatrix<Complex> cMatDenseColSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _dense, _col);
135 LargeMatrix<Complex> cMatDenseDualSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _dense, _dual);
136 LargeMatrix<Complex> cMatDenseSymSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, _dense, _symmetric);
137
138 //////////////////////////////////////////////////////////////////////////
139 //
140 // Real Large Matrix (CS Type)
141 //
142 /////////////////////////////////////////////////////////////////////////
143 // Symmetric
144 LargeMatrix<Real> rMatCsRowSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _cs, _row);
145 LargeMatrix<Real> rMatCsColSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _cs, _col);
146 LargeMatrix<Real> rMatCsDualSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _cs, _dual);
147 LargeMatrix<Real> rMatCsSymSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, _cs, _symmetric);
148
149 //////////////////////////////////////////////////////////////////////////
150 //
151 // Complex Large Matrix (Cs Type)
152 //
153 /////////////////////////////////////////////////////////////////////////
154 // Symmetric
155 LargeMatrix<Complex> cMatCsRowSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _cs, _row);
156 LargeMatrix<Complex> cMatCsColSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _cs, _col);
157 LargeMatrix<Complex> cMatCsDualSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _cs, _dual);
158 LargeMatrix<Complex> cMatCsSymSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, _cs, _symmetric);
159
160 //////////////////////////////////////////////////////////////////////////
161 //
162 // Real Large Matrix (Skyline Type)
163 //
164 /////////////////////////////////////////////////////////////////////////
165 // Symmetric
166 LargeMatrix<Real> rMatSkylineDualSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _skyline, _dual);
167 LargeMatrix<Real> rMatSkylineSymSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, _skyline, _symmetric);
168
169 //////////////////////////////////////////////////////////////////////////
170 //
171 // Complex Large Matrix (Skyline Type)
172 //
173 /////////////////////////////////////////////////////////////////////////
174 // Symmetric
175 LargeMatrix<Complex> cMatSkylineDualSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _skyline, _dual);
176 LargeMatrix<Complex> cMatSkylineSymSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, _skyline, _symmetric);
177
178 //////////////////////////////////////////////////////////////////////////
179 //
180 // Real Large Matrix factorized and used as a preconditioner
181 //
182 /////////////////////////////////////////////////////////////////////////
183 // LargeMatrix<Real> rMatSkylineSymSymLdlt(inputsPathTo(rMatrixDataSymPosDef), _dense, rowNum, _skyline, _symmetric);
184 // ldltFactorize(rMatSkylineSymSymLdlt);
185 // Preconditioner<LargeMatrix<Real> > pcLdltSymSkylineSym(rMatSkylineSymSymLdlt, _ldltPrec, 1.5);
186 //
187 // LargeMatrix<Real> rMatCsSymSymPc(inputsPathTo(rMatrixDataSym), _dense, rowNum, _cs, _symmetric);
188 // Preconditioner<LargeMatrix<Real> > pcSorCsSym(rMatCsSymSymPc, _ssorPrec, 1.5);
189 //
190 // //////////////////////////////////////////////////////////////////////////
191 // //
192 // // Complex Large Matrix to factorize (Skyline Type) as a preconditioner
193 // //
194 // /////////////////////////////////////////////////////////////////////////
195 // LargeMatrix<Complex> cMatSkylineSymSymLdlStar(inputsPathTo(cMatrixDataSymPosDef), _dense, rowNum, _skyline, _selfAdjoint);
196 // ldlstarFactorize(cMatSkylineSymSymLdlStar);
197 // Preconditioner<LargeMatrix<Complex> > pcLdlStar(cMatSkylineSymSymLdlStar, _ldlstarPrec, 1.7);
198
199 /////////////////////////////////////////////////////////////////////////////
200 //-----------------------------------------------------------------------------
201 // A. Test with Real value
202 //-----------------------------------------------------------------------------
203 // I. Solver without a preconditioner
204 //
205 //----------------------------------
206 //---------------------------------
207 // Test with DENSE STORAGE
208 // Test with symmetric matrices
209 rvB = rMatDenseRowSym *rvXe;
210 out << "Real vector B is: " << rvB << std::endl;
211 rvX = cgSolverRow(rMatDenseRowSym, rvB, rvX0, _real);
212 out << "The f CgSolver result with real dense row sym matrix is: " << rvX << std::endl;
213 rvX = cgSolverCol(rMatDenseColSym, rvB, rvX0, _real);
214 out << "The f CgSolver result with real dense col sym matrix is: " << rvX << std::endl;
215 rvX = cgSolverDual(rMatDenseDualSym, rvB, rvX0, _real);
216 out << "The f CgSolver result with real dense dual sym matrix is: " << rvX << std::endl;
217 rvX = cgSolverSym(rMatDenseSymSym, rvB, rvX0, _real);
218 out << "The f CgSolver result with real dense sym sym matrix is: " << rvX << std::endl << std::endl;
219
220 // Test with CS STORAGE
221 // Test with symmetric matrices
222 rvB = rMatCsRowSym *rvXe;
223 out << "Real vector B is: " << rvB << std::endl;
224 rvX = cgSolverRow(rMatCsRowSym, rvB, rvX0, _real);
225 out << "The f CgSolver result with real Cs row sym matrix is: " << rvX << std::endl;
226 rvX = cgSolverCol(rMatCsColSym, rvB, rvX0, _real);
227 out << "The f CgSolver result with real Cs col sym matrix is: " << rvX << std::endl;
228 rvX = cgSolverDual(rMatCsDualSym, rvB, rvX0, _real);
229 out << "The f CgSolver result with real Cs dual sym matrix is: " << rvX << std::endl;
230 rvX = cgSolverSym(rMatCsSymSym, rvB, rvX0, _real);
231 out << "The f CgSolver result with real Cs sym sym matrix is: " << rvX << std::endl << std::endl;
232
233 //----------------------------------
234 // Test with SKYLINE STORAGE
235 // Test with symmetric matrices
236 rvB = rMatSkylineDualSym *rvXe;
237 out << "Real vector B is: " << rvB << std::endl;
238 rvX = cgSolverDual(rMatSkylineDualSym, rvB, rvX0, _real);
239 out << "The f CgSolver result with real Skyline dual sym matrix is: " << rvX << std::endl;
240 rvX = cgSolverSym(rMatSkylineSymSym, rvB, rvX0, _real);
241 out << "The f CgSolver result with real Skyline sym sym matrix is: " << rvX << std::endl << std::endl;
242 /*
243 // II. Solver with a preconditioner
244 //-------------------------------
245 // Test with DENSE STORAGE
246 // Test with symmetric matrices
247 rvX = cgSolverRow(rMatDenseRowSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
248 out << "The f CgSolver result with real Dense row matrix and a ldlt factorized preconditioner is: " << rvX << std::endl;
249 rvX = cgSolverCol(rMatDenseColSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
250 out << "The f CgSolver result with real dense col matrix and a ldlt factorized preconditioner is: " << rvX << std::endl;
251 rvX = cgSolverDual(rMatDenseDualSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
252 out << "The f CgSolver result with real dense dual matrix and a ldlt factorized preconditioner is: " << rvX << std::endl;
253 rvX = cgSolverSym(rMatDenseSymSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
254 out << "The f CgSolver result with real dense sym matrix and a a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
255
256 rvX = cgSolverRow(rMatDenseRowSym, rvB, rvX0, pcSorCsSym, _real);
257 out << "The f CgSolver result with real Dense row matrix and a sor preconditioner is: " << rvX << std::endl;
258 rvX = cgSolverCol(rMatDenseColSym, rvB, rvX0, pcSorCsSym, _real);
259 out << "The f CgSolver result with real dense col matrix and a sor preconditioner is: " << rvX << std::endl;
260 rvX = cgSolverDual(rMatDenseDualSym, rvB, rvX0, pcSorCsSym, _real);
261 out << "The f CgSolver result with real dense dual matrix and a sor preconditioner is: " << rvX << std::endl;
262 rvX = cgSolverSym(rMatDenseSymSym, rvB, rvX0, pcSorCsSym, _real);
263 out << "The f CgSolver result with real dense sym matrix and a a sor preconditioner is: " << rvX << std::endl << std::endl;
264
265 // Test with CS STORAGE
266 // Test with symmetric matrices
267 rvX = cgSolverRow(rMatCsRowSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
268 out << "The f CgSolver result with real Cs row matrix and a ldlt factorized preconditioner is: " << rvX << std::endl;
269 rvX = cgSolverCol(rMatCsColSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
270 out << "The f CgSolver result with real Cs col matrix and a ldlt factorized preconditioner is: " << rvX << std::endl;
271 rvX = cgSolverDual(rMatCsDualSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
272 out << "The f CgSolver result with real Cs dual matrix and a ldlt factorized preconditioner is: " << rvX << std::endl;
273 rvX = cgSolverSym(rMatCsSymSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
274 out << "The f CgSolver result with real Cs sym matrix and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
275
276 rvX = cgSolverRow(rMatCsRowSym, rvB, rvX0, pcSorCsSym, _real);
277 out << "The f CgSolver result with real Cs row matrix and a sor preconditioner is: " << rvX << std::endl;
278 rvX = cgSolverCol(rMatCsColSym, rvB, rvX0, pcSorCsSym, _real);
279 out << "The f CgSolver result with real Cs col matrix and a sor preconditioner is: " << rvX << std::endl;
280 rvX = cgSolverDual(rMatCsDualSym, rvB, rvX0, pcSorCsSym, _real);
281 out << "The f CgSolver result with real Cs dual matrix and a sor preconditioner is: " << rvX << std::endl;
282 rvX = cgSolverSym(rMatCsSymSym, rvB, rvX0, pcSorCsSym, _real);
283 out << "The f CgSolver result with real Cs sym matrix and a sor preconditioner is: " << rvX << std::endl << std::endl;
284
285 //----------------------------------
286 // Test with SKYLINE STORAGE
287 // Test with symmetric matrices
288 rvX = cgSolverDual(rMatSkylineDualSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
289 out << "The f CgSolver result with real Skyline dual matrix and a ldlt factorized preconditioner is: " << rvX << std::endl;
290 rvX = cgSolverSym(rMatSkylineSymSym, rvB, rvX0, pcLdltSymSkylineSym, _real);
291 out << "The f CgSolver result with real Skyline sym matrix and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
292
293 rvX = cgSolverDual(rMatSkylineDualSym, rvB, rvX0, pcSorCsSym, _real);
294 out << "The f CgSolver result with real Skyline dual matrix and a sor preconditioner is: " << rvX << std::endl;
295 rvX = cgSolverSym(rMatSkylineSymSym, rvB, rvX0, pcSorCsSym, _real);
296 out << "The f CgSolver result with real Skyline sym matrix and a sor preconditioner is: " << rvX << std::endl << std::endl;
297 */
298 //-----------------------------------------------------------------------------
299 // B. Test with Complex values (Complex x Complex = Complex)
300 //-----------------------------------------------------------------------------
301 // I. Solver without a preconditioner
302 // Test with DENSE STORAGE
303 // Test with symmetric matrices
304 cvB = cMatDenseRowSym *cvXe;
305 out << "Complex vector B is: " << cvB << std::endl;
306 cvX = cgSolverRowCplx(cMatDenseRowSym, cvB, cvX0, xlifepp::_complex);
307 out << "The f CgSolver result with complex dense row is: " << cvX << std::endl;
308 cvX = cgSolverColCplx(cMatDenseColSym, cvB, cvX0, xlifepp::_complex);
309 out << "The f CgSolver result with complex dense col is: " << cvX << std::endl;
310 cvX = cgSolverDualCplx(cMatDenseDualSym, cvB, cvX0, xlifepp::_complex);
311 out << "The f CgSolver result with complex dense dual is: " << cvX << std::endl;
312 cvX = cgSolverSymCplx(cMatDenseSymSym, cvB, cvX0, xlifepp::_complex);
313 out << "The f CgSolver result with complex dense sym is: " << cvX << std::endl << std::endl;
314
315 //---------------------------------
316 // Test with CS STORAGE
317 // Test with symmetric matrices
318 cvB = cMatCsRowSym *cvXe;
319 out << "Complex vector B is: " << cvB << std::endl;
320 cvX = cgSolverRowCplx(cMatCsRowSym, cvB, cvX0, xlifepp::_complex);
321 out << "The f CgSolver result with complex Cs row is: " << cvX << std::endl;
322 cvX = cgSolverColCplx(cMatCsColSym, cvB, cvX0, xlifepp::_complex);
323 out << "The f CgSolver result with complex Cs col is: " << cvX << std::endl;
324 cvX = cgSolverDualCplx(cMatCsDualSym, cvB, cvX0, xlifepp::_complex);
325 out << "The f CgSolver result with complex Cs dual is: " << cvX << std::endl;
326 cvX = cgSolverSymCplx(cMatCsSymSym, cvB, cvX0, xlifepp::_complex);
327 out << "The f CgSolver result with complex Cs sym is: " << cvX << std::endl << std::endl;
328
329 //--------------------------------
330 // Test with SKYLINE STORAGE
331 // Test with symmetric matrices
332 cvB = cMatSkylineDualSym *cvXe;
333 out << "Complex vector B is: " << cvB << std::endl;
334 cvX = cgSolverDualCplx(cMatSkylineDualSym, cvB, cvX0, xlifepp::_complex);
335 out << "The f CgSolver result with complex Skyline dual is: " << cvX << std::endl;
336 cvX = cgSolverSymCplx(cMatSkylineSymSym, cvB, cvX0, xlifepp::_complex);
337 out << "The f CgSolver result with complex Skyline sym is: " << cvX << std::endl << std::endl;
338 /*
339 // For the complex solver, haven't found a "good" preconditioner to return
340 // a correct result :(
341 // II. Solver with a preconditioner
342 // Test with DENSE STORAGE
343 // Test with symmetric matrices
344 cvX = cgSolverRowCplx(cMatDenseRowSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
345 out << "The f CgSolver result with complex dense row and a ldl* factorized preconditioner is: " << cvX << std::endl;
346 cvX = cgSolverColCplx(cMatDenseColSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
347 out << "The f CgSolver result with complex dense col and a ldl* factorized preconditioner is: " << cvX << std::endl;
348 cvX = cgSolverDualCplx(cMatDenseDualSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
349 out << "The f CgSolver result with complex dense dual and a ldl* factorized preconditioner is: " << cvX << std::endl;
350 cvX = cgSolverSymCplx(cMatDenseSymSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
351 out << "The f CgSolver result with complex dense sym and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
352
353 //---------------------------------
354 // Test with CS STORAGE
355 // Test with symmetric matrices
356 cvX = cgSolverRowCplx(cMatCsRowSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
357 out << "The f CgSolver result with complex Cs row and a ldl* factorized preconditioner is: " << cvX << std::endl;
358 cvX = cgSolverColCplx(cMatCsColSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
359 out << "The f CgSolver result with complex Cs col and a ldl* factorized preconditioner is: " << cvX << std::endl;
360 cvX = cgSolverDualCplx(cMatCsDualSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
361 out << "The f CgSolver result with complex Cs dual and a ldl* factorized preconditioner is: " << cvX << std::endl;
362 cvX = cgSolverSymCplx(cMatCsSymSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
363 out << "The f CgSolver result with complex Cs sym and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
364
365 //--------------------------------
366 // Test with SKYLINE STORAGE
367 // Test with symmetric matrices
368 cvX = cgSolverDualCplx(cMatSkylineDualSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
369 out << "The f CgSolver result with complex Skyline dual and a ldl* factorized preconditioner is: " << cvX << std::endl;
370 cvX = cgSolverSymCplx(cMatSkylineSymSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
371 out << "The f CgSolver result with complex Skyline sym and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
372 */
373 //------------------------------------------------------------------------------------
374 // save results in a file or compare results with some references value in a file
375 //------------------------------------------------------------------------------------
376 trace_p->pop();
377 if (check) { return diffResFile(out, rootname); }
378 else { return saveResToFile(out, rootname); }
379 }
380
381 }
382