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_LargeMatrixSkylineStorage.cpp
19 \author E. Lunéville
20 \since 04 jul 2012
21 \date 12 jul 2012
22
23 Low level tests of LargeMatrix class and SkylineStorage classes.
24 Almost functionalities are checked.
25 This function may either creates a reference file storing the results (check=false)
26 or compares results to those stored in the reference file (check=true)
27 It returns reporting informations in a string
28 */
29
30 #include "xlife++-libs.h"
31 #include "testUtils.hpp"
32
33 #include <iostream>
34 #include <fstream>
35 #include <vector>
36
37 using namespace xlifepp;
38
39 namespace unit_LargeMatrixSkylineStorage {
40
unit_LargeMatrixSkylineStorage(bool check)41 String unit_LargeMatrixSkylineStorage(bool check)
42 {
43 String rootname = "unit_LargeMatrixSkylineStorage";
44 trace_p->push(rootname);
45 std::stringstream out; // string stream receiving results
46 out.precision(testPrec);
47 verboseLevel(10);
48
49 Matrix<Real> m22(2, 1, -5.);
50
51 //EF numbering (regular Q1 of square)
52 Number n = 3;
53 std::vector< std::vector<Number> > elts((n - 1) * (n - 1), std::vector<Number>(4));
54 for(Number j = 1; j <= n - 1; j++)
55 for(Number i = 1; i <= n - 1; i++)
56 {
57 Number e = (j - 1) * (n - 1) + i - 1, p = (j - 1) * n + i;
58 elts[e][0] = p; elts[e][1] = p + 1; elts[e][2] = p + n; elts[e][3] = p + n + 1;
59 }
60 Number n2 = n * n;
61 std::vector< std::vector<Number> >::iterator it;
62 std::vector<Real> x1(n * n, 0.);
63 for(Number i = 0; i < n * n; i++) { x1[i] = i; }
64 out << "x1 = " << x1;
65
66 std::vector<Number> rs(4);rs[0]=1;rs[1]=2;rs[2]=8;rs[3]=9;
67
68 out << "\n----------------------------------- test sym skyline storage -----------------------------------\n";
69 SymSkylineStorage* skysym = new SymSkylineStorage(n * n, elts, elts);
70 skysym->visual(out);
71 LargeMatrix<Real> As1(skysym, 0., _symmetric);
72 //fill matrix
73 for(Number i = 1; i <= n * n; i++)
74 for(Number j = 1; j <= i; j++)
75 if(As1.pos(i, j) > 0) { As1(i, j) = 10 * i + j; }
76 out << "symsky real large matrix As1 : " << As1;
77 As1.viewStorage(out);
78 out << "product As1*x1 : " << As1* x1;
79 out << "\nproduct x1*As1 : " << x1* As1;
80 As1.saveToFile("Adense.mat", _dense);
81 out << "\nsave As1 to file Adense.mat in dense format\n";
82 As1.saveToFile("Acoo.mat", _coo);
83 out << "save As1 to file Acoo.mat in coordinates format\n";
84 LargeMatrix<Real> As1reload(n2, n2, _skyline, _symmetric, 0.);
85 As1reload.loadFromFile("Adense.mat", _dense);
86 out << "matrix As1reload from Adense : " << As1reload << "\n";
87 As1reload.loadFromFile("Acoo.mat", _coo);
88 out << "matrix As1reload from Acoo : " << As1reload << "\n";
89 SymSkylineStorage skysymb=SymSkylineStorage(*skysym);
90 out << "matrix storage updated with dense submatrix indices [1,2,8,9]x[1,2,8,9] :\n";
91 skysymb.visual(out);
92 skysymb.addSubMatrixIndices(rs,rs);
93 skysymb.visual(out);
94 out << "row 1 (col & adrs) : " <<As1.getRow(1) << "\n";
95 out << "row 2 (col & adrs) : " <<As1.getRow(2) << "\n";
96 out << "row 3 (col & adrs) : " <<As1.getRow(3) << "\n";
97 out << "row 9 (col & adrs) : " <<As1.getRow(9) << "\n";
98 out << "col 1 (col & adrs) : " <<As1.getCol(1) << "\n";
99 out << "col 2 (col & adrs) : " <<As1.getCol(2) << "\n";
100 out << "col 3 (col & adrs) : " <<As1.getCol(3) << "\n";
101 out << "col 9 (col & adrs) : " <<As1.getCol(9) << "\n";
102
103 out << "\n---------------test sym skyline storage and non symmetric matrix -------------------\n";
104 skysym->visual(out);
105 LargeMatrix<Real> Ans1(skysym, 0., _noSymmetry);
106 //fill matrix
107 for(Number i = 1; i <= n * n; i++)
108 for(Number j = 1; j <= i; j++)
109 if(Ans1.pos(i, j) > 0) {Ans1(i, j) = 10 * i + j; Ans1(j, i) = Ans1(i, j);}
110 out << "symsky real large non symmetric matrix Ans1 : " << Ans1;
111 Ans1.viewStorage(out);
112 out << "product Ans1*x1 : " << Ans1* x1;
113 out << "\nproduct x1*Ans1 : " << x1* Ans1;
114 Ans1.saveToFile("Adense.mat", _dense);
115 out << "\nsave Ans1 to file Adense.mat in dense format\n";
116 Ans1.saveToFile("Acoo.mat", _coo);
117 out << "save Ans1 to file Acoo.mat in coordinates format\n";
118 LargeMatrix<Real> Ans1reload(n2, n2, _skyline, _sym, 0.);
119 Ans1reload.loadFromFile("Adense.mat", _dense);
120 out << "matrix Ans1reload from Adense : " << Ans1reload << "\n";
121 Ans1reload.loadFromFile("Acoo.mat", _coo);
122 out << "matrix Ans1reload from Acoo : " << Ans1reload << "\n";
123
124 out << "\n\n----------------------------------- test dual skyline storage -----------------------------------\n";
125 DualSkylineStorage* skydual = new DualSkylineStorage(n * n, n * n, elts, elts);
126 skydual->visual(out);
127 LargeMatrix<Real> Ad1(skydual, 0.);
128 //fill matrix Aij='ij', and symmetric in fact
129 for(Number i = 1; i <= n * n; i++)
130 for(Number j = 1; j <= i; j++)
131 if(Ad1.pos(i, j) > 0) {Ad1(i, j) = 10 * i + j; Ad1(j, i) = Ad1(i, j);}
132 out << "dualsky real large matrix Ad1 : " << Ad1;
133 Ad1.viewStorage(out);
134 out << "product Ad1*x1 : " << Ad1* x1;
135 out << "\nproduct x1*Ad1 : " << x1* Ad1;
136 Ad1.saveToFile("Adense.mat", _dense);
137 out << "\nsave Ad1 to file Adense.mat in dense format\n";
138 Ad1.saveToFile("Acoo.mat", _coo);
139 out << "save Ad1 to file Acoo.mat in coordinates format\n";
140 LargeMatrix<Real> Ad1reload(n2, n2, _skyline, _dual, 0.);
141 Ad1reload.loadFromFile("Adense.mat", _dense);
142 out << "matrix Ad1reload from Adense : " << Ad1reload << "\n";
143 Ad1reload.loadFromFile("Acoo.mat", _coo);
144 out << "matrix Ad1reload from Acoo : " << Ad1reload << "\n";
145
146 DualSkylineStorage skydualb=DualSkylineStorage(*skydual);
147 out << "matrix storage updated with dense submatrix indices [1,2,8,9]x[1,2,8,9] :\n";
148 skydualb.visual(out);
149 skydualb.addSubMatrixIndices(rs,rs);
150 skydualb.visual(out);
151
152 out << "row 1 (col & adrs) : " <<Ad1.getRow(1) << "\n";
153 out << "row 2 (col & adrs) : " <<Ad1.getRow(2) << "\n";
154 out << "row 3 (col & adrs) : " <<Ad1.getRow(3) << "\n";
155 out << "row 9 (col & adrs) : " <<Ad1.getRow(9) << "\n";
156 out << "col 1 (col & adrs) : " <<Ad1.getCol(1) << "\n";
157 out << "col 2 (col & adrs) : " <<Ad1.getCol(2) << "\n";
158 out << "col 3 (col & adrs) : " <<Ad1.getCol(3) << "\n";
159 out << "col 9 (col & adrs) : " <<Ad1.getCol(9) << "\n";
160
161 const Number rowNum = 3;
162 const Number colNum = 3;
163
164 const std::string rMatrixDataSym("matrix3x3Sym.data");
165 const std::string rMatrixDataSymPosDef("matrix3x3SymPosDef.data");
166 const std::string rMatrixDataNoSym("matrix3x3NoSym.data");
167 const std::string cMatrixDataSymPosDef("cmatrix3x3SymPosDef.data");
168
169 LargeMatrix<Real> rMatSkylineSymSym(inputsPathTo(rMatrixDataSym), _dense, rowNum, _skyline, _symmetric);
170 LargeMatrix<Real> rMatSkylineSymNoSym(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, _skyline, _noSymmetry);
171 LargeMatrix<Real> rMatSkylineDualNoSym(inputsPathTo(rMatrixDataNoSym), _dense, rowNum, colNum, _skyline, _dual);
172 LargeMatrix<Real> rMatSkylineSymSymPosDef(inputsPathTo(rMatrixDataSymPosDef), _dense, rowNum, _skyline, _symmetric);
173 LargeMatrix<Real> rResSymSkylineLdlt(rMatSkylineSymSymPosDef);
174 LargeMatrix<Real> rResSymSkylineLu(rMatSkylineSymNoSym);
175 LargeMatrix<Real> rResDualSkylineLu(rMatSkylineDualNoSym);
176 Vector<Real> rVecB(3, 1.);
177 Vector<Real> rVecX(3, 1.);
178 LargeMatrix<Complex> cMatSkylineSymSymPosDef(inputsPathTo(cMatrixDataSymPosDef), _dense, rowNum, _skyline, _selfAdjoint);
179 LargeMatrix<Complex> cResLdlStar(cMatSkylineSymSymPosDef);
180 Vector<Complex> cVecB(3, 1.);
181 Vector<Complex> cVecX(3, 1.);
182
183 out << "\n----------------------------------- test LDLt and LU factorization -----------------------------------\n";
184 // I. Real Matrix
185 // SYMSKYLINE STORAGE
186 ldltFactorize(rResSymSkylineLdlt);
187 out << "The result of LDLt factorizing sym skyline matrix is " << rResSymSkylineLdlt << std::endl;
188 luFactorize(rResSymSkylineLu);
189 out << "The result LU factorizing sym skyline matrix is " << rResSymSkylineLu << std::endl;
190
191 // DUALSKYLINE STORAGE
192 luFactorize(rResDualSkylineLu);
193 out << "The result of LU factorizing dual skyline matrix is " << rResDualSkylineLu << std::endl;
194
195 // II. Complex Matrix
196 // SYMSKYLINE STORAGE
197 ldlstarFactorize(cResLdlStar);
198 out << "The result ldlstar factorized matrix is " << cResLdlStar << std::endl;
199
200 out << "----------------------------------- test solver with factorized matrix -----------------------------------\n";
201 rResSymSkylineLdlt.ldltSolve(rVecB, rVecX);
202 out << "The result LDLt solver of sym skyline matrix is " << rVecX << std::endl;
203 cResLdlStar.ldlstarSolve(cVecB, cVecX);
204 out << "The result LDL* solver of sym skyline matrix is " << cVecX << std::endl;
205 rResSymSkylineLu.luSolve(rVecB, rVecX);
206 out << "The result LU solver of sym skyline matrix is " << rVecX << std::endl;
207 rResDualSkylineLu.luSolve(rVecB, rVecX);
208 out << "The result LU solver of dual skyline matrix is " << rVecX << std::endl;
209
210 out << "----------------------------------- test diagonalMatrix -----------------------------------\n";
211 out << "Diagonal DualSkyline matrix : " << diagonalMatrix(Ad1, 3.) << std::endl;
212 out << "Diagonal SymSkyline matrix : " << diagonalMatrix(As1, 3.) << std::endl;
213
214 out << "----------------------------------- test multMatrixScalar -----------------------------------\n";
215 Real sigmaR = 3.0;
216 Complex sigmaC = Complex(3.0, 3.0);
217 out << "product Ad1*sigmaR : " << Ad1* sigmaR << std::endl;
218 out << "product sigmaR*Ad1 : " << sigmaR* Ad1 << std::endl;
219 out << "product Ad1*sigmaC : " << Ad1* sigmaC << std::endl;
220 out << "product sigmaC*Ad1 : " << sigmaC* Ad1 << std::endl;
221 out << "product As1*sigmaR : " << As1* sigmaR << std::endl;
222 out << "product sigmaR*As1 : " << sigmaR* As1 << std::endl;
223 out << "product As1*sigmaC : " << As1* sigmaC << std::endl;
224 out << "product sigmaC*As1 : " << sigmaC* As1 << std::endl;
225
226 out << "----------------------------------- test addMatrixMatrix -----------------------------------\n";
227 LargeMatrix<Real> rSkylineDual(Ad1); LargeMatrix<Real> rSkylineSym(As1);
228 out << "Sum of two dual skyline matrices : " << rSkylineDual + Ad1 << std::endl;
229 out << "Sum of two sym skyline matrices : " << rSkylineSym + As1 << std::endl;
230
231 out << "----------------------------- test factorize matrix * vector -------------------------------\n";
232 out<<"symSkyline factorization LLt : A*x-Af*x = "<<Reals(rMatSkylineSymSymPosDef*rVecB)-Reals(rResSymSkylineLdlt*rVecB)<<eol;
233 out<<"symSkyline factorization LU : A*x-Af*x = "<<Reals(rMatSkylineDualNoSym*rVecB)-Reals(rResSymSkylineLu*rVecB)<<eol;
234 out<<"dualSkyline factorization LU : A*x-Af*x = "<<Reals(rMatSkylineDualNoSym*rVecB)-Reals(rResDualSkylineLu*rVecB)<<eol;
235
236 //------------------------------------------------------------------------------------
237 // save results in a file or compare results with some references value in a file
238 //------------------------------------------------------------------------------------
239 trace_p->pop();
240 if (check) { return diffResFile(out, rootname); }
241 else { return saveResToFile(out, rootname); }
242 }
243
244 }
245