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_SsorSolver.cpp
19    \author ManhHa NGUYEN
20    \since 20 Oct 2012
21    \date 21 Oct 2012
22 
23     Low level tests of SssorSolver.
24     Only tested with CS storage. Almost functionalities are checked.
25     The correctness of this solver depends so much on the value of omega
26     This function may either creates a reference file storing the results (check=false)
27     or compares results to those stored in the reference file (check=true)
28     Test order:
29     1. Real matrix and real vector
30     2. Complex matrix and complex vector
31 */
32 
33 #include "xlife++-libs.h"
34 #include "testUtils.hpp"
35 
36 #include <iostream>
37 #include <fstream>
38 #include <vector>
39 
40 using namespace xlifepp;
41 
42 namespace unit_SsorSolver {
43 
unit_SsorSolver(bool check)44 String unit_SsorSolver(bool check)
45 {
46   String rootname = "unit_SsorSolver";
47   trace_p->push(rootname);
48   std::stringstream out; // string stream receiving results
49   out.precision(testPrec);
50 
51   verboseLevel(3);
52 
53   const int rowNum = 3;
54   const int colNum = 3;
55   const Real omega = 1.08;
56 
57   const std::string rMatrixDataSym("matrix3x3Sym.data");
58   const std::string rMatrixDataSkewSym("matrix3x3SkewSym.data");
59   const std::string rMatrixDataNoSym("matrix3x3NoSym.data");
60   const std::string rMatrixDataSymPosDef("matrix3x3SymPosDef.data");
61 
62   const std::string cMatrixDataSym("cmatrix3x3Sym.data");
63   const std::string cMatrixDataNoSym("cmatrix3x3NoSym.data");
64   const std::string cMatrixDataSymSelfAjoint("cmatrix3x3SymSelfAjoint.data");
65   const std::string cMatrixDataSymSkewAjoint("cmatrix3x3SymSkewAjoint.data");
66   const std::string cMatrixDataSymPosDef("cmatrix3x3SymPosDef.data");
67 
68   SsorSolver ssorSolver(omega, 1.e-10, 1000);
69 
70   //////////////////////////////////////////////////////////////////////////
71   //  Vectors
72   /////////////////////////////////////////////////////////////////////////
73  // Real Vector B
74   Vector<Real> rvB(rowNum, 1.);
75 
76   // Real initial value
77   Vector<Real> rvX0(rowNum, 0.);
78   out << "Real vector x0 is: " << rvX0 << std::endl;
79 
80   // Vector real unknown
81   Vector<Real> rvX(rowNum);
82   Vector<Real> rvXe(rowNum);
83   for(number_t i=0;i<rowNum; i++) rvXe(i+1)=Real(i+1);
84   out << "Real vector Xe is: " << rvXe <<  std::endl;
85 
86   //-------------------------
87   // Vector complex B
88   Vector<Complex> cvB(rowNum, 1.);
89   //out << "The Complex vector B is: " << cvB <<  std::endl;
90 
91   //  Vector initial value
92   Vector<Complex> cvX0(rowNum, 0.);
93   out << "The Complex vector x0 is: " << cvX0 << std::endl;
94 
95   // Vector complex unknown
96   Vector<Complex> cvX(rowNum);
97   Vector<Complex> cvXe(rowNum);
98   for(number_t i=0;i<rowNum; i++) cvXe(i+1)=Complex(1.,1.)*Real(i+1);
99   out << "Complex vector Xe is: " << cvXe <<  std::endl;
100 
101   //////////////////////////////////////////////////////////////////////////
102   //
103   //  Real Large Matrix (Dense Type)
104   //
105   /////////////////////////////////////////////////////////////////////////
106   // No symmetric
107   LargeMatrix<Real> rMatDenseRow(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _dense, _row);
108   LargeMatrix<Real> rMatDenseCol(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _dense, _col);
109   LargeMatrix<Real> rMatDenseDual(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _dense, _dual);
110   LargeMatrix<Real> rMatDenseSym(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _dense, _sym);
111   // Symmetric
112   LargeMatrix<Real> rMatDenseRowSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _dense, _row);
113   LargeMatrix<Real> rMatDenseColSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _dense, _col);
114   LargeMatrix<Real> rMatDenseDualSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _dense, _dual);
115   LargeMatrix<Real> rMatDenseSymSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, _dense, _symmetric);
116   // Skew Symmetric
117   LargeMatrix<Real> rMatDenseRowSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _dense, _row);
118   LargeMatrix<Real> rMatDenseColSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _dense, _col);
119   LargeMatrix<Real> rMatDenseDualSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _dense, _dual);
120   LargeMatrix<Real> rMatDenseSymSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, _dense, _skewSymmetric);
121 
122   //////////////////////////////////////////////////////////////////////////
123   //
124   //  Complex Large Matrix (Dense Type)
125   //
126   /////////////////////////////////////////////////////////////////////////
127   // No symmetric
128   LargeMatrix<Complex> cMatDenseRow(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _dense, _row);
129   LargeMatrix<Complex> cMatDenseCol(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _dense, _col);
130   LargeMatrix<Complex> cMatDenseDual(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _dense, _dual);
131   LargeMatrix<Complex> cMatDenseSym(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _dense, _sym);
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   // Self Ajoint Symmetric
138   LargeMatrix<Complex> cMatDenseRowSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _dense, _row);
139   LargeMatrix<Complex> cMatDenseColSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _dense, _col);
140   LargeMatrix<Complex> cMatDenseDualSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _dense, _dual);
141   LargeMatrix<Complex> cMatDenseSymSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, _dense, _selfAdjoint);
142   // Skew Ajoint Symmetric
143   LargeMatrix<Complex> cMatDenseRowSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _dense, _row);
144   LargeMatrix<Complex> cMatDenseColSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _dense, _col);
145   LargeMatrix<Complex> cMatDenseDualSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _dense, _dual);
146   LargeMatrix<Complex> cMatDenseSymSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, _dense, _skewAdjoint);
147 
148   //////////////////////////////////////////////////////////////////////////
149   //
150   //  Real Large Matrix (CS Type)
151   //
152   /////////////////////////////////////////////////////////////////////////
153   // No symmetric
154   LargeMatrix<Real> rMatCsRow(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _cs, _row);
155   LargeMatrix<Real> rMatCsCol(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _cs, _col);
156   LargeMatrix<Real> rMatCsDual(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _cs, _dual);
157   LargeMatrix<Real> rMatCsSym(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _cs, _sym);
158   // Symmetric
159   LargeMatrix<Real> rMatCsRowSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _cs, _row);
160   LargeMatrix<Real> rMatCsColSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _cs, _col);
161   LargeMatrix<Real> rMatCsDualSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _cs, _dual);
162   LargeMatrix<Real> rMatCsSymSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, _cs, _symmetric);
163   // Skew Symmetric
164   LargeMatrix<Real> rMatCsRowSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _cs, _row);
165   LargeMatrix<Real> rMatCsColSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _cs, _col);
166   LargeMatrix<Real> rMatCsDualSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _cs, _dual);
167   LargeMatrix<Real> rMatCsSymSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, _cs, _skewSymmetric);
168 
169   //////////////////////////////////////////////////////////////////////////
170   //
171   //  Complex Large Matrix (Cs Type)
172   //
173   /////////////////////////////////////////////////////////////////////////
174   // No symmetric
175   LargeMatrix<Complex> cMatCsRow(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _cs, _row);
176   LargeMatrix<Complex> cMatCsCol(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _cs, _col);
177   LargeMatrix<Complex> cMatCsDual(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _cs, _dual);
178   LargeMatrix<Complex> cMatCsSym(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _cs, _sym);
179   // Symmetric
180   LargeMatrix<Complex> cMatCsRowSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _cs, _row);
181   LargeMatrix<Complex> cMatCsColSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _cs, _col);
182   LargeMatrix<Complex> cMatCsDualSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _cs, _dual);
183   LargeMatrix<Complex> cMatCsSymSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, _cs, _symmetric);
184   // Self Ajoint Symmetric
185   LargeMatrix<Complex> cMatCsRowSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _cs, _row);
186   LargeMatrix<Complex> cMatCsColSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _cs, _col);
187   LargeMatrix<Complex> cMatCsDualSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _cs, _dual);
188   LargeMatrix<Complex> cMatCsSymSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, _cs, _selfAdjoint);
189   // Skew Ajoint Symmetric
190   LargeMatrix<Complex> cMatCsRowSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _cs, _row);
191   LargeMatrix<Complex> cMatCsColSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _cs, _col);
192   LargeMatrix<Complex> cMatCsDualSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _cs, _dual);
193   LargeMatrix<Complex> cMatCsSymSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, _cs, _skewAdjoint);
194 
195   //////////////////////////////////////////////////////////////////////////
196   //
197   //  Real Large Matrix (Skyline Type)
198   //
199   /////////////////////////////////////////////////////////////////////////
200   // No symmetric
201   LargeMatrix<Real> rMatSkylineDual(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _dual);
202   LargeMatrix<Real> rMatSkylineSym(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _sym);
203   // Symmetric
204   LargeMatrix<Real> rMatSkylineDualSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _skyline, _dual);
205   LargeMatrix<Real> rMatSkylineSymSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, _skyline, _symmetric);
206   // Skew Symmetric
207   LargeMatrix<Real> rMatSkylineDualSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _skyline, _dual);
208   LargeMatrix<Real> rMatSkylineSymSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, _skyline, _skewSymmetric);
209 
210   //////////////////////////////////////////////////////////////////////////
211   //
212   //  Complex Large Matrix (Skyline Type)
213   //
214   /////////////////////////////////////////////////////////////////////////
215   // No symmetric
216   LargeMatrix<Complex> cMatSkylineDual(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _dual);
217   LargeMatrix<Complex> cMatSkylineSym(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _sym);
218   // Symmetric
219   LargeMatrix<Complex> cMatSkylineDualSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _skyline, _dual);
220   LargeMatrix<Complex> cMatSkylineSymSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, _skyline, _symmetric);
221   // Self Ajoint Symmetric
222   LargeMatrix<Complex> cMatSkylineDualSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _skyline, _dual);
223   LargeMatrix<Complex> cMatSkylineSymSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, _skyline, _selfAdjoint);
224   // Skew Ajoint Symmetric
225   LargeMatrix<Complex> cMatSkylineDualSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _skyline, _dual);
226   LargeMatrix<Complex> cMatSkylineSymSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, _skyline, _skewAdjoint);
227 
228   //-----------------------------------------------------------------------------
229   // A. Test with Real value
230   //-----------------------------------------------------------------------------
231   //
232   // Test with DENSE STORAGE
233   // Test with symmetric matrices
234   rvB = rMatDenseRowSym *rvXe;
235   out << "Real vector B is: " << rvB <<  std::endl;
236   rvX = ssorSolver(rMatDenseRowSym, rvB, rvX0);
237   out << "The SssrSolver result with real dense row and sym matrix is: " << rvX << std::endl;
238   ssorSolver.resetSolver();
239   rvX = ssorSolver(rMatDenseColSym, rvB, rvX0);
240   out << "The SssrSolver result with real dense col and sym matrix is: " << rvX << std::endl;
241   ssorSolver.resetSolver();
242   rvX = ssorSolver(rMatDenseDualSym, rvB, rvX0);
243   out << "The SssrSolver result with real dense dual and sym matrix is: " << rvX << std::endl;
244   ssorSolver.resetSolver();
245   rvX = ssorSolver(rMatDenseSymSym, rvB, rvX0);
246   out << "The SssrSolver result with real dense sym and sym matrix is: " << rvX << std::endl << std::endl;
247   ssorSolver.resetSolver();
248 
249   // Test with non symmetric matrices
250   rvB = rMatDenseRow *rvXe;
251   out << "Real vector B is: " << rvB <<  std::endl;
252   rvX = ssorSolver(rMatDenseRow, rvB, rvX0);
253   out << "The SssrSolver result with real dense row and non-sym matrix is: " << rvX << std::endl;
254   ssorSolver.resetSolver();
255   rvX = ssorSolver(rMatDenseCol, rvB, rvX0);
256   out << "The SssrSolver result with real dense col and non-sym matrix is: " << rvX << std::endl;
257   ssorSolver.resetSolver();
258   rvX = ssorSolver(rMatDenseDual, rvB, rvX0);
259   out << "The SssrSolver result with real dense dual and non-sym matrix is: " << rvX << std::endl;
260   ssorSolver.resetSolver();
261   rvX = ssorSolver(rMatDenseSym, rvB, rvX0);
262   out << "The SssrSolver result with real dense sym and non-sym matrix is: " << rvX << std::endl << std::endl;
263 
264   //---------------------------------
265   // Test with CS STORAGE
266   // Test with symmetric matrices
267   rvB = rMatCsRowSym *rvXe;
268   out << "Real vector B is: " << rvB <<  std::endl;
269   rvX = ssorSolver(rMatCsRowSym, rvB, rvX0);
270   out << "The SssrSolver result with real Cs row and sym matrix is: " << rvX << std::endl;
271   ssorSolver.resetSolver();
272   rvX = ssorSolver(rMatCsColSym, rvB, rvX0);
273   out << "The SssrSolver result with real Cs col and sym matrix is: " << rvX << std::endl;
274   ssorSolver.resetSolver();
275   rvX = ssorSolver(rMatCsDualSym, rvB, rvX0);
276   out << "The SssrSolver result with real Cs dual and sym matrix is: " << rvX << std::endl;
277   ssorSolver.resetSolver();
278   rvX = ssorSolver(rMatCsSymSym, rvB, rvX0);
279   out << "The SssrSolver result with real Cs sym and sym matrix is: " << rvX << std::endl << std::endl;
280   ssorSolver.resetSolver();
281 
282   // Test with non symmetric matrices
283    rvB = rMatCsRow *rvXe;
284   out << "Real vector B is: " << rvB <<  std::endl;
285   rvX = ssorSolver(rMatCsRow, rvB, rvX0);
286   out << "The SssrSolver result with real Cs row and non-sym matrix is: " << rvX << std::endl;
287   ssorSolver.resetSolver();
288   rvX = ssorSolver(rMatCsCol, rvB, rvX0);
289   out << "The SssrSolver result with real Cs col and non-sym matrix is: " << rvX << std::endl;
290   ssorSolver.resetSolver();
291   rvX = ssorSolver(rMatCsDual, rvB, rvX0);
292   out << "The SssrSolver result with real Cs dual and non-sym matrix is: " << rvX << std::endl;
293   ssorSolver.resetSolver();
294   rvX = ssorSolver(rMatCsSym, rvB, rvX0);
295   out << "The SssrSolver result with real Cs sym and non-sym matrix is: " << rvX << std::endl << std::endl;
296   ssorSolver.resetSolver();
297 
298   //------------------------------
299   // Test with SKYLINE STORAGE
300   // Test with symmetric matrices
301   rvB = rMatSkylineDualSym *rvXe;
302   out << "Real vector B is: " << rvB <<  std::endl;
303   rvX = ssorSolver(rMatSkylineDualSym, rvB, rvX0);
304   out << "The SssrSolver result with real Skyline dual and sym matrix is: " << rvX << std::endl;
305   ssorSolver.resetSolver();
306   rvX = ssorSolver(rMatSkylineSymSym, rvB, rvX0);
307   out << "The SssrSolver result with real Skyline sym and sym matrix is: " << rvX << std::endl << std::endl;
308   ssorSolver.resetSolver();
309 
310   // Test with non symmetric matrices
311   rvB = rMatSkylineDual *rvXe;
312   out << "Real vector B is: " << rvB <<  std::endl;
313   rvX = ssorSolver(rMatSkylineDual, rvB, rvX0);
314   out << "The SssrSolver result with real Skyline dual and non-sym matrix is: " << rvX << std::endl;
315   rvX = ssorSolver(rMatSkylineSym, rvB, rvX0);
316   out << "The SssrSolver result with real Skyline sym and non-sym matrix is: " << rvX << std::endl << std::endl;
317 
318   //-----------------------------------------------------------------------------
319   // B. Test with Complex values (Complex x Complex = Complex)
320   //-----------------------------------------------------------------------------
321   //---------------------------------
322   // Test with CS STORAGE
323   // Test with symmetric matrices
324   cvB = cMatCsSymSym *cvXe;
325   out << "Complex vector B is: " << cvB <<  std::endl;
326   cvX = ssorSolver(cMatCsSymSym, cvB, cvX0);
327   out << "The SssorSolver result with complex Cs sym is: " << cvX << std::endl << std::endl;
328 
329   //------------------------------------------------------------------------------------
330   // save results in a file or compare results with some references value in a file
331   //------------------------------------------------------------------------------------
332   trace_p->pop();
333   if (check) { return diffResFile(out, rootname); }
334   else { return saveResToFile(out, rootname); }
335 }
336 
337 }
338