b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2matrix_compressed.C
1//------------------------------------------------------------------------
2// b2matrix_compressed.C --
3//
4// A framework for generic linear algebraic (matrix) computations.
5//
6// written by Mathias Doreille
7// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
8//
9// (c) 2004-2012 SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
13// Linder Höhe, 51147 Köln
14//
15// All Rights Reserved. Proprietary source code. The contents of
16// this file may not be disclosed to third parties, copied or
17// duplicated in any form, in whole or in part, without the prior
18// written permission of SMR.
19//------------------------------------------------------------------------
20
21#include "b2matrix_compressed.H"
22
23namespace b2000 {
24
25template <>
26void b2linalg::Matrix<double, b2linalg::Mcompressed_col_constref>::LUFactorization(
27 b2linalg::Matrix<double, b2linalg::Mcompressed_col>& trans_L,
28 b2linalg::Matrix<double, b2linalg::Mcompressed_col>& U, Index& P, Index& Q,
29 Vector<double, Vdense>& R) {
30 if (size1() == 0 || size2() == 0) {
31 trans_L.resize(size2(), size1(), 0);
32 U.resize(size1(), size2(), 0);
33 P.resize(size1());
34 Q.resize(size2());
35 R.resize(size1());
36 return;
37 }
38#ifdef HAVE_UMFPACK
39 double Control[UMFPACK_CONTROL];
40 umfpack_dl_defaults(Control);
41 // Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC;
42 Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE;
43 // Control[UMFPACK_FIXQ] = 1;
44 // Control[UMFPACK_DROPTOL] = 1e-100;
45 // Control[UMFPACK_PRL] = 6;
46 // Control[UMFPACK_PIVOT_TOLERANCE] = 1.0;
47 // Control[UMFPACK_DROPTOL] = 0.00001;
48 // Control[UMFPACK_AGGRESSIVE] = 0;
49 void* Symbolic;
50 int status;
51 if ((status = umfpack_dl_symbolic(
52 size1(), size2(), reinterpret_cast<long*>(const_cast<size_t*>(si)),
53 reinterpret_cast<long*>(const_cast<size_t*>(index)), const_cast<double*>(m),
54 &Symbolic, Control, 0))
55 != UMFPACK_OK) {
56 Exception() << "umfpack_dl_symbolic return the error code " << status << "." << THROW;
57 }
58 void* Numeric;
59 double info[UMFPACK_INFO];
60 if ((status = umfpack_dl_numeric(
61 reinterpret_cast<long*>(const_cast<size_t*>(si)),
62 reinterpret_cast<long*>(const_cast<size_t*>(index)), const_cast<double*>(m),
63 Symbolic, &Numeric, Control, info))
64 < UMFPACK_OK) {
65 Exception() << "umfpack_dl_numeric return the error code " << status << "." << THROW;
66 }
67
68 /*{
69 umfpack_dl_report_status(Control, status);
70 umfpack_dl_report_control(Control);
71 umfpack_dl_report_info(Control, info);
72 umfpack_dl_report_numeric(Numeric, Control);
73 }*/
74
75 umfpack_dl_free_symbolic(&Symbolic);
76 long lnz, unz, nrow, ncol, nz_udiag;
77 if ((status = umfpack_dl_get_lunz(&lnz, &unz, &nrow, &ncol, &nz_udiag, Numeric))
78 != UMFPACK_OK) {
79 Exception() << "umfpack_dl_get_lunz return the error code " << status << "." << THROW;
80 }
81 long diag = std::min(nrow, ncol);
82 trans_L.resize(diag, nrow, lnz);
83 U.resize(diag, ncol, unz);
84 P.resize(nrow);
85 Q.resize(ncol);
86 R.resize(nrow);
87 long do_recip;
88 if ((status = umfpack_dl_get_numeric(
89 reinterpret_cast<long*>(&trans_L.si[0]), reinterpret_cast<long*>(&trans_L.index[0]),
90 &trans_L.m[0], reinterpret_cast<long*>(&U.si[0]),
91 reinterpret_cast<long*>(&U.index[0]), &U.m[0], reinterpret_cast<long*>(&P[0]),
92 reinterpret_cast<long*>(&Q[0]), 0, &do_recip, &R[0], Numeric))
93 != UMFPACK_OK) {
94 Exception() << "umfpack_dl_get_numeric return the error code " << status << "." << THROW;
95 }
96
97 umfpack_dl_free_numeric(&Numeric);
98
99 if (nz_udiag != diag) { UMFPACK_clean_incomplete_factorization(trans_L, U, nz_udiag, P, Q, R); }
100
101 /*{
102 std::cout << "I = (" << std::endl;
103 for(int i = 0; i != nrow; ++i) {
104 std::cout << "(";
105 for(int j = 0; j != ncol; ++j)
106 std::cout << (*this)(P[i], Q[j]) << ", ";
107 std::cout << ")," << std::endl;
108 }
109 std::cout << ")" << std::endl;
110
111 std::cout << "L = (" << std::endl;
112 for(int i = 0; i != nrow; ++i) {
113 std::cout << "(";
114 for(int j = 0; j != nz_udiag; ++j)
115 std::cout << trans_L(j, i) << ", ";
116 std::cout << ")," << std::endl;
117 }
118 std::cout << ")" << std::endl;
119
120 std::cout << "U = (" << std::endl;
121 for(int i = 0; i != nz_udiag; ++i) {
122 std::cout << "(";
123 for(int j = 0; j != ncol; ++j)
124 std::cout << U(i, j) << ", ";
125 std::cout << ")," << std::endl;
126 }
127 std::cout << ")" << std::endl;
128 }*/
129
130 if (do_recip == 0) {
131 for (size_t i = 0; i != R.size(); ++i) { R[i] = 1.0 / R[i]; }
132 }
133#else
134 UnimplementedError() << "The umfpack library was not avaiable" << THROW;
135#endif
136}
137
138template <>
139void b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col_constref>::LUFactorization(
140 b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>& trans_L,
141 b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>& U, Index& P, Index& Q,
142 b2linalg::Vector<double, b2linalg::Vdense>& R) {
143 if (size1() == 0 || size2() == 0) {
144 trans_L.resize(size2(), size1(), 0);
145 U.resize(size1(), size2(), 0);
146 P.resize(size1());
147 Q.resize(size2());
148 R.resize(size1());
149 return;
150 }
151
152#ifdef HAVE_UMFPACK
153
154 double Control[UMFPACK_CONTROL];
155 umfpack_zl_defaults(Control);
156 // Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC;
157 Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE;
158 // Control[UMFPACK_FIXQ] = -1;
159 // Control[UMFPACK_PRL] = 6;
160 // Control[UMFPACK_PIVOT_TOLERANCE] = 1.0;
161 // Control[UMFPACK_DROPTOL] = 0.00001;
162 // Control[UMFPACK_AGGRESSIVE] = 0;
163 void* Symbolic;
164 int status;
165 if ((status = umfpack_zl_symbolic(
166 size1(), size2(), reinterpret_cast<long*>(const_cast<size_t*>(si)),
167 reinterpret_cast<long*>(const_cast<size_t*>(index)),
168 reinterpret_cast<double*>(const_cast<std::complex<double>*>(m)), 0, &Symbolic,
169 Control, 0))
170 != UMFPACK_OK) {
171 Exception() << "umfpack_zl_symbolic return the error code " << status << "." << THROW;
172 }
173 void* Numeric;
174 double info[UMFPACK_INFO];
175 if ((status = umfpack_zl_numeric(
176 reinterpret_cast<long*>(const_cast<size_t*>(si)),
177 reinterpret_cast<long*>(const_cast<size_t*>(index)),
178 reinterpret_cast<double*>(const_cast<std::complex<double>*>(m)), 0, Symbolic,
179 &Numeric, Control, info))
180 < UMFPACK_OK) {
181 Exception() << "umfpack_zl_numeric return the error code " << status << "." << THROW;
182 }
183
184 /*{
185 umfpack_zl_report_status(Control, status);
186 umfpack_zl_report_control(Control);
187 umfpack_zl_report_info(Control, info);
188 umfpack_zl_report_numeric(Numeric, Control);
189 }*/
190
191 umfpack_zl_free_symbolic(&Symbolic);
192 long lnz, unz, nrow, ncol, nz_udiag;
193 if ((status = umfpack_zl_get_lunz(&lnz, &unz, &nrow, &ncol, &nz_udiag, Numeric))
194 != UMFPACK_OK) {
195 Exception() << "umfpack_zl_get_lunz return the error code " << status << "." << THROW;
196 }
197 long diag = std::min(nrow, ncol);
198 trans_L.resize(diag, nrow, lnz);
199 U.resize(diag, ncol, unz);
200 P.resize(nrow);
201 Q.resize(ncol);
202 R.resize(nrow);
203 long do_recip;
204 if ((status = umfpack_zl_get_numeric(
205 reinterpret_cast<long*>(&trans_L.si[0]), reinterpret_cast<long*>(&trans_L.index[0]),
206 reinterpret_cast<double*>(&trans_L.m[0]), 0, reinterpret_cast<long*>(&U.si[0]),
207 reinterpret_cast<long*>(&U.index[0]), reinterpret_cast<double*>(&U.m[0]), 0,
208 reinterpret_cast<long*>(&P[0]), reinterpret_cast<long*>(&Q[0]), 0, 0, &do_recip,
209 &R[0], Numeric))
210 != UMFPACK_OK) {
211 Exception() << "umfpack_zl_get_numeric return the error code " << status << "." << THROW;
212 }
213
214 umfpack_zl_free_numeric(&Numeric);
215
216 if (nz_udiag != diag) { UMFPACK_clean_incomplete_factorization(trans_L, U, nz_udiag, P, Q, R); }
217
218 /*{
219 std::cout << "I = (" << std::endl;
220 for(int i = 0; i != nrow; ++i) {
221 std::cout << "(";
222 for(int j = 0; j != ncol; ++j)
223 std::cout << (*this)(P[i], Q[j]) << ", ";
224 std::cout << ")," << std::endl;
225 }
226 std::cout << ")" << std::endl;
227
228 std::cout << "L = (" << std::endl;
229 for(int i = 0; i != nrow; ++i) {
230 std::cout << "(";
231 for(int j = 0; j != nz_udiag; ++j)
232 std::cout << trans_L(j, i) << ", ";
233 std::cout << ")," << std::endl;
234 }
235 std::cout << ")" << std::endl;
236
237 std::cout << "U = (" << std::endl;
238 for(int i = 0; i != nz_udiag; ++i) {
239 std::cout << "(";
240 for(int j = 0; j != ncol; ++j)
241 std::cout << U(i, j) << ", ";
242 std::cout << ")," << std::endl;
243 }
244 std::cout << ")" << std::endl;
245 }*/
246
247 if (do_recip == 0) {
248 for (size_t i = 0; i != R.size(); ++i) { R[i] = 1.0 / R[i]; }
249 }
250#else
251 UnimplementedError() << "The umfpack library is not available." << THROW;
252#endif
253}
254
255template <>
256void b2linalg::Matrix<b2000::csda<double>, b2linalg::Mcompressed_col_constref>::LUFactorization(
257 b2linalg::Matrix<b2000::csda<double>, b2linalg::Mcompressed_col>& trans_L,
258 b2linalg::Matrix<b2000::csda<double>, b2linalg::Mcompressed_col>& U, Index& P, Index& Q,
259 b2linalg::Vector<double, b2linalg::Vdense>& R) {
260 b2linalg::Matrix<std::complex<double>, Mcompressed_col_constref>* c_this =
261 reinterpret_cast<b2linalg::Matrix<std::complex<double>, Mcompressed_col_constref>*>(this);
262 c_this->LUFactorization(
263 reinterpret_cast<b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>&>(
264 trans_L),
265 reinterpret_cast<b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col>&>(U),
266 P, Q, R);
267}
268
269template <typename T>
270void b2linalg::Matrix<T, b2linalg::Mcompressed_col_constref>::
271 UMFPACK_clean_incomplete_factorization(
272 b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_L,
273 b2linalg::Matrix<T, b2linalg::Mcompressed_col>& U, const size_t nb_udiag,
274 b2linalg::Index& P, b2linalg::Index& Q, b2linalg::Vector<double, b2linalg::Vdense>& R) {
275 Index perm(U.size1());
276 {
277 size_t ptr_d = 0;
278 size_t ptr_0 = nb_udiag;
279 for (size_t i = 0; i != U.size1(); ++i) {
280 if (U.si[i + 1] - U.si[i] == 0 || U.index[U.si[i + 1] - 1] != i) {
281 perm[ptr_0++] = i;
282 } else {
283 perm[ptr_d++] = i;
284 }
285 }
286 }
287 Index inv_perm = perm.make_dual();
288
289 // U clean
290 {
291 size_t i = 0;
292 for (size_t j = 0; j != U.size2(); ++j) {
293 for (const size_t i_end = U.si[j + 1]; i != i_end; ++i) {
294 size_t& tmp = U.index[i];
295 tmp = inv_perm[tmp];
296 }
297 }
298
299 std::vector<size_t> si_tmp(U.size1() + 1);
300 si_tmp[0] = 0;
301 std::vector<size_t> index_tmp(U.si[U.size1()]);
302 std::vector<T> value_tmp(index_tmp.size());
303 for (size_t j = 0; j != perm.size(); ++j) {
304 const size_t perm_j = perm[j];
305 const size_t begin = U.si[perm_j];
306 const size_t end = U.si[perm_j + 1];
307 std::copy(
308 U.index.begin() + begin, U.index.begin() + end, index_tmp.begin() + si_tmp[j]);
309 std::copy(U.m.begin() + begin, U.m.begin() + end, value_tmp.begin() + si_tmp[j]);
310 si_tmp[j + 1] = si_tmp[j] + (end - begin);
311 }
312 std::copy(si_tmp.begin(), si_tmp.end(), U.si.begin());
313 std::copy(index_tmp.begin(), index_tmp.end(), U.index.begin());
314 std::copy(value_tmp.begin(), value_tmp.end(), U.m.begin());
315 }
316 U.s1 = nb_udiag;
317
318 // trans L clean
319 {
320 std::vector<size_t> si_tmp(trans_L.size2() + 1);
321 si_tmp[0] = 0;
322
323 {
324 size_t i = 0;
325 size_t i_w = i;
326 for (size_t j = 0; j != trans_L.size2();) {
327 ++j;
328 for (const size_t i_end = trans_L.si[j]; i != i_end; ++i) {
329 const size_t tmp = inv_perm[trans_L.index[i]];
330 if (tmp < nb_udiag) {
331 trans_L.index[i_w] = tmp;
332 trans_L.m[i_w] = trans_L.m[i];
333 ++i_w;
334 }
335 }
336 si_tmp[j] = i_w;
337 }
338 }
339
340 std::vector<size_t> index_tmp(si_tmp[trans_L.size2()]);
341 std::vector<T> value_tmp(index_tmp.size());
342 for (size_t j = 0; j != perm.size(); ++j) {
343 const size_t perm_j = perm[j];
344 const size_t begin = si_tmp[perm_j];
345 const size_t end = si_tmp[perm_j + 1];
346 std::copy(
347 trans_L.index.begin() + begin, trans_L.index.begin() + end,
348 index_tmp.begin() + trans_L.si[j]);
349 std::copy(
350 trans_L.m.begin() + begin, trans_L.m.begin() + end,
351 value_tmp.begin() + trans_L.si[j]);
352 trans_L.si[j + 1] = trans_L.si[j] + (end - begin);
353 }
354 std::copy(si_tmp.begin() + perm.size(), si_tmp.end(), trans_L.si.begin() + perm.size());
355 {
356 const size_t i = si_tmp[perm.size()];
357 const size_t i_end = si_tmp.back();
358 std::copy(
359 trans_L.index.begin() + i, trans_L.index.begin() + i_end, index_tmp.begin() + i);
360 std::copy(trans_L.m.begin() + i, trans_L.m.begin() + i_end, value_tmp.begin() + i);
361 }
362 trans_L.index.swap(index_tmp);
363 trans_L.m.swap(value_tmp);
364 }
365 trans_L.s1 = nb_udiag;
366
367 // P, Q permutation
368 {
369 std::vector<size_t> tmp(perm.size());
370 for (size_t i = 0; i != perm.size(); ++i) { tmp[i] = P[perm[i]]; }
371 std::copy(tmp.begin(), tmp.begin() + perm.size(), P.begin());
372 for (size_t i = 0; i != perm.size(); ++i) { tmp[i] = Q[perm[i]]; }
373 std::copy(tmp.begin(), tmp.begin() + perm.size(), Q.begin());
374 }
375}
376
377} // namespace b2000
THROW
#define THROW
Definition b2exception.H:198
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::UnimplementedError
GenericException< UnimplementedError_name > UnimplementedError
Definition b2exception.H:314