eigen_solver.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Non-inline functions for the eigensolvers
27
28#ifdef OOMPH_HAS_MPI
29#include "mpi.h"
30#endif
31
32// Include cfortran.h and the header for the FORTRAN ARPACK routines
33#include "cfortran.h"
34#include "arpack.h"
35#include "lapack_qz.h"
36
37// Oomph-lib headers
38#include "eigen_solver.h"
39#include "linear_solver.h"
40#include "problem.h"
41
42
43namespace oomph
44{
45 //===============================================================
46 /// Constructor, set default values and set the initial
47 /// linear solver to be superlu
48 //===============================================================
50 : EigenSolver(),
51 Spectrum(1),
52 NArnoldi(30),
53 Small(true),
54 Compute_eigenvectors(true)
55 {
57 }
58
59 //===============================================================
60 /// Destructor, delete the default linear solver
61 //===============================================================
63 {
65 }
66
67
68 //==========================================================================
69 /// Use ARPACK to solve an eigen problem that is assembled by elements in
70 /// a mesh in a Problem object.
71 //==========================================================================
72 void ARPACK::solve_eigenproblem(Problem* const& problem_pt,
73 const int& n_eval,
74 Vector<std::complex<double>>& eigenvalue,
75 Vector<DoubleVector>& eigenvector)
76 {
77 bool Verbose = false;
78
79 // Control parameters
80 int ido = 0; // Reverse communication flag
81 char bmat[2]; // Specifies the type of matrix on the RHS
82 // Must be 'I' for standard eigenvalue problem
83 // 'G' for generalised eigenvalue problem
84 int n; // Dimension of the problem
85 char which[3]; // Set which eigenvalues are required.
86 int nev; // Number of eigenvalues desired
87 double tol = 0.0; // Stopping criteria
88 int ncv; // Number of columns of the matrix V (Dimension of Arnoldi
89 // subspace)
90 int iparam[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; // Control parameters
91 // Pointers to starting locations in the workspace vectors
92 int ipntr[14] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
93 int info =
94 0; // Setting to zero gives random initial guess for vector to start
95 // Arnoldi iteration
96
97 // Set up the sizes of the matrix
98 n = problem_pt->ndof(); // Total size of matrix
99 nev = n_eval; // Number of desired eigenvalues
100 ncv = NArnoldi < n ? NArnoldi : n; // Number of Arnoldi vectors allowed
101 // Maximum possible value is max
102 // dimension of matrix
103
104 // If we don't have enough Arnoldi vectors to compute the desired number
105 // of eigenvalues then complain
106 if (nev > ncv)
107 {
108 std::ostringstream warning_stream;
109 warning_stream << "Number of requested eigenvalues " << nev << "\n"
110 << "is greater than the number of Arnoldi vectors:" << ncv
111 << "\n";
112 // Increase the number of Arnoldi vectors
113 ncv = nev + 10;
114 if (ncv > n)
115 {
116 ncv = n;
117 }
118
119 warning_stream << "Increasing number of Arnoldi vectors to " << ncv
120 << "\n but you may want to increase further using\n"
121 << "ARPACK::narnoldi()\n"
122 << "which will also get rid of this warning.\n";
123
125 warning_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
126 }
127
128 // Allocate some workspace
129 int lworkl = 3 * ncv * ncv + 6 * ncv;
130
131 // Array that holds the final set of Arnoldi basis vectors
132 double** v = new double*;
133 *v = new double[ncv * n];
134 // Leading dimension of v (n)
135 int ldv = n;
136
137 // Residual vector
138 double* resid = new double[n];
139 // Workspace vector
140 double* workd = new double[3 * n];
141 // Workspace vector
142 double* workl = new double[lworkl];
143
144 bmat[0] = 'G';
145 // If we require eigenvalues to the left of the shift
146 if (Small)
147 {
148 which[0] = 'S';
149 }
150 // Otherwise we require eigenvalues to the right of the shift
151 else
152 {
153 which[0] = 'L';
154 }
155
156 // Which part of the eigenvalue interests us
157 switch (Spectrum)
158 {
159 // Imaginary part
160 case -1:
161 which[1] = 'I';
162 break;
163
164 // Magnitude
165 case 0:
166 which[1] = 'M';
167 break;
168
169 // Real part
170 case 1:
171 which[1] = 'R';
172 break;
173
174 default:
175 std::ostringstream error_stream;
176 error_stream
177 << "Spectrum is set to an invalid value. It must be 0, 1 or -1\n"
178 << "Not " << Spectrum << std::endl;
179
180 throw OomphLibError(
181 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
182 }
183
184 // %---------------------------------------------------%
185 // | This program uses exact shifts with respect to |
186 // | the current Hessenberg matrix (IPARAM(1) = 1). |
187 // | IPARAM(3) specifies the maximum number of Arnoldi |
188 // | iterations allowed. Mode 3 of DNAUPD is used |
189 // | (IPARAM(7) = 3). All these options can be changed |
190 // | by the user. For details see the documentation in |
191 // | DNAUPD. |
192 // %---------------------------------------------------%
193
194 int ishifts = 1;
195 int maxitr = 300;
196 int mode = 3; // M is symetric and semi-definite
197
198 iparam[0] = ishifts; // Exact shifts wrt Hessenberg matrix H
199 iparam[2] = maxitr; // Maximum number of allowed iterations
200 iparam[3] = 1; // Bloacksize to be used in the recurrence
201 iparam[6] = mode; // Mode is shift and invert in real arithmetic
202
203 // Real and imaginary parts of the shift
204 double sigmar = Sigma_real, sigmai = 0.0;
205
206 // TEMPORARY
207 // only use non distributed matrice and vectors
208 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n, false);
209 this->build_distribution(dist);
210
211 // Before we get into the Arnoldi loop solve the shifted matrix problem
212 // Allocated Row compressed matrices for the mass matrix and shifted main
213 // matrix
214 CRDoubleMatrix M(this->distribution_pt()), AsigmaM(this->distribution_pt());
215
216 // Assemble the matrices
217 problem_pt->get_eigenproblem_matrices(M, AsigmaM, sigmar);
218
219 // Allocate storage for the vectors to be used in matrix vector products
220 DoubleVector rhs(this->distribution_pt(), 0.0);
221 DoubleVector x(this->distribution_pt(), 0.0);
222
223
224 bool LOOP_FLAG = true;
225 bool First = true;
226
227 // Enable resolves for the linear solver
229
230 // Do not report the time taken
232
233 do
234 {
235 // %---------------------------------------------%
236 // | Repeatedly call the routine DNAUPD and take |
237 // | actions indicated by parameter IDO until |
238 // | either convergence is indicated or maxitr |
239 // | has been exceeded. |
240 // %---------------------------------------------%
241 //
242
243 DNAUPD(ido,
244 bmat,
245 n,
246 which,
247 nev,
248 tol,
249 resid,
250 ncv,
251 v,
252 ldv,
253 iparam,
254 ipntr,
255 workd,
256 workl,
257 lworkl,
258 info);
259
260 if (ido == -1)
261 {
262 // %-------------------------------------------%
263 // | Perform matrix vector multiplication |
264 // | y <--- OP*x |
265 // | The user should supply his/her own |
266 // | matrix vector multiplication routine here |
267 // | that takes workd(ipntr(1)) as the input |
268 // | vector, and return the matrix vector |
269 // | product to workd(ipntr(2)). |
270 // %-------------------------------------------%
271 //
272
273 // Firstly multiply by mx
274 for (int i = 0; i < n; i++)
275 {
276 x[i] = workd[ipntr[0] - 1 + i];
277 }
278 M.multiply(x, rhs);
279
280 // Now solve the system
281 DoubleVector temp(rhs);
282 if (First)
283 {
284 Linear_solver_pt->solve(&AsigmaM, temp, rhs);
285 First = false;
286 }
287 else
288 {
289 Linear_solver_pt->resolve(temp, rhs);
290 }
291 temp.clear();
292 // AsigmaM.lubksub(rhs);
293 // Put the solution into the workspace
294 for (int i = 0; i < n; i++)
295 {
296 workd[ipntr[1] - 1 + i] = rhs[i];
297 }
298
299 continue;
300 }
301 else if (ido == 1)
302 {
303 // Already done multiplication by Mx
304 // Need to load the rhs vector
305 for (int i = 0; i < n; i++)
306 {
307 rhs[i] = workd[ipntr[2] - 1 + i];
308 }
309 // Now solve the system
310 // AsigmaM.lubksub(rhs);
311 DoubleVector temp(rhs);
312 if (First)
313 {
314 Linear_solver_pt->solve(&AsigmaM, temp, rhs);
315 First = false;
316 }
317 else
318 {
319 Linear_solver_pt->resolve(temp, rhs);
320 }
321 // Put the solution into the workspace
322 for (int i = 0; i < n; i++)
323 {
324 workd[ipntr[1] - 1 + i] = rhs[i];
325 }
326 continue;
327 }
328 else if (ido == 2)
329 {
330 // Need to multiply by mass matrix
331 // Vector<double> x(n);
332 for (int i = 0; i < n; i++)
333 {
334 x[i] = workd[ipntr[0] - 1 + i];
335 }
336 M.multiply(x, rhs);
337
338 for (int i = 0; i < n; i++)
339 {
340 workd[ipntr[1] - 1 + i] = rhs[i];
341 }
342 continue;
343 }
344
345 if (info < 0)
346 {
347 oomph_info << "ERROR" << std::endl;
348 oomph_info << "Error with _naupd, info = '" << info << std::endl;
349 if (info == -7)
350 {
351 oomph_info << "Not enough memory " << std::endl;
352 }
353 }
354
355 // %-------------------------------------------%
356 // | No fatal errors occurred. |
357 // | Post-Process using DNEUPD. |
358 // | |
359 // | Computed eigenvalues may be extracted. |
360 // | |
361 // | Eigenvectors may also be computed now if |
362 // | desired. (indicated by rvec = .true.) |
363 // %-------------------------------------------%
364 //
365 LOOP_FLAG = false;
366 } while (LOOP_FLAG);
367
368 // Report any problems
369 if (info == 1)
370 {
371 oomph_info << "Maximum number of iterations reached." << std::endl;
372 }
373 else if (info == 3)
374 {
375 oomph_info << "No shifts could be applied during implicit Arnoldi "
376 << "update, try increasing NCV. " << std::endl;
377 }
378
379 // Disable resolves for the linear solver
381
382 // Compute Ritz or Schur vectors, if desired
383 int rvec = Compute_eigenvectors;
384 // Specify the form of the basis computed
385 char howmany[2];
386 howmany[0] = 'A';
387 // Find out the number of converged eigenvalues
388 int nconv = iparam[4];
389
390 // Note that there is an error (feature) in ARPACK that
391 // means in certain cases, if we request eigenvectors,
392 // consequent Schur factorization of the matrix spanned
393 // by the Arnoldi vectors leads to more converged eigenvalues
394 // than reported by DNAPUD. This is a pain because it
395 // causes a segmentation fault and in every case that I've found
396 // the eigenvalues/vectors are not those desired.
397 // At the moment, I'm just going to leave it, but I note
398 // the problem here to remind myself about it.
399
400 // Array to select which Ritz vectors are computed
401 int select[ncv];
402 // Array that holds the real part of the eigenvalues
403 double* eval_real = new double[nconv + 1];
404 // Array that holds the imaginary part of the eigenvalues
405 double* eval_imag = new double[nconv + 1];
406
407 // Workspace
408 double* workev = new double[3 * ncv];
409 // Array to hold the eigenvectors
410 double** z = new double*;
411 *z = new double[(nconv + 1) * n];
412
413 // Error flag
414 int ierr;
415
416 DNEUPD(rvec,
417 howmany,
418 select,
419 eval_real,
420 eval_imag,
421 z,
422 ldv,
423 sigmar,
424 sigmai,
425 workev,
426 bmat,
427 n,
428 which,
429 nev,
430 tol,
431 resid,
432 ncv,
433 v,
434 ldv,
435 iparam,
436 ipntr,
437 workd,
438 workl,
439 lworkl,
440 ierr);
441 //
442 // %-----------------------------------------------%
443 // | The real part of the eigenvalue is returned |
444 // | in the first column of the two dimensional |
445 // | array D, and the imaginary part is returned |
446 // | in the second column of D. The corresponding |
447 // | eigenvectors are returned in the first NEV |
448 // | columns of the two dimensional array V if |
449 // | requested. Otherwise, an orthogonal basis |
450 // | for the invariant subspace corresponding to |
451 // | the eigenvalues in D is returned in V. |
452 // %-----------------------------------------------%
453
454 if (ierr != 0)
455 {
456 oomph_info << "Error with _neupd, info = " << ierr << std::endl;
457 }
458 else
459 // Print the output anyway
460 {
461 // Resize the eigenvalue output vector
462 eigenvalue.resize(nconv);
463 for (int j = 0; j < nconv; j++)
464 {
465 // Turn the output from ARPACK into a complex number
466 std::complex<double> eigen(eval_real[j], eval_imag[j]);
467 // Add the eigenvalues to the output vector
468 eigenvalue[j] = eigen;
469 }
470
472 {
473 // Load the eigenvectors
474 eigenvector.resize(nconv);
475 for (int j = 0; j < nconv; j++)
476 {
477 eigenvector[j].build(this->distribution_pt(), 0.0);
478 for (int i = 0; i < n; i++)
479 {
480 eigenvector[j][i] = z[0][j * n + i];
481 }
482 }
483 }
484 else
485 {
486 eigenvector.resize(0);
487 }
488
489 // Report the information
490 if (Verbose)
491 {
492 oomph_info << " ARPACK " << std::endl;
493 oomph_info << " ====== " << std::endl;
494 oomph_info << " Size of the matrix is " << n << std::endl;
495 oomph_info << " The number of Ritz values requested is " << nev
496 << std::endl;
497 oomph_info << " The number of Arnoldi vectors generated (NCV) is "
498 << ncv << std::endl;
499 oomph_info << " What portion of the spectrum: " << which << std::endl;
500 oomph_info << " The number of converged Ritz values is " << nconv
501 << std::endl;
503 << " The number of Implicit Arnoldi update iterations taken is "
504 << iparam[2] << std::endl;
505 oomph_info << " The number of OP*x is " << iparam[8] << std::endl;
506 oomph_info << " The convergence criterion is " << tol << std::endl;
507 }
508 }
509
510 // Clean up the allocated memory
511 delete[] * z;
512 *z = 0;
513 delete z;
514 z = 0;
515 delete[] workev;
516 workev = 0;
517 delete[] eval_imag;
518 eval_imag = 0;
519 delete[] eval_real;
520 eval_real = 0;
521
522 // Clean up the memory
523 delete[] workl;
524 workl = 0;
525 delete[] workd;
526 workd = 0;
527 delete[] resid;
528 resid = 0;
529 delete[] * v;
530 *v = 0;
531 delete v;
532 v = 0;
533 }
534
535
536 //==========================================================================
537 /// Use LAPACK to solve an eigen problem that is assembled by elements in
538 /// a mesh in a Problem object.
539 //==========================================================================
540 void LAPACK_QZ::solve_eigenproblem(Problem* const& problem_pt,
541 const int& n_eval,
542 Vector<std::complex<double>>& eigenvalue,
543 Vector<DoubleVector>& eigenvector)
544 {
545 // Some character identifiers for use in the LAPACK routine
546 // Do not calculate the left eigenvectors
547 char no_eigvecs[2] = "N";
548 // Do caculate the eigenvectors
549 char eigvecs[2] = "V";
550
551 // Get the dimension of the matrix
552 int n = problem_pt->ndof(); // Total size of matrix
553
554 // Initialise a padding integer
555 int padding = 0;
556 // If the dimension of the matrix is even, then pad the arrays to
557 // make the size odd. This somehow sorts out a strange run-time behaviour
558 // identified by Rich Hewitt.
559 if (n % 2 == 0)
560 {
561 padding = 1;
562 }
563
564 // Get the real and imaginary parts of the shift
565 double sigmar = Sigma_real; // sigmai = 0.0;
566
567 // Actual size of matrix that will be allocated
568 int padded_n = n + padding;
569
570 // Storage for the matrices in the packed form required by the LAPACK
571 // routine
572 double* M = new double[padded_n * padded_n];
573 double* A = new double[padded_n * padded_n];
574
575 // TEMPORARY
576 // only use non-distributed matrices and vectors
577 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n, false);
578 this->build_distribution(dist);
579
580 // Enclose in a separate scope so that memory is cleaned after assembly
581 {
582 // Allocated Row compressed matrices for the mass matrix and shifted main
583 // matrix
584 CRDoubleMatrix temp_M(this->distribution_pt()),
585 temp_AsigmaM(this->distribution_pt());
586
587 // Assemble the matrices
588 problem_pt->get_eigenproblem_matrices(temp_M, temp_AsigmaM, sigmar);
589
590 // Now convert these matrices into the appropriate packed form
591 unsigned index = 0;
592 for (int i = 0; i < n; ++i)
593 {
594 for (int j = 0; j < n; ++j)
595 {
596 M[index] = temp_M(j, i);
597 A[index] = temp_AsigmaM(j, i);
598 ++index;
599 }
600 // If necessary pad the columns with zeros
601 if (padding)
602 {
603 M[index] = 0.0;
604 A[index] = 0.0;
605 ++index;
606 }
607 }
608 // No need to pad the final row because it is never accessed by the
609 // routine.
610 }
611
612 // Temporary eigenvalue storage
613 double* alpha_r = new double[n];
614 double* alpha_i = new double[n];
615 double* beta = new double[n];
616 // Temporary eigenvector storage
617 double* vec_left = new double[1];
618 double* vec_right = new double[n * n];
619
620 // Workspace for the LAPACK routine
621 std::vector<double> work(1, 0.0);
622 // Info flag for the LAPACK routine
623 int info = 0;
624
625 // Call FORTRAN LAPACK to get the required workspace
626 // Note the use of the padding flag
627 LAPACK_DGGEV(no_eigvecs,
628 eigvecs,
629 n,
630 &A[0],
631 padded_n,
632 &M[0],
633 padded_n,
634 alpha_r,
635 alpha_i,
636 beta,
637 vec_left,
638 1,
639 vec_right,
640 n,
641 &work[0],
642 -1,
643 info);
644
645 // Get the amount of requires workspace
646 int required_workspace = (int)work[0];
647 // If we need it
648 work.resize(required_workspace);
649
650 // call FORTRAN LAPACK again with the optimum workspace
651 LAPACK_DGGEV(no_eigvecs,
652 eigvecs,
653 n,
654 &A[0],
655 padded_n,
656 &M[0],
657 padded_n,
658 alpha_r,
659 alpha_i,
660 beta,
661 vec_left,
662 1,
663 vec_right,
664 n,
665 &work[0],
666 required_workspace,
667 info);
668
669 // Now resize storage for the eigenvalues and eigenvectors
670 // We get them all!
671 eigenvalue.resize(n);
672 eigenvector.resize(n);
673
674 // Now convert the output into our format
675 for (int i = 0; i < n; ++i)
676 {
677 // N.B. This is naive and dangerous according to the documentation
678 // beta could be very small giving over or under flow
679 // Remember to fix the shift back again
680 eigenvalue[i] = std::complex<double>(sigmar + alpha_r[i] / beta[i],
681 alpha_i[i] / beta[i]);
682
683 // Resize the eigenvector storage
684 eigenvector[i].build(this->distribution_pt(), 0.0);
685 // Load up the eigenvector (assume that it's real)
686 for (int k = 0; k < n; ++k)
687 {
688 eigenvector[i][k] = vec_right[i * n + k];
689 }
690 }
691
692 // Delete all allocated storage
693 delete[] vec_right;
694 delete[] vec_left;
695 delete[] beta;
696 delete[] alpha_r;
697 delete[] alpha_i;
698 delete[] A;
699 delete[] M;
700 }
701
702
703 //==========================================================================
704 /// Use LAPACK to solve a complex eigen problem specified by the given
705 /// matrices.
706 //==========================================================================
708 const ComplexMatrixBase& A,
709 const ComplexMatrixBase& M,
710 Vector<std::complex<double>>& eigenvalue,
711 Vector<Vector<std::complex<double>>>& eigenvector)
712 {
713 // Some character identifiers for use in the LAPACK routine
714 // Do not calculate the left eigenvectors
715 char no_eigvecs[2] = "N";
716 // Do caculate the eigenvectors
717 char eigvecs[2] = "V";
718
719 // Get the dimension of the matrix
720 int n = A.nrow(); // Total size of matrix
721
722 // Storage for the matrices in the packed form required by the LAPACK
723 // routine
724 double* M_linear = new double[2 * n * n];
725 double* A_linear = new double[2 * n * n];
726
727 // Now convert the matrices into the appropriate packed form
728 unsigned index = 0;
729 for (int i = 0; i < n; ++i)
730 {
731 for (int j = 0; j < n; ++j)
732 {
733 M_linear[index] = M(j, i).real();
734 A_linear[index] = A(j, i).real();
735 ++index;
736 M_linear[index] = M(j, i).imag();
737 A_linear[index] = A(j, i).imag();
738 ++index;
739 }
740 }
741
742 // Temporary eigenvalue storage
743 double* alpha = new double[2 * n];
744 double* beta = new double[2 * n];
745 // Temporary eigenvector storage
746 double* vec_left = new double[2];
747 double* vec_right = new double[2 * n * n];
748
749 // Workspace for the LAPACK routine
750 std::vector<double> work(2, 0.0);
751 std::vector<double> rwork(8 * n, 0.0);
752
753 // Info flag for the LAPACK routine
754 int info = 0;
755
756 // Call FORTRAN LAPACK to get the required workspace
757 // Note the use of the padding flag
758 LAPACK_ZGGEV(no_eigvecs,
759 eigvecs,
760 n,
761 &A_linear[0],
762 n,
763 &M_linear[0],
764 n,
765 alpha,
766 beta,
767 vec_left,
768 1,
769 vec_right,
770 n,
771 &work[0],
772 -1,
773 &rwork[0],
774 info);
775
776 // Get the amount of requires workspace
777 int required_workspace = (int)work[0];
778 // If we need it
779 work.resize(2 * required_workspace);
780
781 // call FORTRAN LAPACK again with the optimum workspace
782 LAPACK_ZGGEV(no_eigvecs,
783 eigvecs,
784 n,
785 &A_linear[0],
786 n,
787 &M_linear[0],
788 n,
789 alpha,
790 beta,
791 vec_left,
792 1,
793 vec_right,
794 n,
795 &work[0],
796 required_workspace,
797 &rwork[0],
798 info);
799
800 // Now resize storage for the eigenvalues and eigenvectors
801 // We get them all!
802 eigenvalue.resize(n);
803 eigenvector.resize(n);
804
805 // Now convert the output into our format
806 for (int i = 0; i < n; ++i)
807 {
808 // We have temporary complex numbers
809 std::complex<double> num(alpha[2 * i], alpha[2 * i + 1]);
810 std::complex<double> den(beta[2 * i], beta[2 * i + 1]);
811
812 // N.B. This is naive and dangerous according to the documentation
813 // beta could be very small giving over or under flow
814 eigenvalue[i] = num / den;
815
816 // Resize the eigenvector storage
817 eigenvector[i].resize(n);
818 // Load up the eigenvector (assume that it's real)
819 for (int k = 0; k < n; ++k)
820 {
821 eigenvector[i][k] = std::complex<double>(
822 vec_right[2 * i * n + 2 * k], vec_right[2 * i * n + 2 * k + 1]);
823 }
824 }
825
826 // Delete all allocated storage
827 delete[] vec_right;
828 delete[] vec_left;
829 delete[] beta;
830 delete[] alpha;
831 delete[] A_linear;
832 delete[] M_linear;
833 }
834} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
int Spectrum
Integer to set whether the real, imaginary or magnitude is required to be small or large.
Definition: eigen_solver.h:115
int NArnoldi
Number of Arnoldi vectors to compute.
Definition: eigen_solver.h:119
LinearSolver * Linear_solver_pt
Pointer to a linear solver.
Definition: eigen_solver.h:107
virtual ~ARPACK()
Destructor, delete the linear solver.
Definition: eigen_solver.cc:62
bool Small
Boolean to set which part of the spectrum left (default) or right of the shifted value.
Definition: eigen_solver.h:124
void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector)
Solve the eigen problem.
Definition: eigen_solver.cc:72
LinearSolver * Default_linear_solver_pt
Pointer to a default linear solver.
Definition: eigen_solver.h:110
ARPACK()
Constructor.
Definition: eigen_solver.cc:49
bool Compute_eigenvectors
Boolean to indicate whether or not to compute the eigenvectors.
Definition: eigen_solver.h:127
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
Abstract base class for matrices of complex doubles – adds abstract interfaces for solving,...
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void clear()
wipes the DoubleVector
Base class for all EigenProblem solves. This simply defines standard interfaces so that different sol...
Definition: eigen_solver.h:61
double Sigma_real
Double value that represents the real part of the shift in shifted eigensolvers.
Definition: eigen_solver.h:65
void find_eigenvalues(const ComplexMatrixBase &A, const ComplexMatrixBase &M, Vector< std::complex< double > > &eigenvalue, Vector< Vector< std::complex< double > > > &eigenvector)
Find the eigenvalues of a generalised eigenvalue problem specified by .
void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector)
Solve the eigen problem.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
virtual void solve(Problem *const &problem_pt, DoubleVector &result)=0
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
virtual void enable_resolve()
Enable resolve (i.e. store matrix and/or LU decomposition, say) Virtual so it can be overloaded to pe...
virtual void resolve(const DoubleVector &rhs, DoubleVector &result)
Resolve the system defined by the last assembled jacobian and the rhs vector. Solution is returned in...
void disable_doc_time()
Disable documentation of solve times.
virtual void disable_resolve()
Disable resolve (i.e. store matrix and/or LU decomposition, say) This function simply resets an inter...
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
virtual void get_eigenproblem_matrices(CRDoubleMatrix &mass_matrix, CRDoubleMatrix &main_matrix, const double &shift=0.0)
Get the matrices required by a eigensolver. If the shift parameter is non-zero the second matrix will...
Definition: problem.cc:8368
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1674
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1246
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...