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_GmresSolver.cpp
19     \author ManhHa NGUYEN
20     \since 13 Sep 2012
21     \date 06 Oct 2012
22 
23     Low level tests of GmresSolver.
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 information in a string
28 */
29 
30 #include "xlife++-libs.h"
31 #include "testUtils.hpp"
32 
33 #include <iostream>
34 #include <fstream>
35 #include <vector>
36 
37 using namespace xlifepp;
38 
39 namespace unit_GmresSolver {
40 
unit_GmresSolver(bool check)41 String unit_GmresSolver(bool check)
42 {
43   String rootname = "unit_GmresSolver";
44   trace_p->push(rootname);
45   std::stringstream out; // string stream receiving results
46   out.precision(testPrec);
47 
48   const int rowNum = 3;
49   const int colNum = 3;
50 
51   const std::string rMatrixDataSym("matrix3x3Sym.data");
52   const std::string rMatrixDataSkewSym("matrix3x3SkewSym.data");
53   const std::string rMatrixDataNoSym("matrix3x3NoSym.data");
54   const std::string rMatrixDataSymPosDef("matrix3x3SymPosDef.data");
55 
56   const std::string cMatrixDataSym("cmatrix3x3Sym.data");
57   const std::string cMatrixDataNoSym("cmatrix3x3NoSym.data");
58   const std::string cMatrixDataSymSelfAjoint("cmatrix3x3SymSelfAjoint.data");
59   const std::string cMatrixDataSymSkewAjoint("cmatrix3x3SymSkewAjoint.data");
60   const std::string cMatrixDataSymPosDef("cmatrix3x3SymPosDef.data");
61 
62   const Number krylovDim = 3;
63   const Real epsilon = 1.e-10;
64   const Number noIteration = 100;
65   const Number vb = 0;
66   GmresSolver gmresRow(krylovDim, epsilon, noIteration, vb);
67   GmresSolver gmresCol(krylovDim, epsilon, noIteration, vb);
68   GmresSolver gmresDual(krylovDim, epsilon, noIteration, vb);
69   GmresSolver gmresSym(krylovDim, epsilon, noIteration, vb);
70 
71   GmresSolver gmresRowCplx(krylovDim, epsilon, noIteration, vb);
72   GmresSolver gmresColCplx(krylovDim, epsilon, noIteration, vb);
73   GmresSolver gmresDualCplx(krylovDim, epsilon, noIteration, vb);
74   GmresSolver gmresSymCplx(krylovDim, epsilon, noIteration, vb);
75 
76   //////////////////////////////////////////////////////////////////////////
77   //
78   //  Vectors
79   //
80   /////////////////////////////////////////////////////////////////////////
81  // Real Vector B
82   Vector<Real> rvB(rowNum, 1.);
83   //out << "Real vector B is: " << rvB <<  std::endl;
84 
85   // Real initial value
86   Vector<Real> rvX0(rowNum, 0.);
87   out << "Real vector x0 is: " << rvX0 << std::endl;
88 
89   // Vector real unknown
90   Vector<Real> rvX(rowNum);
91   Vector<Real> rvXe(rowNum);
92   for(number_t i=0;i<rowNum; i++) rvXe(i+1)=Real(i+1);
93   out << "Real vector Xe is: " << rvXe <<  std::endl;
94 
95   //-------------------------
96   // Vector complex B
97   Vector<Complex> cvB(rowNum, 1.);
98   //out << "The Complex vector B is: " << cvB <<  std::endl;
99 
100   //  Vector initial value
101   Vector<Complex> cvX0(rowNum, 0.);
102   out << "The Complex vector x0 is: " << cvX0 << std::endl;
103 
104   // Vector complex unknown
105   Vector<Complex> cvX(rowNum);
106   Vector<Complex> cvXe(rowNum);
107   for(number_t i=0;i<rowNum; i++) cvXe(i+1)=Complex(1.,1.)*Real(i+1);
108    out << "Complex vector Xe is: " << cvXe <<  std::endl;
109 
110   //////////////////////////////////////////////////////////////////////////
111   //
112   //  Real Large Matrix (Dense Type)
113   //
114   /////////////////////////////////////////////////////////////////////////
115   // No symmetric
116   LargeMatrix<Real> rMatDenseRow(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _dense, _row);
117   LargeMatrix<Real> rMatDenseCol(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _dense, _col);
118   LargeMatrix<Real> rMatDenseDual(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _dense, _dual);
119   LargeMatrix<Real> rMatDenseSym(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _dense, _sym);
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   // Skew Symmetric
128   LargeMatrix<Real> rMatDenseRowSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _dense, _row);
129   LargeMatrix<Real> rMatDenseColSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _dense, _col);
130   LargeMatrix<Real> rMatDenseDualSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _dense, _dual);
131   LargeMatrix<Real> rMatDenseSymSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, _dense, _skewSymmetric);
132 
133   //////////////////////////////////////////////////////////////////////////
134   //
135   //  Complex Large Matrix (Dense Type)
136   //
137   /////////////////////////////////////////////////////////////////////////
138   // No symmetric
139   LargeMatrix<Complex> cMatDenseRow(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _dense, _row);
140   LargeMatrix<Complex> cMatDenseCol(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _dense, _col);
141   LargeMatrix<Complex> cMatDenseDual(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _dense, _dual);
142   LargeMatrix<Complex> cMatDenseSym(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _dense, _sym);
143 
144   // Symmetric
145   LargeMatrix<Complex> cMatDenseRowSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _dense, _row);
146   LargeMatrix<Complex> cMatDenseColSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _dense, _col);
147   LargeMatrix<Complex> cMatDenseDualSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _dense, _dual);
148   LargeMatrix<Complex> cMatDenseSymSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, _dense, _symmetric);
149 
150   // Self Ajoint Symmetric
151   LargeMatrix<Complex> cMatDenseRowSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _dense, _row);
152   LargeMatrix<Complex> cMatDenseColSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _dense, _col);
153   LargeMatrix<Complex> cMatDenseDualSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _dense, _dual);
154   LargeMatrix<Complex> cMatDenseSymSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, _dense, _selfAdjoint);
155 
156   // Skew Ajoint Symmetric
157   LargeMatrix<Complex> cMatDenseRowSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _dense, _row);
158   LargeMatrix<Complex> cMatDenseColSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _dense, _col);
159   LargeMatrix<Complex> cMatDenseDualSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _dense, _dual);
160   LargeMatrix<Complex> cMatDenseSymSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, _dense, _skewAdjoint);
161 
162   //////////////////////////////////////////////////////////////////////////
163   //
164   //  Real Large Matrix (CS Type)
165   //
166   /////////////////////////////////////////////////////////////////////////
167   // No symmetric
168   LargeMatrix<Real> rMatCsRow(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _cs, _row);
169   LargeMatrix<Real> rMatCsCol(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _cs, _col);
170   LargeMatrix<Real> rMatCsDual(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _cs, _dual);
171   LargeMatrix<Real> rMatCsSym(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _cs, _sym);
172 
173   // Symmetric
174   LargeMatrix<Real> rMatCsRowSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _cs, _row);
175   LargeMatrix<Real> rMatCsColSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _cs, _col);
176   LargeMatrix<Real> rMatCsDualSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _cs, _dual);
177   LargeMatrix<Real> rMatCsSymSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, _cs, _symmetric);
178 
179   // Skew Symmetric
180   LargeMatrix<Real> rMatCsRowSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _cs, _row);
181   LargeMatrix<Real> rMatCsColSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _cs, _col);
182   LargeMatrix<Real> rMatCsDualSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _cs, _dual);
183   LargeMatrix<Real> rMatCsSymSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, _cs, _skewSymmetric);
184 
185   //////////////////////////////////////////////////////////////////////////
186   //
187   //  Complex Large Matrix (Cs Type)
188   //
189   /////////////////////////////////////////////////////////////////////////
190   // No symmetric
191   LargeMatrix<Complex> cMatCsRow(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _cs, _row);
192   LargeMatrix<Complex> cMatCsCol(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _cs, _col);
193   LargeMatrix<Complex> cMatCsDual(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _cs, _dual);
194   LargeMatrix<Complex> cMatCsSym(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _cs, _sym);
195 
196   // Symmetric
197   LargeMatrix<Complex> cMatCsRowSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _cs, _row);
198   LargeMatrix<Complex> cMatCsColSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _cs, _col);
199   LargeMatrix<Complex> cMatCsDualSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _cs, _dual);
200   LargeMatrix<Complex> cMatCsSymSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, _cs, _symmetric);
201 
202   // Self Ajoint Symmetric
203   LargeMatrix<Complex> cMatCsRowSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _cs, _row);
204   LargeMatrix<Complex> cMatCsColSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _cs, _col);
205   LargeMatrix<Complex> cMatCsDualSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _cs, _dual);
206   LargeMatrix<Complex> cMatCsSymSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, _cs, _selfAdjoint);
207 
208   // Skew Ajoint Symmetric
209   LargeMatrix<Complex> cMatCsRowSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _cs, _row);
210   LargeMatrix<Complex> cMatCsColSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _cs, _col);
211   LargeMatrix<Complex> cMatCsDualSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _cs, _dual);
212   LargeMatrix<Complex> cMatCsSymSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, _cs, _skewAdjoint);
213 
214   //////////////////////////////////////////////////////////////////////////
215   //
216   //  Real Large Matrix (Skyline Type)
217   //
218   /////////////////////////////////////////////////////////////////////////
219   // No symmetric
220   LargeMatrix<Real> rMatSkylineDual(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _dual);
221   LargeMatrix<Real> rMatSkylineSym(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _sym);
222 
223   // Symmetric
224   LargeMatrix<Real> rMatSkylineDualSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, colNum, _skyline, _dual);
225   LargeMatrix<Real> rMatSkylineSymSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, _skyline, _symmetric);
226 
227   // Skew Symmetric
228   LargeMatrix<Real> rMatSkylineDualSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, colNum, _skyline, _dual);
229   LargeMatrix<Real> rMatSkylineSymSkewSym(inputsPathTo(rMatrixDataSkewSym), _dense, rowNum, _skyline, _skewSymmetric);
230 
231   //////////////////////////////////////////////////////////////////////////
232   //
233   //  Complex Large Matrix (Skyline Type)
234   //
235   /////////////////////////////////////////////////////////////////////////
236   // No symmetric
237   LargeMatrix<Complex> cMatSkylineDual(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _dual);
238   LargeMatrix<Complex> cMatSkylineSym(inputsPathTo(cMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _sym);
239 
240   // Symmetric
241   LargeMatrix<Complex> cMatSkylineDualSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, colNum, _skyline, _dual);
242   LargeMatrix<Complex> cMatSkylineSymSym(inputsPathTo(cMatrixDataSym), _dense, rowNum, _skyline, _symmetric);
243 
244   // Self Ajoint Symmetric
245   LargeMatrix<Complex> cMatSkylineDualSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, colNum, _skyline, _dual);
246   LargeMatrix<Complex> cMatSkylineSymSymSelfAjoint(inputsPathTo(cMatrixDataSymSelfAjoint), _dense, rowNum, _skyline, _selfAdjoint);
247 
248   // Skew Ajoint Symmetric
249   LargeMatrix<Complex> cMatSkylineDualSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint), _dense, rowNum, colNum, _skyline, _dual);
250   LargeMatrix<Complex> cMatSkylineSymSymSkewAjoint(inputsPathTo(cMatrixDataSymSkewAjoint),	 _dense, rowNum, _skyline, _skewAdjoint);
251 
252   //////////////////////////////////////////////////////////////////////////
253   //
254   //  Real Large Matrix to factorize (Skyline Type) as a preconditioner
255   //
256   /////////////////////////////////////////////////////////////////////////
257 //  LargeMatrix<Real> rMat2Factorize(inputsPathTo(rMatrixDataSymPosDef), _dense, rowNum, _skyline, _symmetric);
258 //  ldltFactorize(rMat2Factorize);
259 //  Preconditioner<LargeMatrix<Real> > pcLdlt(rMat2Factorize, _ldltPrec, 1.5);
260 //
261 //  LargeMatrix<Real> rMat2LuFac(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, _skyline, _noSymmetry);
262 //  luFactorize(rMat2LuFac);
263 //  Preconditioner<LargeMatrix<Real> > pcLu(rMat2LuFac, _ldltPrec, 1.5);
264 //
265 //  LargeMatrix<Real> rMatCsSymSymPc(inputsPathTo(rMatrixDataSym), _dense, rowNum, _cs, _symmetric);
266 //  Preconditioner<LargeMatrix<Real> > pcSorCsSym(rMatCsSymSymPc, _ssorPrec, 1.5);
267 //
268 //  //////////////////////////////////////////////////////////////////////////
269 //  //
270 //  //  Complex Large Matrix to factorize (Skyline Type) as a preconditioner
271 //  //
272 //  /////////////////////////////////////////////////////////////////////////
273 //  LargeMatrix<Complex> cMat2Factorize(inputsPathTo(cMatrixDataSymPosDef), _dense, rowNum, _skyline, _selfAdjoint);
274 //  ldlstarFactorize(cMat2Factorize);
275 //  Preconditioner<LargeMatrix<Complex> > pcLdlStar(cMat2Factorize, _ldlstarPrec, 1.5);
276 
277   ////////////////////////////////////////////////////////////////////////////////
278   //
279   //------------------------------------------------------------------------------
280   // Test with Real value
281   //-----------------------------------------------------------------------------
282   //
283   // I. Solver without a preconditioner
284   // Test with DENSE STORAGE
285   // Test with symmetric matrices
286   rvB = rMatDenseRowSym *rvXe;
287   out << "Real vector B is: " << rvB <<  std::endl;
288   rvX = gmresRow(rMatDenseRowSym, rvB, rvX0, _real);
289   out << "The GmresSolver result with real dense row and sym matrix  is: " << rvX << std::endl;
290   rvX = gmresCol(rMatDenseColSym, rvB, rvX0, _real);
291   out << "The GmresSolver result with real dense col and sym matrix is: " << rvX << std::endl;
292   rvX = gmresDual(rMatDenseDualSym, rvB, rvX0, _real);
293   out << "The GmresSolver result with real dense dual and sym matrix is: " << rvX << std::endl;
294   rvX = gmresSym(rMatDenseSymSym, rvB, rvX0, _real);
295   out << "The GmresSolver result with real dense sym and sym matrix is: " << rvX << std::endl << std::endl;
296 
297   // Test with skew symmetric matrices
298   rvB = rMatDenseRowSkewSym *rvXe;
299   out << "Real vector B is: " << rvB <<  std::endl;
300   rvX = gmresRow(rMatDenseRowSkewSym, rvB, rvX0, _real);
301   out << "The GmresSolver result with real dense row and skew sym matrix is: " << rvX << std::endl;
302   rvX = gmresCol(rMatDenseColSkewSym, rvB, rvX0, _real);
303   out << "The GmresSolver result with real dense col and skew sym matrix is: " << rvX << std::endl;
304   rvX = gmresDual(rMatDenseDualSkewSym, rvB, rvX0, _real);
305   out << "The GmresSolver result with real dense dual and skew sym matrix is: " << rvX << std::endl;
306   rvX = gmresSym(rMatDenseSymSkewSym, rvB, rvX0, _real);
307   out << "The GmresSolver result with real dense sym and skew sym matrix is: " << rvX << std::endl << std::endl;
308 
309   // Test with non-symmetric matrices
310   rvB = rMatDenseRow *rvXe;
311   out << "Real vector B is: " << rvB <<  std::endl;
312   rvX = gmresRow(rMatDenseRow, rvB, rvX0, _real);
313   out << "The GmresSolver result with real dense row and non sym matrix is: " << rvX << std::endl;
314   rvX = gmresCol(rMatDenseCol, rvB, rvX0, _real);
315   out << "The GmresSolver result with real dense col and non sym matrix is: " << rvX << std::endl;
316   rvX = gmresDual(rMatDenseDual, rvB, rvX0, _real);
317   out << "The GmresSolver result with real dense dual and non sym matrix is: " << rvX << std::endl;
318   rvX = gmresSym(rMatDenseSym, rvB, rvX0, _real);
319   out << "The GmresSolver result with real dense sym and non sym matrix is: " << rvX << std::endl << std::endl;
320 
321   //---------------------------------
322   // Test with CS STORAGE
323   // Test with symmetric matrices
324   rvB = rMatCsRowSym *rvXe;
325   out << "Real vector B is: " << rvB <<  std::endl;
326   rvX = gmresRow(rMatCsRowSym, rvB, rvX0, _real);
327   out << "The GmresSolver result with real Cs row and sym matrix is: " << rvX << std::endl;
328   rvX = gmresCol(rMatCsColSym, rvB, rvX0, _real);
329   out << "The GmresSolver result with real Cs col and sym matrix is: " << rvX << std::endl;
330   rvX = gmresDual(rMatCsDualSym, rvB, rvX0, _real);
331   out << "The GmresSolver result with real Cs dual and sym matrix is: " << rvX << std::endl;
332   rvX = gmresSym(rMatCsSymSym, rvB, rvX0, _real);
333   out << "The GmresSolver result with real Cs sym and sym matrix is: " << rvX << std::endl << std::endl;
334 
335   // Test with skew symmetric matrices
336   rvB = rMatCsRowSkewSym *rvXe;
337   out << "Real vector B is: " << rvB <<  std::endl;
338   rvX = gmresRow(rMatCsRowSkewSym, rvB, rvX0, _real);
339   out << "The GmresSolver result with real Cs row and skew sym matrix is: " << rvX << std::endl;
340   rvX = gmresCol(rMatCsColSkewSym, rvB, rvX0, _real);
341   out << "The GmresSolver result with real Cs col and skew sym matrix is: " << rvX << std::endl;
342   rvX = gmresDual(rMatCsDualSkewSym, rvB, rvX0, _real);
343   out << "The GmresSolver result with real Cs dual and skew sym matrix is: " << rvX << std::endl;
344   rvX = gmresSym(rMatCsSymSkewSym, rvB, rvX0, _real);
345   out << "The GmresSolver result with real Cs sym and skew sym matrix is: " << rvX << std::endl << std::endl;
346 
347   // Test with non symmetric matrices
348   rvB = rMatCsRow *rvXe;
349   out << "Real vector B is: " << rvB <<  std::endl;
350   rvX = gmresRow(rMatCsRow, rvB, rvX0, _real);
351   out << "The GmresSolver result with real Cs row and non sym matrix is: " << rvX << std::endl;
352   rvX = gmresCol(rMatCsCol, rvB, rvX0, _real);
353   out << "The GmresSolver result with real Cs col and non sym matrix is: " << rvX << std::endl;
354   rvX = gmresDual(rMatCsDual, rvB, rvX0, _real);
355   out << "The GmresSolver result with real Cs dual and non sym matrix is: " << rvX << std::endl;
356   rvX = gmresSym(rMatCsSym, rvB, rvX0, _real);
357   out << "The GmresSolver result with real Cs sym and non sym matrix is: " << rvX << std::endl << std::endl;
358 
359   //------------------------------
360   // Test with SKYLINE STORAGE
361   // Test with symmetric matrices
362   rvB = rMatSkylineDualSym *rvXe;
363   out << "Real vector B is: " << rvB <<  std::endl;
364   rvX = gmresDual(rMatSkylineDualSym, rvB, rvX0, _real);
365   out << "The GmresSolver result with real Skyline dual and sym matrix is: " << rvX << std::endl;
366   rvX = gmresSym(rMatSkylineSymSym, rvB, rvX0, _real);
367   out << "The GmresSolver result with real Skyline sym and sym matrix is: " << rvX << std::endl << std::endl;
368 
369   // Test with skew symmetric matrices
370   rvB = rMatSkylineDualSkewSym *rvXe;
371   out << "Real vector B is: " << rvB <<  std::endl;
372   rvX = gmresDual(rMatSkylineDualSkewSym, rvB, rvX0, _real);
373   out << "The GmresSolver result with real Skyline dual and skew sym matrix is: " << rvX << std::endl;
374   rvX = gmresSym(rMatSkylineSymSkewSym, rvB, rvX0, _real);
375   out << "The GmresSolver result with real Skyline sym and skew sym matrix is: " << rvX << std::endl << std::endl;
376 
377   // Test with non symmetric matrices
378   rvB = rMatSkylineDual *rvXe;
379   out << "Real vector B is: " << rvB <<  std::endl;
380   rvX = gmresDual(rMatSkylineDual, rvB, rvX0, _real);
381   out << "The GmresSolver result with real Skyline dual and skew non matrix is: " << rvX << std::endl;
382   rvX = gmresSym(rMatSkylineSym, rvB, rvX0, _real);
383   out << "The GmresSolver result with real Skyline non and skew sym matrix is: " << rvX << std::endl << std::endl;
384 
385   //
386   // -----------------------------------------------------------------------------
387   // Test with Complex values (Complex x Complex = Complex)
388   //-----------------------------------------------------------------------------
389   //
390   // Test with DENSE STORAGE
391   // Test with symmetric matrices
392   cvB = cMatDenseRowSym *cvXe;
393   out << "Complex vector B is: " << cvB <<  std::endl;
394   cvX = gmresRowCplx(cMatDenseRowSym, cvB, cvX0, xlifepp::_complex);
395   out << "The GmresSolver result with complex dense row and sym matrix is: " << cvX << std::endl;
396   cvX = gmresColCplx(cMatDenseColSym, cvB, cvX0, xlifepp::_complex);
397   out << "The GmresSolver result with complex dense col and sym matrix  is: " << cvX << std::endl;
398   cvX = gmresDualCplx(cMatDenseDualSym, cvB, cvX0, xlifepp::_complex);
399   out << "The GmresSolver result with complex dense dual and sym matrix is: " << cvX << std::endl;
400   cvX = gmresSymCplx(cMatDenseSymSym, cvB, cvX0, xlifepp::_complex);
401   out << "The GmresSolver result with complex dense sym and sym matrix is: " << cvX << std::endl << std::endl;
402 
403   // Test with Sym SelfAjoint symmetric matrices
404   cvB = cMatDenseRowSymSelfAjoint *cvXe;
405   out << "Complex vector B is: " << cvB <<  std::endl;
406   cvX = gmresRowCplx(cMatDenseRowSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
407   out << "The GmresSolver result with complex dense row and selfajoint sym matrix is: " << cvX << std::endl;
408   cvX = gmresColCplx(cMatDenseColSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
409   out << "The GmresSolver result with complex dense col and selfajoint sym matrix is: " << cvX << std::endl;
410   cvX = gmresDualCplx(cMatDenseDualSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
411   out << "The GmresSolver result with complex dense dual and selfajoint sym matrix is: " << cvX << std::endl;
412   cvX = gmresSymCplx(cMatDenseSymSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
413   out << "The GmresSolver result with complex dense sym and selfajoint sym matrix is: " << cvX << std::endl << std::endl;
414 
415   // Test with Skew Ajoint symmetric matrices
416   cvB = cMatDenseRowSymSkewAjoint *cvXe;
417   out << "Complex vector B is: " << cvB <<  std::endl;
418   cvX = gmresRowCplx(cMatDenseRowSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
419   out << "The GmresSolver result with complex dense row and skewajoint sym matrix is: " << cvX << std::endl;
420   cvX = gmresColCplx(cMatDenseColSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
421   out << "The GmresSolver result with complex dense col and skewajoint sym matrix is: " << cvX << std::endl;
422   cvX = gmresDualCplx(cMatDenseDualSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
423   out << "The GmresSolver result with complex dense dual and skewajoint sym matrix is: " << cvX << std::endl;
424   cvX = gmresSymCplx(cMatDenseSymSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
425   out << "The GmresSolver result with complex dense sym and skewajoint sym matrix is: " << cvX << std::endl << std::endl;
426 
427   // Test with no symmetric matrices
428   cvB = cMatDenseRow *cvXe;
429   out << "Complex vector B is: " << cvB <<  std::endl;
430   cvX = gmresRowCplx(cMatDenseRow, cvB, cvX0, xlifepp::_complex);
431   out << "The GmresSolver result with complex dense row and non-sym matrix is: " << cvX << std::endl;
432   cvX = gmresColCplx(cMatDenseCol, cvB, cvX0, xlifepp::_complex);
433   out << "The GmresSolver result with complex dense col and non-sym matrix is: " << cvX << std::endl;
434   cvX = gmresDualCplx(cMatDenseDual, cvB, cvX0, xlifepp::_complex);
435   out << "The GmresSolver result with complex dense dual and non-sym matrix is: " << cvX << std::endl;
436   cvX = gmresSymCplx(cMatDenseSym, cvB, cvX0, xlifepp::_complex);
437   out << "The GmresSolver result with complex dense sym and non-sym matrix is: " << cvX << std::endl << std::endl;
438 
439   //---------------------------------
440   // Test with CS STORAGE
441   // Test with symmetric matrices
442   cvB = cMatCsRowSym *cvXe;
443   out << "Complex vector B is: " << cvB <<  std::endl;
444   cvX = gmresRowCplx(cMatCsRowSym, cvB, cvX0, xlifepp::_complex);
445   out << "The GmresSolver result with complex Cs row and sym matrix is: " << cvX << std::endl;
446   cvX = gmresColCplx(cMatCsColSym, cvB, cvX0, xlifepp::_complex);
447   out << "The GmresSolver result with complex Cs col and sym matrix is: " << cvX << std::endl;
448   cvX = gmresDualCplx(cMatCsDualSym, cvB, cvX0, xlifepp::_complex);
449   out << "The GmresSolver result with complex Cs dual and sym matrix is: " << cvX << std::endl;
450   cvX = gmresSymCplx(cMatCsSymSym, cvB, cvX0, xlifepp::_complex);
451   out << "The GmresSolver result with complex Cs sym and sym matrix is: " << cvX << std::endl << std::endl;
452 
453   // Test with selfAjoint symmetric matrices
454   cvB = cMatCsRowSymSelfAjoint *cvXe;
455   out << "Complex vector B is: " << cvB <<  std::endl;
456   cvX = gmresRowCplx(cMatCsRowSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
457   out << "The GmresSolver result with complex Cs row and selfajoint sym matrix is: " << cvX << std::endl;
458   cvX = gmresColCplx(cMatCsColSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
459   out << "The GmresSolver result with complex Cs col and selfajoint sym matrix is: " << cvX << std::endl;
460   cvX = gmresDualCplx(cMatCsDualSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
461   out << "The GmresSolver result with complex Cs dual and selfajoint sym matrix is: " << cvX << std::endl;
462   cvX = gmresSymCplx(cMatCsSymSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
463   out << "The GmresSolver result with complex Cs sym and selfajoint sym matrix is: " << cvX << std::endl << std::endl;
464 
465   // Test with skewAjoint symmetric matrices
466   cvB = cMatCsRowSymSkewAjoint*cvXe;
467   out << "Complex vector B is: " << cvB <<  std::endl;
468   cvX = gmresRowCplx(cMatCsRowSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
469   out << "The GmresSolver result with complex Cs row and skewajoint sym matrix is: " << cvX << std::endl;
470   cvX = gmresColCplx(cMatCsColSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
471   out << "The GmresSolver result with complex Cs col and skewajoint sym matrix is: " << cvX << std::endl;
472   cvX = gmresDualCplx(cMatCsDualSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
473   out << "The GmresSolver result with complex Cs dual and skewajoint sym matrix is: " << cvX << std::endl;
474   cvX = gmresSymCplx(cMatCsSymSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
475   out << "The GmresSolver result with complex Cs sym and skewajoint sym matrix is: " << cvX << std::endl << std::endl;
476 
477   // Test with no symmetric matrices
478   cvB = cMatCsRow*cvXe;
479   out << "Complex vector B is: " << cvB <<  std::endl;
480   cvX = gmresRowCplx(cMatCsRow, cvB, cvX0, xlifepp::_complex);
481   out << "The GmresSolver result with complex Cs row and non-sym matrix is: " << cvX << std::endl;
482   cvX = gmresColCplx(cMatCsCol, cvB, cvX0, xlifepp::_complex);
483   out << "The GmresSolver result with complex Cs col and non-sym matrix is: " << cvX << std::endl;
484   cvX = gmresDualCplx(cMatCsDual, cvB, cvX0, xlifepp::_complex);
485   out << "The GmresSolver result with complex Cs dual and non-sym matrix is: " << cvX << std::endl;
486   cvX = gmresSymCplx(cMatCsSym, cvB, cvX0, xlifepp::_complex);
487   out << "The GmresSolver result with complex Cs sym and non-sym matrix is: " << cvX << std::endl << std::endl;
488 
489   //---------------------------------
490   // Test with SKYLINE STORAGE
491   // Test with symmetric matrices
492   cvB = cMatSkylineDualSym*cvXe;
493   out << "Complex vector B is: " << cvB <<  std::endl;
494   cvX = gmresDualCplx(cMatSkylineDualSym, cvB, cvX0, xlifepp::_complex);
495   out << "The GmresSolver result with complex Skyline dual and sym matrix is: " << cvX << std::endl;
496   cvX = gmresSymCplx(cMatSkylineSymSym, cvB, cvX0, xlifepp::_complex);
497   out << "The GmresSolver result with complex Skyline sym and sym matrix is: " << cvX << std::endl << std::endl;
498 
499   // Test with self ajoint sym matrices
500   cvB = cMatSkylineDualSymSelfAjoint*cvXe;
501   out << "Complex vector B is: " << cvB <<  std::endl;
502   cvX = gmresDualCplx(cMatSkylineDualSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
503   out << "The GmresSolver result with complex Skyline dual and selfajoint sym matrix is: " << cvX << std::endl;
504   cvX = gmresSymCplx(cMatSkylineSymSymSelfAjoint, cvB, cvX0, xlifepp::_complex);
505   out << "The GmresSolver result with complex Skyline sym and selfajoint sym matrix is: " << cvX << std::endl << std::endl;
506 
507   // Test with skew ajoint matrices
508   cvB = cMatSkylineDualSymSkewAjoint*cvXe;
509   out << "Complex vector B is: " << cvB <<  std::endl;
510   cvX = gmresDualCplx(cMatSkylineDualSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
511   out << "The GmresSolver result with complex Skyline dual and skewajoint sym matrix is: " << cvX << std::endl;
512   cvX = gmresSymCplx(cMatSkylineSymSymSkewAjoint, cvB, cvX0, xlifepp::_complex);
513   out << "The GmresSolver result with complex Skyline sym and skewajoint sym matrix is: " << cvX << std::endl << std::endl;
514 
515   // Test with no symmetric matrices
516   cvB = cMatSkylineDual*cvXe;
517   out << "Complex vector B is: " << cvB <<  std::endl;
518   cvX = gmresDualCplx(cMatSkylineDual, cvB, cvX0, xlifepp::_complex);
519   out << "The GmresSolver result with complex Skyline dual and non-sym matrix is: " << cvX << std::endl;
520   cvX = gmresSymCplx(cMatSkylineSym, cvB, cvX0, xlifepp::_complex);
521   out << "The GmresSolver result with complex Skyline sym and non-sym matrix is: " << cvX << std::endl << std::endl;
522 
523 /*
524   // II. Solver with a preconditioner
525   // Test with DENSE STORAGE
526   // Test with symmetric matrices
527   gmresRow.resetSolver();
528   gmresCol.resetSolver();
529   gmresDual.resetSolver();
530   gmresSym.resetSolver();
531   rvX = gmresRow(rMatDenseRowSym, rvB, rvX0, pcLdlt, _real);
532   out << "The GmresSolver result with real dense row sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
533   rvX = gmresCol(rMatDenseColSym, rvB, rvX0, pcLdlt, _real);
534   out << "The GmresSolver result with real dense col sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
535   rvX = gmresDual(rMatDenseDualSym, rvB, rvX0, pcLdlt, _real);
536   out << "The GmresSolver result with real dense dual sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
537   rvX = gmresSym(rMatDenseSymSym, rvB, rvX0, pcLdlt, _real);
538   out << "The GmresSolver result with real dense sym sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
539 
540   gmresRow.resetSolver();
541   gmresCol.resetSolver();
542   gmresDual.resetSolver();
543   gmresSym.resetSolver();
544   rvX = gmresRow(rMatDenseRowSym, rvB, rvX0, pcLu, _real);
545   out << "The GmresSolver result with real dense row sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
546   rvX = gmresCol(rMatDenseColSym, rvB, rvX0, pcLu, _real);
547   out << "The GmresSolver result with real dense col sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
548   rvX = gmresDual(rMatDenseDualSym, rvB, rvX0, pcLu, _real);
549   out << "The GmresSolver result with real dense dual sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
550   rvX = gmresSym(rMatDenseSymSym, rvB, rvX0, pcLu, _real);
551   out << "The GmresSolver result with real dense sym sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl << std::endl;
552 
553   gmresRow.resetSolver();
554   gmresCol.resetSolver();
555   gmresDual.resetSolver();
556   gmresSym.resetSolver();
557   rvX = gmresRow(rMatDenseRowSym, rvB, rvX0, pcSorCsSym, _real);
558   out << "The GmresSolver result with real dense row sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
559   rvX = gmresCol(rMatDenseColSym, rvB, rvX0, pcSorCsSym, _real);
560   out << "The GmresSolver result with real dense col sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
561   rvX = gmresDual(rMatDenseDualSym, rvB, rvX0, pcSorCsSym, _real);
562   out << "The GmresSolver result with real dense dual sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
563   rvX = gmresSym(rMatDenseSymSym, rvB, rvX0, pcSorCsSym, _real);
564   out << "The GmresSolver result with real dense sym sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl << std::endl;
565 
566   // Test with skew symmetric matrices
567   gmresRow.resetSolver();
568   gmresCol.resetSolver();
569   gmresDual.resetSolver();
570   gmresSym.resetSolver();
571   rvX = gmresRow(rMatDenseRowSkewSym, rvB, rvX0, pcLdlt, _real);
572   out << "The GmresSolver result with real dense row skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
573   rvX = gmresCol(rMatDenseColSkewSym, rvB, rvX0, pcLdlt, _real);
574   out << "The GmresSolver result with real dense col skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
575   rvX = gmresDual(rMatDenseDualSkewSym, rvB, rvX0, pcLdlt, _real);
576   out << "The GmresSolver result with real dense dual skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
577   rvX = gmresSym(rMatDenseSymSkewSym, rvB, rvX0, pcLdlt, _real);
578   out << "The GmresSolver result with real dense sym skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
579 
580   gmresRow.resetSolver();
581   gmresCol.resetSolver();
582   gmresDual.resetSolver();
583   gmresSym.resetSolver();
584   rvX = gmresRow(rMatDenseRowSkewSym, rvB, rvX0, pcLu, _real);
585   out << "The GmresSolver result with real dense row skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
586   rvX = gmresCol(rMatDenseColSkewSym, rvB, rvX0, pcLu, _real);
587   out << "The GmresSolver result with real dense col skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
588   rvX = gmresDual(rMatDenseDualSkewSym, rvB, rvX0, pcLu, _real);
589   out << "The GmresSolver result with real dense dual skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
590   rvX = gmresSym(rMatDenseSymSkewSym, rvB, rvX0, pcLu, _real);
591   out << "The GmresSolver result with real dense sym skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl << std::endl;
592 
593   gmresRow.resetSolver();
594   gmresCol.resetSolver();
595   gmresDual.resetSolver();
596   gmresSym.resetSolver();
597   rvX = gmresRow(rMatDenseRowSkewSym, rvB, rvX0, pcSorCsSym, _real);
598   out << "The GmresSolver result with real dense row skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
599   rvX = gmresCol(rMatDenseColSkewSym, rvB, rvX0, pcSorCsSym, _real);
600   out << "The GmresSolver result with real dense col skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
601   rvX = gmresDual(rMatDenseDualSkewSym, rvB, rvX0, pcSorCsSym, _real);
602   out << "The GmresSolver result with real dense dual skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
603   rvX = gmresSym(rMatDenseSymSkewSym, rvB, rvX0, pcSorCsSym, _real);
604   out << "The GmresSolver result with real dense sym skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl << std::endl;
605 
606   // Test with non symmetric matrices
607   gmresRow.resetSolver();
608   gmresCol.resetSolver();
609   gmresDual.resetSolver();
610   gmresSym.resetSolver();
611   rvX = gmresRow(rMatDenseRow, rvB, rvX0, pcLdlt, _real);
612   out << "The GmresSolver result with real dense row non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
613   rvX = gmresCol(rMatDenseCol, rvB, rvX0, pcLdlt, _real);
614   out << "The GmresSolver result with real dense col non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
615   rvX = gmresDual(rMatDenseDual, rvB, rvX0, pcLdlt, _real);
616   out << "The GmresSolver result with real dense dual non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
617   rvX = gmresSym(rMatDenseSym, rvB, rvX0, pcLdlt, _real);
618   out << "The GmresSolver result with real dense sym non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
619 
620   gmresRow.resetSolver();
621   gmresCol.resetSolver();
622   gmresDual.resetSolver();
623   gmresSym.resetSolver();
624   rvX = gmresRow(rMatDenseRow, rvB, rvX0, pcLu, _real);
625   out << "The GmresSolver result with real dense row non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
626   rvX = gmresCol(rMatDenseCol, rvB, rvX0, pcLu, _real);
627   out << "The GmresSolver result with real dense col non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
628   rvX = gmresDual(rMatDenseDual, rvB, rvX0, pcLu, _real);
629   out << "The GmresSolver result with real dense dual non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
630   rvX = gmresSym(rMatDenseSym, rvB, rvX0, pcLu, _real);
631   out << "The GmresSolver result with real dense sym non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl << std::endl;
632 
633   gmresRow.resetSolver();
634   gmresCol.resetSolver();
635   gmresDual.resetSolver();
636   gmresSym.resetSolver();
637   rvX = gmresRow(rMatDenseRow, rvB, rvX0, pcSorCsSym, _real);
638   out << "The GmresSolver result with real dense row non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
639   rvX = gmresCol(rMatDenseCol, rvB, rvX0, pcSorCsSym, _real);
640   out << "The GmresSolver result with real dense col non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
641   rvX = gmresDual(rMatDenseDual, rvB, rvX0, pcSorCsSym, _real);
642   out << "The GmresSolver result with real dense dual non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
643   rvX = gmresSym(rMatDenseSym, rvB, rvX0, pcSorCsSym, _real);
644   out << "The GmresSolver result with real dense sym non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl << std::endl;
645 
646   //---------------------------------
647   // Test with CS STORAGE
648   // Test with symmetric matrices
649   gmresRow.resetSolver();
650   gmresCol.resetSolver();
651   gmresDual.resetSolver();
652   gmresSym.resetSolver();
653   rvX = gmresRow(rMatCsRowSym, rvB, rvX0, pcLdlt, _real);
654   out << "The GmresSolver result with real Cs row sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
655   rvX = gmresCol(rMatCsColSym, rvB, rvX0, pcLdlt, _real);
656   out << "The GmresSolver result with real Cs col sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
657   rvX = gmresDual(rMatCsDualSym, rvB, rvX0, pcLdlt, _real);
658   out << "The GmresSolver result with real Cs dual sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
659   rvX = gmresSym(rMatCsSymSym, rvB, rvX0, pcLdlt, _real);
660   out << "The GmresSolver result with real Cs sym sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
661 
662   gmresRow.resetSolver();
663   gmresCol.resetSolver();
664   gmresDual.resetSolver();
665   gmresSym.resetSolver();
666   rvX = gmresRow(rMatCsRowSym, rvB, rvX0, pcLu, _real);
667   out << "The GmresSolver result with real Cs row sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
668   rvX = gmresCol(rMatCsColSym, rvB, rvX0, pcLu, _real);
669   out << "The GmresSolver result with real Cs col sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
670   rvX = gmresDual(rMatCsDualSym, rvB, rvX0, pcLu, _real);
671   out << "The GmresSolver result with real Cs dual sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
672   rvX = gmresSym(rMatCsSymSym, rvB, rvX0, pcLu, _real);
673   out << "The GmresSolver result with real Cs sym sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl << std::endl;
674 
675   gmresRow.resetSolver();
676   gmresCol.resetSolver();
677   gmresDual.resetSolver();
678   gmresSym.resetSolver();
679   rvX = gmresRow(rMatCsRowSym, rvB, rvX0, pcSorCsSym, _real);
680   out << "The GmresSolver result with real Cs row sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
681   rvX = gmresCol(rMatCsColSym, rvB, rvX0, pcSorCsSym, _real);
682   out << "The GmresSolver result with real Cs col sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
683   rvX = gmresDual(rMatCsDualSym, rvB, rvX0, pcSorCsSym, _real);
684   out << "The GmresSolver result with real Cs dual sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
685   rvX = gmresSym(rMatCsSymSym, rvB, rvX0, pcSorCsSym, _real);
686   out << "The GmresSolver result with real Cs sym sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl << std::endl;
687 
688   // Test with skew symmetric matrices
689   gmresRow.resetSolver();
690   gmresCol.resetSolver();
691   gmresDual.resetSolver();
692   gmresSym.resetSolver();
693   rvX = gmresRow(rMatCsRowSkewSym, rvB, rvX0, pcLdlt, _real);
694   out << "The GmresSolver result with real Cs row skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
695   rvX = gmresCol(rMatCsColSkewSym, rvB, rvX0, pcLdlt, _real);
696   out << "The GmresSolver result with real Cs col skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
697   rvX = gmresDual(rMatCsDualSkewSym, rvB, rvX0, pcLdlt, _real);
698   out << "The GmresSolver result with real Cs dual skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
699   rvX = gmresSym(rMatCsSymSkewSym, rvB, rvX0, pcLdlt, _real);
700   out << "The GmresSolver result with real Cs sym skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
701 
702   gmresRow.resetSolver();
703   gmresCol.resetSolver();
704   gmresDual.resetSolver();
705   gmresSym.resetSolver();
706   rvX = gmresRow(rMatCsRowSkewSym, rvB, rvX0, pcLu, _real);
707   out << "The GmresSolver result with real Cs row skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
708   rvX = gmresCol(rMatCsColSkewSym, rvB, rvX0, pcLu, _real);
709   out << "The GmresSolver result with real Cs col skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
710   rvX = gmresDual(rMatCsDualSkewSym, rvB, rvX0, pcLu, _real);
711   out << "The GmresSolver result with real Cs dual skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
712   rvX = gmresSym(rMatCsSymSkewSym, rvB, rvX0, pcLu, _real);
713   out << "The GmresSolver result with real Cs sym skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl << std::endl;
714 
715   gmresRow.resetSolver();
716   gmresCol.resetSolver();
717   gmresDual.resetSolver();
718   gmresSym.resetSolver();
719   rvX = gmresRow(rMatCsRowSkewSym, rvB, rvX0, pcSorCsSym, _real);
720   out << "The GmresSolver result with real Cs row skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
721   rvX = gmresCol(rMatCsColSkewSym, rvB, rvX0, pcSorCsSym, _real);
722   out << "The GmresSolver result with real Cs col skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
723   rvX = gmresDual(rMatCsDualSkewSym, rvB, rvX0, pcSorCsSym, _real);
724   out << "The GmresSolver result with real Cs dual skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
725   rvX = gmresSym(rMatCsSymSkewSym, rvB, rvX0, pcSorCsSym, _real);
726   out << "The GmresSolver result with real Cs sym skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl << std::endl;
727 
728   // Test with non symmetric matrices
729   gmresRow.resetSolver();
730   gmresCol.resetSolver();
731   gmresDual.resetSolver();
732   gmresSym.resetSolver();
733   rvX = gmresRow(rMatCsRow, rvB, rvX0, pcLdlt, _real);
734   out << "The GmresSolver result with real Cs row non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
735   rvX = gmresCol(rMatCsCol, rvB, rvX0, pcLdlt, _real);
736   out << "The GmresSolver result with real Cs col non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
737   rvX = gmresDual(rMatCsDual, rvB, rvX0, pcLdlt, _real);
738   out << "The GmresSolver result with real Cs dual non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
739   rvX = gmresSym(rMatCsSym, rvB, rvX0, pcLdlt, _real);
740   out << "The GmresSolver result with real Cs sym non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
741 
742   gmresRow.resetSolver();
743   gmresCol.resetSolver();
744   gmresDual.resetSolver();
745   gmresSym.resetSolver();
746   rvX = gmresRow(rMatCsRow, rvB, rvX0, pcLu, _real);
747   out << "The GmresSolver result with real Cs row non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
748   rvX = gmresCol(rMatCsCol, rvB, rvX0, pcLu, _real);
749   out << "The GmresSolver result with real Cs col non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
750   rvX = gmresDual(rMatCsDual, rvB, rvX0, pcLu, _real);
751   out << "The GmresSolver result with real Cs dual non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
752   rvX = gmresSym(rMatCsSym, rvB, rvX0, pcLu, _real);
753   out << "The GmresSolver result with real Cs sym non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl << std::endl;
754 
755   gmresRow.resetSolver();
756   gmresCol.resetSolver();
757   gmresDual.resetSolver();
758   gmresSym.resetSolver();
759   rvX = gmresRow(rMatCsRow, rvB, rvX0, pcSorCsSym, _real);
760   out << "The GmresSolver result with real Cs row non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
761   rvX = gmresCol(rMatCsCol, rvB, rvX0, pcSorCsSym, _real);
762   out << "The GmresSolver result with real Cs col non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
763   rvX = gmresDual(rMatCsDual, rvB, rvX0, pcSorCsSym, _real);
764   out << "The GmresSolver result with real Cs dual non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
765   rvX = gmresSym(rMatCsSym, rvB, rvX0, pcSorCsSym, _real);
766   out << "The GmresSolver result with real Cs sym non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl << std::endl;
767 
768   //------------------------------
769   // Test with SKYLINE STORAGE
770   // Test with symmetric matrices
771   gmresDual.resetSolver();
772   gmresSym.resetSolver();
773   rvX = gmresDual(rMatSkylineDualSym, rvB, rvX0, pcLdlt, _real);
774   out << "The GmresSolver result with real Skyline dual sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
775   rvX = gmresSym(rMatSkylineSymSym, rvB, rvX0, pcLdlt, _real);
776   out << "The GmresSolver result with real Skyline sym sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
777 
778   gmresDual.resetSolver();
779   gmresSym.resetSolver();
780   rvX = gmresDual(rMatSkylineDualSym, rvB, rvX0, pcLu, _real);
781   out << "The GmresSolver result with real Skyline dual sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
782   rvX = gmresSym(rMatSkylineSymSym, rvB, rvX0, pcLu, _real);
783   out << "The GmresSolver result with real Skyline sym sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl << std::endl;
784 
785   gmresDual.resetSolver();
786   gmresSym.resetSolver();
787   rvX = gmresDual(rMatSkylineDualSym, rvB, rvX0, pcSorCsSym, _real);
788   out << "The GmresSolver result with real Skyline dual sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
789   rvX = gmresSym(rMatSkylineSymSym, rvB, rvX0, pcSorCsSym, _real);
790   out << "The GmresSolver result with real Skyline sym sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl << std::endl;
791 
792   // Test with skew symmetric matrices
793   gmresDual.resetSolver();
794   gmresSym.resetSolver();
795   rvX = gmresDual(rMatSkylineDualSkewSym, rvB, rvX0, pcLdlt, _real);
796   out << "The GmresSolver result with real Skyline dual skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
797   rvX = gmresSym(rMatSkylineSymSkewSym, rvB, rvX0, pcLdlt, _real);
798   out << "The GmresSolver result with real Skyline sym skew sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
799 
800   gmresDual.resetSolver();
801   gmresSym.resetSolver();
802   rvX = gmresDual(rMatSkylineDualSkewSym, rvB, rvX0, pcLu, _real);
803   out << "The GmresSolver result with real Skyline dual skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
804   rvX = gmresSym(rMatSkylineSymSkewSym, rvB, rvX0, pcLu, _real);
805   out << "The GmresSolver result with real Skyline sym skew sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl << std::endl;
806 
807   gmresDual.resetSolver();
808   gmresSym.resetSolver();
809   rvX = gmresDual(rMatSkylineDualSkewSym, rvB, rvX0, pcSorCsSym, _real);
810   out << "The GmresSolver result with real Skyline dual skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
811   rvX = gmresSym(rMatSkylineSymSkewSym, rvB, rvX0, pcSorCsSym, _real);
812   out << "The GmresSolver result with real Skyline sym skew sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl << std::endl;
813 
814   // Test with non symmetric matrices
815   gmresDual.resetSolver();
816   gmresSym.resetSolver();
817   rvX = gmresDual(rMatSkylineDual, rvB, rvX0, pcLdlt, _real);
818   out << "The GmresSolver result with real Skyline dual non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl;
819   rvX = gmresSym(rMatSkylineSym, rvB, rvX0, pcLdlt, _real);
820   out << "The GmresSolver result with real Skyline sym non-sym matrix  and a ldlt factorized preconditioner is: " << rvX << std::endl << std::endl;
821 
822   gmresDual.resetSolver();
823   gmresSym.resetSolver();
824   rvX = gmresDual(rMatSkylineDual, rvB, rvX0, pcLu, _real);
825   out << "The GmresSolver result with real Skyline dual non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl;
826   rvX = gmresSym(rMatSkylineSym, rvB, rvX0, pcLu, _real);
827   out << "The GmresSolver result with real Skyline sym non-sym matrix  and a lu factorized preconditioner is: " << rvX << std::endl << std::endl;
828 
829   gmresDual.resetSolver();
830   gmresSym.resetSolver();
831   rvX = gmresDual(rMatSkylineDual, rvB, rvX0, pcSorCsSym, _real);
832   out << "The GmresSolver result with real Skyline dual non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl;
833   rvX = gmresSym(rMatSkylineSym, rvB, rvX0, pcSorCsSym, _real);
834   out << "The GmresSolver result with real Skyline sym non-sym matrix  and a sor factorized preconditioner is: " << rvX << std::endl << std::endl;
835 */
836 /*
837   // II. Solver with a preconditioner
838   // Test with DENSE STORAGE
839   // Test with symmetric matrices
840   cvX = gmresRowCplx(cMatDenseRowSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
841   out << "The GmresSolver result with complex dense row sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
842   cvX = gmresColCplx(cMatDenseColSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
843   out << "The GmresSolver result with complex dense col sym matrix and a ldl* factorized preconditioner  is: " << cvX << std::endl;
844   cvX = gmresDualCplx(cMatDenseDualSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
845   out << "The GmresSolver result with complex dense dual sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
846   cvX = gmresSymCplx(cMatDenseSymSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
847   out << "The GmresSolver result with complex dense sym sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
848 
849   // Test with Sym SelfAjoint symmetric matrices
850   gmresRowCplx.resetSolver();
851   gmresColCplx.resetSolver();
852   gmresDualCplx.resetSolver();
853   gmresSymCplx.resetSolver();
854   cvX = gmresRowCplx(cMatDenseRowSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
855   out << "The GmresSolver result with complex dense row selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
856   cvX = gmresColCplx(cMatDenseColSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
857   out << "The GmresSolver result with complex dense col selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
858   cvX = gmresDualCplx(cMatDenseDualSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
859   out << "The GmresSolver result with complex dense dual selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
860   cvX = gmresSymCplx(cMatDenseSymSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
861   out << "The GmresSolver result with complex dense sym selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
862 
863   // Test with Skew Ajoint symmetric matrices
864   gmresRowCplx.resetSolver();
865   gmresColCplx.resetSolver();
866   gmresDualCplx.resetSolver();
867   gmresSymCplx.resetSolver();
868   cvX = gmresRowCplx(cMatDenseRowSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
869   out << "The GmresSolver result with complex dense row skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
870   cvX = gmresColCplx(cMatDenseColSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
871   out << "The GmresSolver result with complex dense col skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
872   cvX = gmresDualCplx(cMatDenseDualSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
873   out << "The GmresSolver result with complex dense dual skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
874   cvX = gmresSymCplx(cMatDenseSymSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
875   out << "The GmresSolver result with complex dense sym skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
876 
877   // Test with no symmetric matrices
878   gmresRowCplx.resetSolver();
879   gmresColCplx.resetSolver();
880   gmresDualCplx.resetSolver();
881   gmresSymCplx.resetSolver();
882   cvX = gmresRowCplx(cMatDenseRow, cvB, cvX0, pcLdlStar, xlifepp::_complex);
883   out << "The GmresSolver result with complex dense row non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
884   cvX = gmresColCplx(cMatDenseCol, cvB, cvX0, pcLdlStar, xlifepp::_complex);
885   out << "The GmresSolver result with complex dense col non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
886   cvX = gmresDualCplx(cMatDenseDual, cvB, cvX0, pcLdlStar, xlifepp::_complex);
887   out << "The GmresSolver result with complex dense dual non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
888   cvX = gmresSymCplx(cMatDenseSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
889   out << "The GmresSolver result with complex dense sym non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
890 
891   //---------------------------------
892   // Test with CS STORAGE
893   // Test with symmetric matrices
894   gmresRowCplx.resetSolver();
895   gmresColCplx.resetSolver();
896   gmresDualCplx.resetSolver();
897   gmresSymCplx.resetSolver();
898   cvX = gmresRowCplx(cMatCsRowSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
899   out << "The GmresSolver result with complex Cs row sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
900   cvX = gmresColCplx(cMatCsColSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
901   out << "The GmresSolver result with complex Cs col sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
902   cvX = gmresDualCplx(cMatCsDualSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
903   out << "The GmresSolver result with complex Cs dual sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
904   cvX = gmresSymCplx(cMatCsSymSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
905   out << "The GmresSolver result with complex Cs sym sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
906 
907   // Test with selfAjoint symmetric matrices
908   gmresRowCplx.resetSolver();
909   gmresColCplx.resetSolver();
910   gmresDualCplx.resetSolver();
911   gmresSymCplx.resetSolver();
912   cvX = gmresRowCplx(cMatCsRowSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
913   out << "The GmresSolver result with complex Cs row selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
914   cvX = gmresColCplx(cMatCsColSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
915   out << "The GmresSolver result with complex Cs col selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
916   cvX = gmresDualCplx(cMatCsDualSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
917   out << "The GmresSolver result with complex Cs dual selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
918   cvX = gmresSymCplx(cMatCsSymSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
919   out << "The GmresSolver result with complex Cs sym selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
920 
921   // Test with skewAjoint symmetric matrices
922   gmresRowCplx.resetSolver();
923   gmresColCplx.resetSolver();
924   gmresDualCplx.resetSolver();
925   gmresSymCplx.resetSolver();
926   cvX = gmresRowCplx(cMatCsRowSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
927   out << "The GmresSolver result with complex Cs row skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
928   cvX = gmresColCplx(cMatCsColSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
929   out << "The GmresSolver result with complex Cs col skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
930   cvX = gmresDualCplx(cMatCsDualSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
931   out << "The GmresSolver result with complex Cs dual skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
932   cvX = gmresSymCplx(cMatCsSymSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
933   out << "The GmresSolver result with complex Cs sym skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
934 
935   // Test with no symmetric matrices
936   gmresRowCplx.resetSolver();
937   gmresColCplx.resetSolver();
938   gmresDualCplx.resetSolver();
939   gmresSymCplx.resetSolver();
940   cvX = gmresRowCplx(cMatCsRow, cvB, cvX0, pcLdlStar, xlifepp::_complex);
941   out << "The GmresSolver result with complex Cs row non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
942   cvX = gmresColCplx(cMatCsCol, cvB, cvX0, pcLdlStar, xlifepp::_complex);
943   out << "The GmresSolver result with complex Cs col non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
944   cvX = gmresDualCplx(cMatCsDual, cvB, cvX0, pcLdlStar, xlifepp::_complex);
945   out << "The GmresSolver result with complex Cs dual non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
946   cvX = gmresSymCplx(cMatCsSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
947   out << "The GmresSolver result with complex Cs sym non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
948 
949   //---------------------------------
950   // Test with SKYLINE STORAGE
951   // Test with symmetric matrices
952   gmresDualCplx.resetSolver();
953   gmresSymCplx.resetSolver();
954   cvX = gmresDualCplx(cMatSkylineDualSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
955   out << "The GmresSolver result with complex Skyline dual sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
956   cvX = gmresSymCplx(cMatSkylineSymSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
957   out << "The GmresSolver result with complex Skyline sym sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
958 
959   // Test with self ajoint sym matrices
960   gmresDualCplx.resetSolver();
961   gmresSymCplx.resetSolver();
962   cvX = gmresDualCplx(cMatSkylineDualSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
963   out << "The GmresSolver result with complex Skyline dual selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
964   cvX = gmresSymCplx(cMatSkylineSymSymSelfAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
965   out << "The GmresSolver result with complex Skyline sym selfajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
966 
967   // Test with skew ajoint matrices
968   gmresDualCplx.resetSolver();
969   gmresSymCplx.resetSolver();
970   cvX = gmresDualCplx(cMatSkylineDualSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
971   out << "The GmresSolver result with complex Skyline dual skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
972   cvX = gmresSymCplx(cMatSkylineSymSymSkewAjoint, cvB, cvX0, pcLdlStar, xlifepp::_complex);
973   out << "The GmresSolver result with complex Skyline sym skewajoint sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
974 
975   // Test with no symmetric matrices
976   gmresDualCplx.resetSolver();
977   gmresSymCplx.resetSolver();
978   cvX = gmresDualCplx(cMatSkylineDual, cvB, cvX0, pcLdlStar, xlifepp::_complex);
979   out << "The GmresSolver result with complex Skyline dual non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl;
980   cvX = gmresSymCplx(cMatSkylineSym, cvB, cvX0, pcLdlStar, xlifepp::_complex);
981   out << "The GmresSolver result with complex Skyline sym non-sym matrix and a ldl* factorized preconditioner is: " << cvX << std::endl << std::endl;
982   */
983   //------------------------------------------------------------------------------------
984   // save results in a file or compare results with some references value in a file
985   //------------------------------------------------------------------------------------
986   trace_p->pop();
987   if (check) { return diffResFile(out, rootname); }
988   else { return saveResToFile(out, rootname); }
989 }
990 
991 }
992