b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2matrix_operation.H
1//------------------------------------------------------------------------
2// b2matrix_operation.H --
3//
4// A framework for generic linear algebraic (matrix) computations.
5//
6// written by Mathias Doreille
7// Thomas Blome <thomas.blome@dlr.de>
8//
9// Copyright (c) 2004-2012 SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// Copyright (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#ifndef __B2MATRIX_OPERATION_H__
22#define __B2MATRIX_OPERATION_H__
23
24#include <filesystem>
25#include <memory>
26
27#include "b2linear_algebra_def.H"
28#include "b2linear_algebra_utils.H"
29
30namespace b2000::b2linalg {
31
32template <typename T, typename STORAGE1, typename STORAGE2>
33class Matrix<T, Sum<STORAGE1, STORAGE2> > {
34public:
35 Matrix(const Matrix<T, STORAGE1>& m1_, const Matrix<T, STORAGE2>& m2_, T scale_ = T(1))
36 : m1(m1_), m2(m2_), scale(scale_), trans(false) {
37 if (m1.size() != m2.size()) {
38 Exception() << "The two matrices do not have the same size." << THROW;
39 }
40 }
41
42 std::pair<size_t, size_t> size() const {
43 if (trans) {
44 return std::pair<size_t, size_t>(m1.size2(), m1.size1());
45 } else {
46 return m1.size();
47 }
48 }
49
50 size_t size1() const {
51 if (trans) {
52 return m1.size2();
53 } else {
54 return m1.size1();
55 }
56 }
57
58 size_t size2() const {
59 if (trans) {
60 return m1.size1();
61 } else {
62 return m1.size2();
63 }
64 }
65
66 T operator()(size_t i, size_t j) const {
67 if (trans) {
68 return scale * (m1(j, i) + m2(j, i));
69 } else {
70 return scale * (m1(i, j) + m2(i, j));
71 }
72 }
73
74 Matrix& operator*(T t) {
75 scale *= t;
76 return *this;
77 }
78
79 Matrix& transpose() {
80 trans = trans ? false : true;
81 return *this;
82 }
83
84private:
85 Matrix<T, STORAGE1> m1;
86 Matrix<T, STORAGE2> m2;
87 T scale;
88 char trans;
89 MVFRIEND;
90};
91
92template <typename T1, typename T2>
93struct MMProd {
94 using const_base = MMProd;
95};
96
97template <typename T, typename STORAGE1, typename STORAGE2>
98class Matrix<T, MMProd<STORAGE1, STORAGE2> > {
99public:
100 Matrix(const Matrix<T, STORAGE1>& m1_, const Matrix<T, STORAGE2>& m2_, T scale_ = T(1))
101 : m1(m1_), m2(m2_), scale(scale_), trans(false) {
102 if (m1.size2() != m2.size1()) {
103 Exception() << "Matrix-matrix inner product: The number of columns of the "
104 "first matrix ("
105 << m1.size1() << " x " << m1.size2()
106 << ") differs from the number of rows of the second matrix (" << m2.size1()
107 << " x " << m2.size2() << ")." << THROW;
108 }
109 }
110
111 std::pair<size_t, size_t> size() const {
112 return std::pair<size_t, size_t>(m1.size1(), m2.size2());
113 }
114
115 size_t size1() const { return m1.size1(); }
116
117 size_t size2() const { return m2.size2(); }
118
119 T operator()(size_t i, size_t j) const {
120 if (trans) {
121 return scale * transpose(m1)[i] * m2[j];
122 } else {
123 return scale * transpose(m1)[j] * m2[i];
124 }
125 }
126
127 Matrix& operator*(T t) {
128 scale *= t;
129 return *this;
130 }
131
132 Matrix& transpose() {
133 trans = trans ? false : true;
134 return *this;
135 }
136
137private:
138 Matrix<T, STORAGE1> m1;
139 Matrix<T, STORAGE2> m2;
140 T scale;
141 char trans;
142 MVFRIEND;
143};
144
145template <typename T1, typename S1, typename T2, typename S2>
146inline void copy_m(const Matrix<T1, S1>& a, Matrix<T2, S2>& b) {
147 b.resize(a.size1(), a.size2());
148 for (size_t i = 0; i != a.size1(); ++i) {
149 for (size_t j = 0; j != a.size2(); ++j) { b(i, j) = a(i, j); }
150 }
151}
152
153template <typename T1, typename S1, typename T2, typename S2>
154inline void copy_vm(const std::vector<Matrix<T1, S1> >& a, std::vector<Matrix<T2, S2> >& b) {
155 b.resize(a.size());
156 for (size_t i = 0; i != a.size(); ++i) { copy_m(a[i], b[i]); }
157}
158
159template <typename T, typename STORAGE>
160Matrix<T, typename STORAGE::const_base> transposed(const Matrix<T, STORAGE>& m) {
161 return Matrix<T, typename STORAGE::const_base>(m).transpose();
162}
163
164template <typename T, typename STORAGE1, typename STORAGE2>
165Matrix<T, Sum<typename STORAGE1::const_base, typename STORAGE2::const_base> > operator+(
166 const Matrix<T, STORAGE1>& m1, const Matrix<T, STORAGE2>& m2) {
167 return Matrix<T, Sum<typename STORAGE1::const_base, typename STORAGE2::const_base> >(m1, m2);
168}
169
170template <typename T, typename STORAGE1, typename STORAGE2>
171Matrix<T, Sum<typename STORAGE1::const_base, typename STORAGE2::const_base> > operator-(
172 const Matrix<T, STORAGE1>& m1, const Matrix<T, STORAGE2>& m2) {
173 return Matrix<T, Sum<typename STORAGE1::const_base, typename STORAGE2::const_base> >(m1, -m2);
174}
175
176template <typename T, typename STORAGE1, typename STORAGE2>
177Matrix<T, MMProd<typename STORAGE1::const_base, typename STORAGE2::const_base> > operator*(
178 const Matrix<T, STORAGE1>& m1, const Matrix<T, STORAGE2>& m2) {
179 return Matrix<T, MMProd<typename STORAGE1::const_base, typename STORAGE2::const_base> >(m1, m2);
180}
181
182template <typename T, typename STORAGE>
183Matrix<T, typename STORAGE::const_base> operator*(T t, const Matrix<T, STORAGE>& m) {
184 return Matrix<T, typename STORAGE::const_base>(m) * t;
185}
186
187template <typename T, typename STORAGE>
188Matrix<T, typename STORAGE::const_base> operator*(const Matrix<T, STORAGE>& m, T t) {
189 return Matrix<T, typename STORAGE::const_base>(m) * t;
190}
191
192template <typename T, typename STORAGE>
193Matrix<T, typename STORAGE::const_base> operator-(const Matrix<T, STORAGE>& m) {
194 return Matrix<T, typename STORAGE::const_base>(m) * T(-1);
195}
196
197template <typename MTYPE, typename INDEX1, typename INDEX2>
198struct Msub_constref {
199 using base = Msub_constref;
200 using const_base = Msub_constref;
201};
202
203template <typename T, typename MTYPE, typename INDEX1, typename INDEX2>
204class Matrix<T, Msub_constref<MTYPE, INDEX1, INDEX2> > {
205public:
206 Matrix(const Matrix<T, MTYPE> m_, const INDEX1& index1_, const INDEX2& index2_)
207 : m(m_), index1(index1_), index2(index2_) {}
208
209 std::pair<size_t, size_t> size() const {
210 return std::pair<size_t, size_t>(
211 index1.is_all() ? m.size1() : index1.size(),
212 index2.is_all() ? m.size2() : index2.size());
213 }
214
215 size_t size1() const { return index1.is_all() ? m.size1() : index1.size(); }
216
217 size_t size2() const { return index2.is_all() ? m.size2() : index2.size(); }
218
219 T operator()(size_t i, size_t j) const {
220 return m(index1.is_all() ? i : index1[i], index2.is_all() ? j : index2[j]);
221 }
222
223private:
224 MVFRIEND;
225 const Matrix<T, MTYPE> m;
226 const INDEX1& index1;
227 const INDEX2& index2;
228};
229
230template <typename MTYPE>
231struct Minverse_constref {
232 using base = Minverse_constref;
233 using const_base = Minverse_constref;
234};
235
236template <typename T, typename MTYPE>
237class Matrix<T, Minverse_constref<MTYPE> > {
238public:
239 Matrix(
240 const Matrix<T, MTYPE>& m_, Matrix<T>& null_space_ = Matrix<T>::null,
241 ssize_t max_null_space_vector_ = 0)
242 : m(m_), null_space(null_space_), max_null_space_vector(max_null_space_vector_) {}
243
244 Matrix(const Matrix& m_)
245 : m(m_.m), null_space(m_.null_space), max_null_space_vector(m_.max_null_space_vector) {}
246
247 std::pair<size_t, size_t> size() const { return m.size(); }
248
249 size_t size1() const { return m.size1(); }
250
251 size_t size2() const { return m.size2(); }
252
253private:
254 const Matrix<T, MTYPE>& m;
255 Matrix<T>& null_space;
256 ssize_t max_null_space_vector;
257 MVFRIEND;
258};
259
260template <typename T, typename MTYPE>
261Matrix<T, Minverse_constref<MTYPE> > inverse(
262 const Matrix<T, MTYPE>& m, Matrix<T, Mrectangle>& null_space = Matrix<T, Mrectangle>::null,
263 ssize_t max_null_space_vector = 0) {
264 return Matrix<T, Minverse_constref<MTYPE> >(m, null_space, max_null_space_vector);
265}
266
267template <typename MTYPE>
268struct Minverse_ext_constref {
269 using base = Minverse_ext_constref;
270 using const_base = Minverse_ext_constref;
271};
272
273template <typename T, typename MTYPE>
274class Matrix<T, Minverse_ext_constref<MTYPE> > {
275public:
276 Matrix(
277 const Matrix<T, MTYPE>& m_, const T* a_, const T* b_, const T* c_,
278 Matrix<T>& null_space_ = Matrix<T>::null, ssize_t max_null_space_vector_ = 0)
279 : m(m_),
280 a(a_),
281 b(b_),
282 c(c_),
283 null_space(null_space_),
284 max_null_space_vector(max_null_space_vector_) {}
285
286 Matrix(const Matrix& m_)
287 : m(m_.m),
288 a(m_.a),
289 b(m_.b),
290 c(m_.c),
291 null_space(m_.null_space),
292 max_null_space_vector(m_.max_null_space_vector) {}
293
294 std::pair<size_t, size_t> size() const { return m.size(); }
295
296 size_t size1() const { return m.size1(); }
297
298 size_t size2() const { return m.size2(); }
299
300private:
301 const Matrix<T, MTYPE>& m;
302 const T* a;
303 const T* b;
304 const T* c;
305 Matrix<T>& null_space;
306 ssize_t max_null_space_vector;
307 MVFRIEND;
308};
309
310template <typename T, typename MTYPE, typename MTYPEA, typename MTYPEB, typename MTYPEC>
311Matrix<T, Minverse_ext_constref<MTYPE> > update_extension_and_inverse(
312 const Matrix<T, MTYPE>& m, const Matrix<T, MTYPEA>& a, const Matrix<T, MTYPEB>& b,
313 const Matrix<T, MTYPEC>& c, Matrix<T, Mrectangle>& null_space = Matrix<T, Mrectangle>::null,
314 ssize_t max_null_space_vector = 0) {
315 return Matrix<T, Minverse_ext_constref<MTYPE> >(
316 m, &a(0, 0), &b(0, 0), &c(0, 0), null_space, max_null_space_vector);
317}
318
319template <typename T, typename MTYPE, typename MTYPEA, typename MTYPEB>
320Matrix<T, Minverse_ext_constref<MTYPE> > update_extension_and_inverse(
321 const Matrix<T, MTYPE>& m, const Vector<T, MTYPEA>& a, const Vector<T, MTYPEB>& b, const T& c,
322 Matrix<T, Mrectangle>& null_space = Matrix<T, Mrectangle>::null,
323 ssize_t max_null_space_vector = 0) {
324 return Matrix<T, Minverse_ext_constref<MTYPE> >(
325 m, &a[0], &b[0], &c, null_space, max_null_space_vector);
326}
327
329
347template <typename T, typename TYPE>
348void print(
349 const Matrix<T, TYPE>& m, std::ostream& out, std::string seperator = ", ",
350 std::string row_sep = "\n",
351 std::pair<std::string, std::string> row_enclosure =
352 std::pair<std::string, std::string>("", ""),
353 std::pair<std::string, std::string> mat_enclosure =
354 std::pair<std::string, std::string>("", "")) {
355 size_t n_rows = m.size1();
356 size_t n_cols = m.size2();
357
358 out << mat_enclosure.first;
359 for (size_t i = 0; i < n_rows; ++i) {
360 out << row_enclosure.first;
361 for (size_t j = 0; j < n_cols - 1; ++j) { out << m(i, j) << seperator; }
362 // finish the row with row_enclosure instead of seperator
363 out << m(i, n_cols - 1) << row_enclosure.second;
364 // add the row seperator except for the last row
365 if (i < n_rows - 1) { out << row_sep; }
366 }
367 out << mat_enclosure.second;
368}
369
371template <typename T, typename TYPE>
372std::ostream& operator<<(std::ostream& out, const Matrix<T, TYPE>& m) {
373 const auto closures = std::pair<std::string, std::string>("( ", " )");
374 print(m, out, ", ", ",\n", closures, closures);
375 return out;
376}
377
379
384template <typename T, typename TYPE>
385void export_csv(
386 const Matrix<T, TYPE>& m, std::string matrixname, std::filesystem::path filepath,
387 std::streamsize precision = 16) {
388 std::ofstream out(filepath);
389 out.precision(precision);
390 out << "# " << matrixname << " with (" << m.size1() << " x " << m.size2() << ") DoF\n";
391 print(m, out);
392 out << std::endl;
393}
394
395template <typename ATYPE, typename BTYPE>
396struct Eigenvector {
397 using base = Eigenvector;
398 using const_base = Eigenvector;
399};
400
401template <typename T, typename ATYPE, typename BTYPE>
402class Matrix<T, Eigenvector<ATYPE, BTYPE> > {
403public:
404 Matrix(
405 const Matrix<T, ATYPE>& a_, const char type_a_, const Matrix<T, BTYPE>& b_,
406 const char type_b_, Vector<T, Vdense_ref> eigenvalue_, const char which_[2], T sigma_,
407 int ncv_, double tol_, int maxit_, const Vector<T, Vdense_constref> starting_eigenvector_)
408 : a(a_),
409 type_a(type_a_),
410 b(b_),
411 type_b(type_b_),
412 eigenvalue(eigenvalue_),
413 sigma(sigma_),
414 ncv(ncv_),
415 tol(tol_),
416 maxit(maxit_),
417 starting_eigenvector(starting_eigenvector_) {
418 if (!b.is_null() && a.size() != b.size()) {
419 Exception() << "Eigenvalue extraction: Matrix A and B do not have same size." << THROW;
420 }
421 if (!starting_eigenvector.is_null() && starting_eigenvector.size() != a.size1()) {
422 Exception() << "Eigenvalue extraction: Eigenvalue start vector size != matrix size."
423 << THROW;
424 }
425 strcpy(which, which_);
426 if (tol < 0) { tol = 0; }
427 }
428
429 std::pair<size_t, size_t> size() const {
430 return std::pair<size_t, size_t>(a.size1(), eigenvalue.size());
431 }
432
433 size_t size1() const { return a.size1(); }
434
435 size_t size2() const { return eigenvalue.size(); }
436
437 void resolve(Matrix<T, Mrectangle_increment_ref>& eigenvector);
438
439private:
440 const Matrix<T, ATYPE>& a;
441 const char type_a;
442 const Matrix<T, BTYPE>& b;
443 const char type_b;
444 Vector<T, Vdense_ref> eigenvalue;
445 char which[3];
446 T sigma;
447 int ncv;
448 double tol;
449 int maxit;
450 Vector<T, Vdense_constref> starting_eigenvector;
451
452 struct OP {
453 virtual ~OP() {}
454 virtual void op(
455 Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {}
456 virtual void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) {}
457 virtual void end(
458 Vector<T, Vdense_ref> eigenvalue, Matrix<T, Mrectangle_increment_ref> eigenvector) {}
459 };
460
461 struct OPS1 : public OP {
462 OPS1(const Matrix<T, ATYPE>& a_) : a(a_) {}
463 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
464 y = a * x;
465 }
466 const Matrix<T, ATYPE>& a;
467 };
468
469 struct OPS1i : public OP {
470 OPS1i(const Matrix<T, ATYPE>& a_) : a(a_) {}
471 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
472 y = inverse(a) * x;
473 }
474 void end(
475 Vector<T, Vdense_ref> eigenvalue, Matrix<T, Mrectangle_increment_ref> eigenvector) {
476 for (size_t i = 0; i != eigenvalue.size(); ++i) {
477 eigenvalue[i] = T(1) / eigenvalue[i];
478 }
479 }
480 const Matrix<T, ATYPE>& a;
481 };
482
483 struct OPS2 : public OP {
484 OPS2(const Matrix<T, ATYPE>& a_, const Matrix<T, BTYPE>& b_) : a(a_), b(b_) {}
485 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
486 y = a * x;
487 y = inverse(b) * y;
488 }
489 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) { y = b * x; }
490 const Matrix<T, ATYPE>& a;
491 const Matrix<T, BTYPE>& b;
492 };
493
494 struct OPS2i : public OP {
495 OPS2i(const Matrix<T, ATYPE>& a_, const Matrix<T, BTYPE>& b_) : a(a_), b(b_) {}
496 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
497 y = b * x;
498 y = inverse(a) * y;
499 }
500 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) { y = a * x; }
501 void end(
502 Vector<T, Vdense_ref> eigenvalue, Matrix<T, Mrectangle_increment_ref> eigenvector) {
503 for (size_t i = 0; i != eigenvalue.size(); ++i) {
504 eigenvalue[i] = T(1) / eigenvalue[i];
505 }
506 }
507 const Matrix<T, ATYPE>& a;
508 const Matrix<T, BTYPE>& b;
509 };
510
511 struct OPS3_0 : public OP {
512 OPS3_0(const Matrix<T, ATYPE>& a_, const Matrix<T, BTYPE>& b_) : a(a_), b(b_) {}
513 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
514 if (Bx.is_null()) {
515 y = b * x;
516 y = inverse(a) * y;
517 } else {
518 y = inverse(a) * Bx;
519 }
520 }
521 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) { y = b * x; }
522 const Matrix<T, ATYPE>& a;
523 const Matrix<T, BTYPE>& b;
524 };
525
526 struct OPS3 : public OP {
527 OPS3(const Matrix<T, ATYPE>& a, const Matrix<T, BTYPE>& b_, T sigma) : b(b_), m(a) {
528 if (b.is_null()) {
529 for (size_t i = 0; i != a.size1(); ++i) { m(i, i) -= sigma; }
530 } else {
531 m += -sigma * b;
532 }
533 }
534 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
535 if (Bx.is_null()) {
536 y = b * x;
537 y = inverse(m) * y;
538 } else {
539 y = inverse(m) * Bx;
540 }
541 }
542 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) { y = b * x; }
543 const Matrix<T, BTYPE>& b;
544 Matrix<T, typename ATYPE::inverse> m;
545 };
546
547 struct OPS4 : public OP {
548 OPS4(const Matrix<T, ATYPE>& a_, const Matrix<T, BTYPE>& b_, T sigma) : a(a_), b(b_), m(a) {
549 if (b.is_null()) {
550 for (size_t i = 0; i != a.size1(); ++i) { m(i, i) -= sigma; }
551 } else {
552 m += -sigma * b;
553 }
554 }
555 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
556 if (Bx.is_null()) {
557 y = b * x;
558 y = inverse(m) * y;
559 } else {
560 y = inverse(m) * Bx;
561 }
562 }
563 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) { y = a * x; }
564 const Matrix<T, ATYPE>& a;
565 const Matrix<T, BTYPE>& b;
566 Matrix<T, typename ATYPE::inverse> m;
567 };
568
569 struct OPS5 : public OP {
570 OPS5(const Matrix<T, ATYPE>& a_, const Matrix<T, BTYPE>& b_, T sigma_)
571 : sigma(sigma_), a(a_), b(b_), m(a) {
572 if (b.is_null()) {
573 for (size_t i = 0; i != a.size1(); ++i) { m(i, i) -= sigma; }
574 } else {
575 m += -sigma * b;
576 }
577 }
578 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
579 if (Bx.is_null()) {
580 y = b * x;
581 y += a * sigma * y;
582 } else {
583 y = a * sigma * Bx;
584 }
585 y = inverse(m) * y;
586 }
587 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) { y = b * x; }
588 T sigma;
589 const Matrix<T, ATYPE>& a;
590 const Matrix<T, BTYPE>& b;
591 Matrix<T, typename ATYPE::inverse> m;
592 };
593
594 struct OPNN1 : public OP {
595 OPNN1(const Matrix<T, ATYPE>& a_, const Matrix<T, BTYPE>& b_) : a(a_), b(b_) {}
596 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
597 y = b * x;
598 y = inverse(a) * y;
599 }
600 void end(
601 Vector<T, Vdense_ref> eigenvalue, Matrix<T, Mrectangle_increment_ref> eigenvector) {
602 for (size_t i = 0; i != eigenvalue.size(); ++i) {
603 eigenvalue[i] = T(1) / eigenvalue[i];
604 }
605 }
606 const Matrix<T, ATYPE>& a;
607 const Matrix<T, BTYPE>& b;
608 };
609
610 struct OPNN1_shift : public OP {
611 OPNN1_shift(const Matrix<T, ATYPE>& a_, const Matrix<T, BTYPE>& b_, T sigma_)
612 : sigma(sigma_), a(a_), b(b_), m(a) {
613 if (b.is_null()) {
614 for (size_t i = 0; i != a.size1(); ++i) { m(i, i) -= sigma; }
615 } else {
616 m += -sigma * b;
617 }
618 }
619 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
620 y = b * x;
621 y = inverse(m) * y;
622 }
623 void end(
624 Vector<T, Vdense_ref> eigenvalue, Matrix<T, Mrectangle_increment_ref> eigenvector) {
625 std::vector<std::pair<double, size_t> > tmp;
626 for (size_t i = 0; i != eigenvalue.size(); ++i) {
627 eigenvalue[i] = sigma + T(1) / eigenvalue[i];
628 tmp.push_back(std::pair<double, size_t>(b2000::real(eigenvalue[i]), tmp.size()));
629 }
630 std::sort(tmp.begin(), tmp.end(), CompareFirstOfPair());
631 Vector<T, Vdense> eigenvalue_tmp(eigenvalue.size());
632 Matrix<T, Mrectangle> eigenvector_tmp(eigenvector.size1(), eigenvector.size2());
633 for (size_t i = 0; i != eigenvalue.size(); ++i) {
634 eigenvalue_tmp[i] = eigenvalue[tmp[i].second];
635 eigenvector_tmp[i] = eigenvector[tmp[i].second];
636 }
637 eigenvalue = eigenvalue_tmp;
638 eigenvector = Matrix<T, Mrectangle_increment_constref>(eigenvector_tmp);
639 }
640 T sigma;
641 const Matrix<T, ATYPE>& a;
642 const Matrix<T, BTYPE>& b;
643 Matrix<T, typename ATYPE::inverse> m;
644 };
645
646 struct OPNN2 : public OP {
647 OPNN2(const Matrix<T, ATYPE>& a_, const Matrix<T, BTYPE>& b_) : a(a_), b(b_) {}
648 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
649 y = a * x;
650 y = inverse(b) * y;
651 y = inverse(b) * y;
652 }
653 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) {
654 y = b * transposed(b) * x;
655 }
656 void end(
657 Vector<T, Vdense_ref> eigenvalue, Matrix<T, Mrectangle_increment_ref> eigenvector) {
658 Matrix<T, Mrectangle> tmp;
659 tmp = eigenvector;
660 eigenvector = transposed(b) * tmp;
661 }
662 const Matrix<T, ATYPE>& a;
663 const Matrix<T, BTYPE>& b;
664 };
665
666 struct OPNN3 : public OP {
667 OPNN3(const Matrix<T, ATYPE>& a, const Matrix<T, BTYPE>& b_, T sigma) : b(b_), m(a) {
668 if (b.is_null()) {
669 for (size_t i = 0; i != a.size1(); ++i) { m(i, i) -= sigma; }
670 } else {
671 m += -sigma * b;
672 }
673 }
674 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
675 if (Bx.is_null()) {
676 y = b * transposed(b) * x;
677 y = inverse(m) * y;
678 y = inverse(b) * y;
679 } else {
680 y = inverse(m) * Bx;
681 y = inverse(b) * y;
682 }
683 }
684 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) {
685 y = b * transposed(b) * x;
686 }
687 void end(
688 Vector<T, Vdense_ref> eigenvalue, Matrix<T, Mrectangle_increment_ref> eigenvector) {
689 Matrix<T, Mrectangle> tmp;
690 tmp = eigenvector;
691 eigenvector = transposed(b) * tmp;
692 }
693 const Matrix<T, BTYPE>& b;
694 Matrix<T, typename ATYPE::inverse> m;
695 };
696
697 struct OPNN3_0 : public OP {
698 OPNN3_0(const Matrix<T, ATYPE>& a_, const Matrix<T, BTYPE>& b_) : a(a_), b(b_) {}
699 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
700 if (Bx.is_null()) {
701 y = b * transposed(b) * x;
702 y = inverse(a) * y;
703 y = inverse(b) * y;
704 } else {
705 y = inverse(a) * Bx;
706 y = inverse(b) * y;
707 }
708 }
709 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) {
710 y = b * transposed(b) * x;
711 }
712 void end(
713 Vector<T, Vdense_ref> eigenvalue, Matrix<T, Mrectangle_increment_ref> eigenvector) {
714 Matrix<T, Mrectangle> tmp;
715 tmp = eigenvector;
716 eigenvector = transposed(b) * tmp;
717 }
718 const Matrix<T, ATYPE>& a;
719 const Matrix<T, BTYPE>& b;
720 };
721};
722
723template <typename T, typename ATYPE, typename BTYPE>
724void Matrix<T, Eigenvector<ATYPE, BTYPE> >::resolve(
725 Matrix<T, Mrectangle_increment_ref>& eigenvector) {
726 size_t size_ = a.size1();
727 std::unique_ptr<OP> op;
728 int mode = 0;
729 const char* bmat;
730 if (b.is_null()) {
731 if (size_ == 1) {
732 eigenvector(0, 0) = 1;
733 eigenvalue[0] = a(0, 0);
734 return;
735 }
736 bmat = "I";
737 if (strcmp(which, "RA") == 0) {
738 strcpy(which, "LA");
739 if (sigma == T(0)) {
740 if (type_a == 'P') {
741 mode = 1;
742 op = std::unique_ptr<OP>(new OPS1i(a));
743 } else {
744 UnimplementedError() << THROW;
745 }
746 } else {
747 mode = 3;
748 op = std::unique_ptr<OP>(new OPS3(a, Matrix<T, BTYPE>::null, sigma));
749 }
750 } else if (strcmp(which, "LA") == 0) {
751 strcpy(which, "SA");
752 if (sigma == T(0)) {
753 if (type_a == 'P') {
754 mode = 1;
755 op = std::unique_ptr<OP>(new OPS1i(a));
756 } else {
757 UnimplementedError() << THROW;
758 }
759 } else {
760 mode = 3;
761 op = std::unique_ptr<OP>(new OPS3(a, Matrix<T, BTYPE>::null, sigma));
762 }
763 } else if (strcmp(which, "SM") == 0) {
764 strcpy(which, "LM");
765 if (sigma == T(0)) {
766 if (type_a == 'P') {
767 mode = 1;
768 op = std::unique_ptr<OP>(new OPS1i(a));
769 } else {
770 UnimplementedError() << THROW;
771 }
772 } else {
773 mode = 3;
774 op = std::unique_ptr<OP>(new OPS3(a, Matrix<T, BTYPE>::null, sigma));
775 }
776 } else {
777 mode = 1;
778 op = std::unique_ptr<OP>(new OPS1(a));
779 }
780 } else {
781 if (size_ == 1) {
782 eigenvector(0, 0) = 1.0 / b(0, 0);
783 eigenvalue[0] = a(0, 0) * eigenvector(0, 0);
784 return;
785 }
786 bmat = "G";
787 if (strcmp(which, "RA") == 0) {
788 strcpy(which, "LA");
789 if (sigma == T(0)) {
790 if ((type_b == 'S' || type_b == 'P') && type_a != 'N') {
791 mode = 3;
792 op = std::unique_ptr<OP>(new OPS3_0(a, b));
793 } else if ((type_a == 'S' || type_a == 'P') && type_b != 'N') {
794 mode = 2;
795 op = std::unique_ptr<OP>(new OPS2i(a, b));
796 } else if (type_a == 'N' && type_b == 'N') {
797 mode = 1;
798 bmat = "I";
799 strcpy(which, "LR");
800 op = std::unique_ptr<OP>(new OPNN1(a, b));
801 } else {
802 UnimplementedError() << THROW;
803 }
804 } else {
805 if ((type_b == 'S' || type_b == 'P') && type_a != 'N') {
806 mode = 3;
807 op = std::unique_ptr<OP>(new OPS3(a, b, sigma));
808 } else if ((type_a == 'S' || type_a == 'P') && type_b != 'N') {
809 mode = 4;
810 op = std::unique_ptr<OP>(new OPS4(a, b, sigma));
811 } else if (type_a == 'N' && type_b == 'N') {
812 mode = 1;
813 bmat = "I";
814 strcpy(which, "LR");
815 op = std::unique_ptr<OP>(new OPNN1_shift(a, b, sigma));
816 } else {
817 UnimplementedError() << THROW;
818 }
819 }
820 } else if (strcmp(which, "LA") == 0) {
821 strcpy(which, "SA");
822 if (sigma == T(0)) {
823 if (type_a == 'P') {
824 mode = 2;
825 op = std::unique_ptr<OP>(new OPS2i(a, b));
826 } else if (type_a == 'N' && type_b == 'N') {
827 mode = 1;
828 bmat = "I";
829 strcpy(which, "SR");
830 op = std::unique_ptr<OP>(new OPNN1(a, b));
831 } else {
832 UnimplementedError() << THROW;
833 }
834 } else {
835 if (type_b == 'S' || type_b == 'P') {
836 mode = 3;
837 op = std::unique_ptr<OP>(new OPS3(a, b, sigma));
838 } else if (type_a == 'S' || type_a == 'P') {
839 mode = 4;
840 op = std::unique_ptr<OP>(new OPS4(a, b, sigma));
841 } else if (type_a == 'N' && type_b == 'N') {
842 mode = 1;
843 bmat = "I";
844 strcpy(which, "SR");
845 op = std::unique_ptr<OP>(new OPNN1_shift(a, b, sigma));
846 } else {
847 UnimplementedError() << THROW;
848 }
849 }
850 } else if (strcmp(which, "SM") == 0) {
851 strcpy(which, "LM");
852 if (sigma == T(0)) {
853 if (type_b == 'S' || type_b == 'P') {
854 mode = 3;
855 op = std::unique_ptr<OP>(new OPS3_0(a, b));
856 } else if (type_a == 'S' || type_a == 'P') {
857 UnimplementedError() << THROW;
858 } else if (type_a == 'N' && type_b == 'N') {
859 mode = 1;
860 bmat = "I";
861 op = std::unique_ptr<OP>(new OPNN1(a, b));
862 } else {
863 UnimplementedError() << THROW;
864 }
865 } else {
866 if (type_b == 'S' || type_b == 'P') {
867 mode = 3;
868 op = std::unique_ptr<OP>(new OPS3(a, b, sigma));
869 } else if (type_a == 'S' || type_a == 'P') {
870 mode = 4;
871 op = std::unique_ptr<OP>(new OPS4(a, b, sigma));
872 } else if (type_a == 'N' && type_b == 'N') {
873 mode = 1;
874 bmat = "I";
875 op = std::unique_ptr<OP>(new OPNN1_shift(a, b, sigma));
876 } else {
877 UnimplementedError() << THROW;
878 }
879 }
880 } else {
881 if (type_b == 'P') {
882 mode = 2;
883 op = std::unique_ptr<OP>(new OPS2(a, b));
884 } else {
885 UnimplementedError() << THROW;
886 }
887 }
888 }
889
890 int ido = 0;
891 Vector<T, Vdense> resid(size_);
892 if (ncv == -1) { ncv = std::min(size_, eigenvalue.size() * 2 + 1); }
893 if (int(eigenvalue.size()) >= ncv) {
894 Exception() << "Eigenvalue extraction: Number of required eigenvalues ("
895 << eigenvalue.size() << ") must be < number of DOFs (" << size_
896 << ") of reduced system." << THROW;
897 }
898 std::vector<T> vv(size_ * ncv);
899 int iparam[11] = {1, 0, 300, 0, 0, 0, 0, 0, 0, 0, 0};
900 if (maxit > 0) { iparam[2] = maxit; }
901 iparam[6] = mode;
902 int ipntr[14];
903 std::vector<T> workd(size_ * 3);
904 size_t lworkl = ncv * (3 * ncv + 8);
905 std::vector<T> workl(lworkl);
906 std::vector<T> rwork(ncv);
907 int info;
908 if (starting_eigenvector.is_null()) {
909 info = 0;
910 } else {
911 info = 1;
912 resid = starting_eigenvector;
913 }
914 do {
915 // arpack::set_debug();
916 arpack::aupd(
917 ido, bmat, size_, which, eigenvalue.size(), tol, &resid[0], ncv, &vv[0], size_,
918 iparam, ipntr, &workd[0], &workl[0], lworkl, &rwork[0], info);
919 switch (ido) {
920 case -1:
921 op->op(
922 Vector<T, Vdense_ref>(size_, &workd[ipntr[0] - 1]),
923 Vector<T, Vdense_ref>(size_, &workd[ipntr[1] - 1]),
924 Vector<T, Vdense_constref>(size_, 0));
925 break;
926 case 1:
927 op->op(
928 Vector<T, Vdense_ref>(size_, &workd[ipntr[0] - 1]),
929 Vector<T, Vdense_ref>(size_, &workd[ipntr[1] - 1]),
930 Vector<T, Vdense_constref>(size_, &workd[ipntr[2] - 1]));
931 break;
932 case 2:
933 op->prod_b(
934 Vector<T, Vdense_constref>(size_, &workd[ipntr[0] - 1]),
935 Vector<T, Vdense_ref>(size_, &workd[ipntr[1] - 1]));
936 break;
937 case 99:
938 break;
939 default:
940 UnimplementedError() << THROW;
941 break;
942 }
943 if (info != 0) {
944 logging::Logger& logger = logging::get_logger("solver.eigenvalue.arpack");
945 switch (info) {
946 case 1:
947 logger << logging::warning
948 << "Eigenvalue extraction: Maximum number of iterations reached.\n"
949 << "All possible eigenvalues (" << int(iparam[4]) << ") have been found."
950 << logging::LOGFLUSH;
951 break;
952 case 3:
953 if (size_t(ncv) == size_) {
954 Exception() << "Eigenvalue extraction: No shifts could be applied during a "
955 "cycle\n"
956 << "of the implicitly restarted Arnoldi iteration." << THROW;
957 }
958 logger << logging::info
959 << "Eigenvalue extraction: No shifts could be applied during a cycle\n"
960 << "of the implicitly restarted Arnoldi iteration. Retry by increasing "
961 "NCV."
962 << logging::LOGFLUSH;
963 ncv = std::min(size_, ncv + eigenvalue.size());
964 vv.resize(size_ * ncv);
965 lworkl = ncv * (3 * ncv + 8);
966 workl.resize(lworkl);
967 rwork.resize(ncv);
968 ido = 0;
969 info = 0;
970 break;
971 case -9999:
972 if (size_t(iparam[4]) < eigenvalue.size() + 2) {
973 Exception() << "Eigenvalue extraction: Could not build an Arnoldi "
974 "factorization\n."
975 << "All possible eigenvalues (" << int(iparam[4])
976 << ") have been found." << THROW;
977 }
978 logger << logging::info
979 << "Eigenvalue extraction: Could not build an Arnoldi factorization.\n"
980 "Retry by reducing the NCV"
981 << logging::LOGFLUSH;
982 ncv = std::min(int(size_), iparam[4]);
983 vv.resize(size_ * ncv);
984 lworkl = ncv * (3 * ncv + 8);
985 workl.resize(lworkl);
986 rwork.resize(ncv);
987 ido = 0;
988 info = 0;
989 break;
990 default:
991 Exception() << "Eigenvalue extraction: arpack eigenvalue extraction failed, "
992 "aupd error="
993 << info << ".\nPossible reason: B matrix is 0." << THROW;
994 break;
995 }
996 }
997 } while (ido != 99);
998 std::vector<int> select(ncv);
999 std::vector<T> D(eigenvalue.size() + 1);
1000 arpack::eupd(
1001 eigenvector.is_null() ? 0 : 1, "A", &select[0], &D[0], eigenvector.m, eigenvector.ldm,
1002 sigma, bmat, size_, which, eigenvalue.size(), tol, &resid[0], ncv, &vv[0], size_, iparam,
1003 ipntr, &workd[0], &workl[0], lworkl, &rwork[0], info);
1004 std::copy(D.begin(), D.end() - 1, eigenvalue.v);
1005 if (info != 0) {
1006 Exception() << "Eigenvalue extraction: arpack eigenvalue extraction failed, eupd error="
1007 << info << "." << THROW;
1008 }
1009 op->end(eigenvalue, eigenvector);
1010}
1011
1012template <typename T, typename ATYPE, typename VTYPE1>
1013Matrix<T, Eigenvector<ATYPE, ATYPE> > eigenvector(
1014 const Matrix<T, ATYPE>& a, char type_a, Vector<T, VTYPE1>& eigenvalue,
1015 const char which[2] = "SM", T sigma = 0, int ncv = -1, double tol = -1, int maxit = -1) {
1016 return Matrix<T, Eigenvector<ATYPE, ATYPE> >(
1017 a, type_a, Matrix<T, ATYPE>::null, ' ', eigenvalue, which, sigma, ncv, tol, maxit,
1018 Vector<T, Vdense_constref>::null);
1019}
1020
1021template <typename T, typename ATYPE, typename VTYPE1, typename VTYPE2>
1022Matrix<T, Eigenvector<ATYPE, ATYPE> > eigenvector(
1023 const Matrix<T, ATYPE>& a, char type_a, Vector<T, VTYPE1>& eigenvalue,
1024 const char which[2] = "SM", T sigma = 0, int ncv = -1, double tol = -1, int maxit = -1,
1025 const Vector<T, VTYPE2>& starting_eigenvector = Vector<T, VTYPE2>::null) {
1026 return Matrix<T, Eigenvector<ATYPE, ATYPE> >(
1027 a, type_a, Matrix<T, ATYPE>::null, ' ', eigenvalue, which, sigma, ncv, tol, maxit,
1028 starting_eigenvector.is_null()
1029 ? Vector<T, Vdense_constref>::null
1030 : static_cast<const Vector<T, Vdense_constref>&>(starting_eigenvector));
1031}
1032
1055template <typename T, typename ATYPE, typename BTYPE, typename VTYPE1>
1056Matrix<T, Eigenvector<ATYPE, BTYPE> > eigenvector(
1057 const Matrix<T, ATYPE>& a, char type_a, const Matrix<T, BTYPE>& b, char type_b,
1058 Vector<T, VTYPE1>& eigenvalue, const char which[2] = "SM", T sigma = 0, int ncv = -1,
1059 double tol = -1, int maxit = -1) {
1060 return Matrix<T, Eigenvector<ATYPE, BTYPE> >(
1061 a, type_a, b, type_b, eigenvalue, which, sigma, ncv, tol, maxit,
1062 Vector<T, Vdense_constref>::null);
1063}
1064
1065template <typename T, typename ATYPE, typename BTYPE, typename VTYPE1>
1066Matrix<T, Eigenvector<ATYPE, BTYPE> > eigenvector(
1067 const Matrix<T, ATYPE>& a, char type_a, const std::vector<Matrix<T, BTYPE> >& b, char type_b,
1068 Vector<T, VTYPE1>& eigenvalue, const char which[2] = "SM", T sigma = 0, int ncv = -1,
1069 double tol = -1, int maxit = -1) {
1070 return Matrix<T, Eigenvector<ATYPE, BTYPE> >(
1071 a, type_a, b, type_b, eigenvalue, which, sigma, ncv, tol, maxit,
1072 Vector<T, Vdense_constref>::null);
1073}
1074
1075template <typename T, typename ATYPE, typename BTYPE, typename VTYPE1, typename VTYPE2>
1076Matrix<T, Eigenvector<ATYPE, BTYPE> > eigenvector(
1077 const Matrix<T, ATYPE>& a, char type_a, const Matrix<T, BTYPE>& b, char type_b,
1078 Vector<T, VTYPE1>& eigenvalue, const char which[2] = "SM", T sigma = 0, int ncv = -1,
1079 double tol = -1, int maxit = -1,
1080 const Vector<T, VTYPE2>& starting_eigenvector = Vector<T, VTYPE2>::null) {
1081 return Matrix<T, Eigenvector<ATYPE, BTYPE> >(
1082 a, type_a, b, type_b, eigenvalue, which, sigma, ncv, tol, maxit,
1083 starting_eigenvector.is_null()
1084 ? Vector<T, Vdense_constref>::null
1085 : static_cast<const Vector<T, Vdense_constref>&>(starting_eigenvector));
1086}
1087
1088template <typename ATYPE, typename BTYPE>
1089struct PolyEigenvector {
1090 using base = PolyEigenvector;
1091 using const_base = PolyEigenvector;
1092};
1093
1094template <typename T, typename ATYPE, typename BTYPE>
1095class Matrix<T, PolyEigenvector<ATYPE, BTYPE> > {
1096public:
1097 Matrix(
1098 const Matrix<T, ATYPE>& a_, const char type_a_, const std::vector<Matrix<T, BTYPE> >& b_,
1099 const char type_b_, Vector<T, Vdense_ref> eigenvalue_r,
1100 Vector<T, Vdense_ref> eigenvalue_i, const char which_[2], T sigma_, int ncv_, double tol_,
1101 int maxit_, const Vector<T, Vdense_constref> starting_eigenvector_)
1102 : a(a_),
1103 type_a(type_a_),
1104 b(b_),
1105 type_b(type_b_),
1106 eigenvaluer(eigenvalue_r),
1107 eigenvaluei(eigenvalue_i),
1108 sigma(sigma_),
1109 ncv(ncv_),
1110 tol(tol_),
1111 maxit(maxit_),
1112 starting_eigenvector(starting_eigenvector_) {
1113 for (size_t i = 0; i != b.size(); ++i) {
1114 if (a.size() != b[i].size()) {
1115 Exception() << "The two matrices are not of the same size." << THROW;
1116 }
1117 }
1118 if (!starting_eigenvector.is_null()
1119 && starting_eigenvector.size() != a.size1() * b.size()) {
1120 Exception() << "The starting eigenvalue vector is not of the same size than the matrix."
1121 << THROW;
1122 }
1123 strcpy(which, which_);
1124 if (tol < 0) { tol = 0; }
1125 }
1126
1127 std::pair<size_t, size_t> size() const {
1128 return std::pair<size_t, size_t>(a.size1(), eigenvaluer.size());
1129 }
1130
1131 size_t size1() const { return a.size1(); }
1132
1133 size_t size2() const { return eigenvaluer.size(); }
1134
1135 void resolve(Matrix<T, Mrectangle_increment_ref>& eigenvector);
1136
1137private:
1138 const Matrix<T, ATYPE>& a;
1139 const char type_a;
1140 const std::vector<Matrix<T, BTYPE> >& b;
1141 const char type_b;
1142 Vector<T, Vdense_ref> eigenvaluer;
1143 Vector<T, Vdense_ref> eigenvaluei;
1144 char which[3];
1145 T sigma;
1146 int ncv;
1147 double tol;
1148 int maxit;
1149 Vector<T, Vdense_constref> starting_eigenvector;
1150
1151 struct OP {
1152 virtual ~OP() {}
1153 virtual void op(
1154 Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {}
1155 virtual void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) {}
1156 virtual void end(
1157 Vector<T, Vdense_ref> eigenvaluer, Vector<T, Vdense_ref> eigenvaluei,
1158 Matrix<T, Mrectangle_increment_ref> eigenvector) {}
1159 };
1160
1161 struct OP1 : public OP {
1162 OP1(const Matrix<T, ATYPE>& a_, const std::vector<Matrix<T, BTYPE> >& b_) : a(a_), b(b_) {}
1163 void op(Vector<T, Vdense_ref> x, Vector<T, Vdense_ref> y, Vector<T, Vdense_constref> Bx) {
1164 size_t s = a.size1();
1165 y(Interval(0, s)) = b[0] * x(Interval(0, s));
1166 for (size_t i = 1; i < b.size(); ++i) {
1167 y(Interval(0, s)) += b[i] * x(Interval(i * s, (i + 1) * s));
1168 }
1169 for (size_t i = 1; i < b.size(); ++i) {
1170 y(Interval(i * s, (i + 1) * s)) = x(Interval((i - 1) * s, i * s));
1171 }
1172 y(Interval(0, s)) = inverse(a) * y(Interval(0, s));
1173 }
1174 void prod_b(Vector<T, Vdense_constref> x, Vector<T, Vdense_ref> y) {
1175 size_t s = a.size1();
1176 y(Interval(0, s)) = a * x(Interval(0, s));
1177 y(Interval(s, y.size())) = x(Interval(s, x.size()));
1178 }
1179 const Matrix<T, ATYPE>& a;
1180 const std::vector<Matrix<T, BTYPE> >& b;
1181 };
1182};
1183
1184template <typename T, typename ATYPE, typename BTYPE>
1185void Matrix<T, PolyEigenvector<ATYPE, BTYPE> >::resolve(
1186 Matrix<T, Mrectangle_increment_ref>& eigenvector) {
1187 size_t size_ = a.size1() * b.size();
1188 if (eigenvaluer.size() > size_ - 2) {
1189 // limitation of arpack
1190 UnimplementedError() << "Eigenvalue extraction: The required number of eigenvalues "
1191 << eigenvaluer.size() << " is too large for the size " << size_
1192 << " of problem." << THROW;
1193 }
1194 std::unique_ptr<OP> op;
1195 int mode = 0;
1196 bool inv = false;
1197 const char* bmat;
1198 bmat = "G";
1199 if (sigma == T(0)) {
1200 if (type_a == 'P') {
1201 mode = 2;
1202 inv = true;
1203 op = std::unique_ptr<OP>(new OP1(a, b));
1204 if (strcmp(which, "SM") == 0) {
1205 strcpy(which, "LM");
1206 } else if (strcmp(which, "RA") == 0) {
1207 strcpy(which, "LR");
1208 } else {
1209 UnimplementedError() << THROW;
1210 }
1211 } else {
1212 UnimplementedError() << THROW;
1213 }
1214 } else {
1215 UnimplementedError() << THROW;
1216 }
1217
1218 int ido = 0;
1219 Vector<T, Vdense> resid(size_);
1220 if (ncv == -1) { ncv = std::min(size_, eigenvaluer.size() * 2 + 1); }
1221 std::vector<T> vv(size_ * (ncv + 1));
1222 int iparam[11] = {1, 0, 300, 1, 0, 0, 0, 0, 0, 0, 0};
1223 if (maxit > 0) { iparam[2] = maxit; }
1224 iparam[6] = mode;
1225 int ipntr[14];
1226 std::vector<T> workd(size_ * 3);
1227 size_t lworkl = ncv * (3 * ncv + 8);
1228 std::vector<T> workl(lworkl);
1229 std::vector<T> rwork(ncv);
1230 int info;
1231 if (starting_eigenvector.is_null()) {
1232 info = 0;
1233 } else {
1234 info = 1;
1235 resid = starting_eigenvector;
1236 }
1237 do {
1238 // arpack::set_debug();
1239 arpack::naupd(
1240 ido, bmat, size_, which, eigenvaluer.size(), tol, &resid[0], ncv, &vv[0], size_,
1241 iparam, ipntr, &workd[0], &workl[0], lworkl, &rwork[0], info);
1242 switch (ido) {
1243 case -1:
1244 op->op(
1245 Vector<T, Vdense_ref>(size_, &workd[ipntr[0] - 1]),
1246 Vector<T, Vdense_ref>(size_, &workd[ipntr[1] - 1]),
1247 Vector<T, Vdense_constref>(size_, 0));
1248 break;
1249 case 1:
1250 op->op(
1251 Vector<T, Vdense_ref>(size_, &workd[ipntr[0] - 1]),
1252 Vector<T, Vdense_ref>(size_, &workd[ipntr[1] - 1]),
1253 Vector<T, Vdense_constref>(size_, &workd[ipntr[2] - 1]));
1254 break;
1255 case 2:
1256 op->prod_b(
1257 Vector<T, Vdense_constref>(size_, &workd[ipntr[0] - 1]),
1258 Vector<T, Vdense_ref>(size_, &workd[ipntr[1] - 1]));
1259 break;
1260 case 99:
1261 break;
1262 default:
1263 UnimplementedError() << THROW;
1264 break;
1265 }
1266 if (info != 0) {
1267 logging::Logger& logger = logging::get_logger("solver.eigenvalue.arpack");
1268 switch (info) {
1269 case 1:
1270 logger << logging::warning
1271 << "Eigenvalue extraction: Maximum number of iterations reached\n."
1272 << "All possible eigenvalues (" << int(iparam[4]) << ") have been found."
1273 << logging::LOGFLUSH;
1274 break;
1275 case 3:
1276 if (size_t(ncv) == size_) {
1277 Exception() << "Eigenvalue extraction: No shifts could be applied during a "
1278 "cycle\n"
1279 << "of the implicitly restarted Arnoldi iteration." << THROW;
1280 }
1281 logger << logging::info
1282 << "Eigenvalue extraction: No shifts could be applied during a cycle of "
1283 "the implicitly\n"
1284 << "restarted Arnoldi iteration. Retry by increasing NCV."
1285 << logging::LOGFLUSH;
1286 ncv = std::min(size_, ncv + eigenvaluer.size());
1287 vv.resize(size_ * ncv);
1288 lworkl = ncv * (3 * ncv + 8);
1289 workl.resize(lworkl);
1290 rwork.resize(ncv);
1291 ido = 0;
1292 info = 0;
1293 break;
1294 case -9999:
1295 if (size_t(iparam[4]) < eigenvaluer.size() + 2) {
1296 Exception() << "Eigenvalue extraction: Could not build an Arnoldi "
1297 "factorization.\n"
1298 << "All possible eigenvalues (" << int(iparam[4])
1299 << ") have been found." << THROW;
1300 }
1301 logger << logging::info
1302 << "Eigenvalue extraction: Could not build an Arnoldi factorization.\n"
1303 " Retry by reducing the NCV"
1304 << logging::LOGFLUSH;
1305 ncv = std::min(int(size_), iparam[4]);
1306 vv.resize(size_ * ncv);
1307 lworkl = ncv * (3 * ncv + 8);
1308 workl.resize(lworkl);
1309 rwork.resize(ncv);
1310 ido = 0;
1311 info = 0;
1312 break;
1313 default:
1314 Exception() << "Eigenvalue extraction: arpack subroutine aupd error=" << info
1315 << ".\nPossible reason: B matrix is 0." << THROW;
1316 break;
1317 }
1318 }
1319 } while (ido != 99);
1320 std::vector<int> select(ncv);
1321 std::vector<T> Dr(eigenvaluer.size() + 1);
1322 std::vector<T> Di(eigenvaluei.size() + 1);
1323 std::vector<T> workev(3 * ncv);
1324 Matrix<T, Mrectangle> eigenvector_tmp(size_, eigenvaluer.size() + 1);
1325 arpack::eupd(
1326 eigenvector.is_null() ? 0 : 1, "A", &select[0], &Dr[0], &Di[0], &eigenvector_tmp.m[0],
1327 size_, sigma, 0, &workev[0], bmat, size_, which, eigenvaluer.size(), tol, &resid[0], ncv,
1328 &vv[0], size_, iparam, ipntr, &workd[0], &workl[0], lworkl, &rwork[0], info);
1329
1330 eigenvector =
1331 eigenvector_tmp(Interval(0, eigenvector.size1()), Interval(0, eigenvector.size2()));
1332 if (inv) {
1333 for (size_t i = 0; i != eigenvaluer.size(); ++i) {
1334 const double d = 1 / (Dr[i] * Dr[i] + Di[i] * Di[i]);
1335 eigenvaluer[i] = Dr[i] * d;
1336 eigenvaluei[i] = -Di[i] * d;
1337 }
1338 } else {
1339 std::copy(Dr.begin(), Dr.end() - 1, eigenvaluer.v);
1340 std::copy(Di.begin(), Di.end() - 1, eigenvaluei.v);
1341 }
1342 if (info != 0) {
1343 Exception() << "Eigenvalue extraction: Error with the arpack subroutine eupd, info = "
1344 << info << "." << THROW;
1345 }
1346 op->end(eigenvaluer, eigenvaluei, eigenvector);
1347}
1348
1371template <typename T, typename ATYPE, typename BTYPE, typename VTYPE1>
1372Matrix<T, PolyEigenvector<ATYPE, BTYPE> > eigenvector(
1373 const Matrix<T, ATYPE>& a, char type_a, const std::vector<Matrix<T, BTYPE> >& b, char type_b,
1374 Vector<T, VTYPE1>& eigenvalue_r, Vector<T, VTYPE1>& eigenvalue_i, const char which[2] = "SM",
1375 T sigma = 0, int ncv = -1, double tol = -1, int maxit = -1) {
1376 return Matrix<T, PolyEigenvector<ATYPE, BTYPE> >(
1377 a, type_a, b, type_b, eigenvalue_r, eigenvalue_i, which, sigma, ncv, tol, maxit,
1378 Vector<T, Vdense_constref>::null);
1379}
1380
1381template <typename T>
1382void copy_real(const Matrix<b2000::csda<T>, Mrectangle_constref>& in, Matrix<T, Mrectangle>& out) {
1383 out.resize(in.size());
1384 const b2000::csda<T>* pin = &in(0, 0);
1385 T* pout = &out(0, 0);
1386 const T* pout_end = pout + out.size1() * out.size2();
1387 for (; pout != pout_end; ++pin, ++pout) { *pout = real(*pin); }
1388}
1389
1390template <typename T>
1391void copy_real(const Matrix<std::complex<T>, Mrectangle_constref>& in, Matrix<T, Mrectangle>& out) {
1392 out.resize(in.size());
1393 const std::complex<T>* pin = &in(0, 0);
1394 T* pout = &out(0, 0);
1395 const T* pout_end = pout + out.size1() * out.size2();
1396 for (; pout != pout_end; ++pin, ++pout) { *pout = real(*pin); }
1397}
1398
1399template <typename T>
1400void copy_imag(const Matrix<b2000::csda<T>, Mrectangle_constref>& in, Matrix<T, Mrectangle>& out) {
1401 out.resize(in.size());
1402 const b2000::csda<T>* pin = &in(0, 0);
1403 T* pout = &out(0, 0);
1404 const T* pout_end = pout + out.size1() * out.size2();
1405 for (; pout != pout_end; ++pin, ++pout) { *pout = imag(*pin); }
1406}
1407
1408template <typename T>
1409void copy_imag(const Matrix<std::complex<T>, Mrectangle_constref>& in, Matrix<T, Mrectangle>& out) {
1410 out.resize(in.size());
1411 const std::complex<T>* pin = &in(0, 0);
1412 T* pout = &out(0, 0);
1413 const T* pout_end = pout + out.size1() * out.size2();
1414 for (; pout != pout_end; ++pin, ++pout) { *pout = imag(*pin); }
1415}
1416
1417template <typename T>
1418void copy(const Matrix<b2000::csda<T>, Mrectangle_constref>& in, Matrix<T, Mrectangle>& out) {
1419 out.resize(in.size());
1420 const b2000::csda<T>* pin = &in(0, 0);
1421 T* pout = &out(0, 0);
1422 const T* pout_end = pout + out.size1() * out.size2();
1423 for (; pout != pout_end; ++pin, ++pout) { *pout = real(*pin); }
1424}
1425
1426template <typename T>
1427void copy(const Matrix<std::complex<T>, Mrectangle_constref>& in, Matrix<T, Mrectangle>& out) {
1428 out.resize(in.size());
1429 const std::complex<T>* pin = &in(0, 0);
1430 T* pout = &out(0, 0);
1431 const T* pout_end = pout + out.size1() * out.size2();
1432 for (; pout != pout_end; ++pin, ++pout) { *pout = real(*pin); }
1433}
1434
1435template <typename T>
1436void copy(const Matrix<T, Mrectangle_constref>& in, Matrix<b2000::csda<T>, Mrectangle>& out) {
1437 out.resize(in.size());
1438 const T* pin = &in(0, 0);
1439 b2000::csda<T>* pout = &out(0, 0);
1440 const b2000::csda<T>* pout_end = pout + out.size1() * out.size2();
1441 for (; pout != pout_end; ++pin, ++pout) { *pout = b2000::csda<T>(*pin); }
1442}
1443
1444template <typename T>
1445void copy(const Matrix<T, Mrectangle_constref>& in, Matrix<std::complex<T>, Mrectangle>& out) {
1446 out.resize(in.size());
1447 const T* pin = &in(0, 0);
1448 std::complex<T>* pout = &out(0, 0);
1449 const std::complex<T>* pout_end = pout + out.size1() * out.size2();
1450 for (; pout != pout_end; ++pin, ++pout) { *pout = std::complex<T>(*pin); }
1451}
1452
1453template <typename T>
1454void copy(const Matrix<T, Mrectangle_constref>& in, Matrix<T, Mrectangle>& out) {
1455 out = in;
1456}
1457
1458} // namespace b2000::b2linalg
1459
1460#endif
THROW
#define THROW
Definition b2exception.H:198
b2000::logging::get_logger
Logger & get_logger(const std::string &logger_name="")
Definition b2logging.H:829
b2000::UnimplementedError
GenericException< UnimplementedError_name > UnimplementedError
Definition b2exception.H:314