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