Content |
Eigenvalues (Non-Symmetric Matrix)
In this example we compute the eigenvalues \(w \in \mathbb{C}\), left eigenvectors \(V_L\) and right eigenvectors \(V_R\) of a non-symmetric Matrix \(A\).
We use lapack::ev which is the FLENS port of LAPACK's dgeev.
Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
typedef double T;
int
main()
{
typedef GeMatrix<FullStorage<T, ColMajor> > Matrix;
typedef DenseVector<Array<T> > Vector;
const int n = 5;
Matrix A(n, n), VL(n, n), VR(n, n);
Vector wr(n), wi(n);
A = 2, 3, -1, 0, 2,
-6, -5, 0, 2, -6,
2, -5, 6, -6, 2,
2, 3, -1, 0, 8,
-6, -5, 10, 2, -6;
cerr << "A = " << A << endl;
Vector work;
// int optSize = ev_wsq(true, true, A);
// Vector work(optSize);
lapack::ev(true, true, A, wr, wi, VL, VR, work);
cerr << "wr = " << wr << endl;
cerr << "wi = " << wi << endl;
cerr << "VL = " << VL << endl;
cerr << "VR = " << VR << endl;
}
Comments on Example Code
With header flens.cxx all of FLENS gets included.
using namespace std;
using namespace flens;
typedef double T;
int
main()
{
Define some convenient typedefs for the matrix/vector types of our system of linear equations.
typedef DenseVector<Array<T> > Vector;
const int n = 5;
Matrix A(n, n), VL(n, n), VR(n, n);
Vector wr(n), wi(n);
A = 2, 3, -1, 0, 2,
-6, -5, 0, 2, -6,
2, -5, 6, -6, 2,
2, 3, -1, 0, 8,
-6, -5, 10, 2, -6;
cerr << "A = " << A << endl;
Vector for workspace. If this vector has zero length then lapack::ev will do a worksize query and also resize work.
You also could do a worksize query manually
// Vector work(optSize);
Call lapack::ev to compute eigenvalues \(w = w_r+i w_i\), left eigenvectors \(V_L\) and right eigenvectors \(V_R\).
cerr << "wr = " << wr << endl;
cerr << "wi = " << wi << endl;
cerr << "VL = " << VL << endl;
cerr << "VR = " << VR << endl;
}
Compile
$shell> cd flens/examples $shell> clang++ -std=c++0x -Wall -I../.. -o lapack-geev lapack-geev.cc
Run
$shell> cd flens/examples $shell> ./lapack-geev A = 2 3 -1 0 2 -6 -5 0 2 -6 2 -5 6 -6 2 2 3 -1 0 8 -6 -5 10 2 -6 wr = -11.3213 5.9563 5.9563 -0.676184 -2.91514 wi = 0 2.74603 -2.74603 0 0 VL = -0.491133 -0.414826 -0.114499 0.852469 -0.404822 -0.167186 0.270514 -0.144594 0.485239 -0.765229 0.383388 -0.489829 0.274855 0.167246 -0.329007 0.352486 0.61711 0 0.0979434 -0.304721 -0.677942 0.097845 0.121857 -0.0166152 0.22236 VR = 0.217859 0.0128981 0.0260741 0.589244 0.42297 -0.538399 0.22579 -0.0243359 -0.53668 -0.758293 0.0759671 -0.497059 -0.211583 -0.0361856 -0.217856 0.53876 -0.483482 0.25492 0.602872 0.445542 -0.6055 -0.5975 0 -0.00153586 -0.0109669
Using QD (Double-Double and Quad-Double Package)
The following description is taken from the QD Library page:
The QD Library supports both a double-double datatype (approx. 32 decimal digits) and a quad-double datatype (approx. 64 decimal digits). The computational library is written in C++. Both C++ and Fortran-90 high-level language interfaces are provided to permit one to use convert an existing C++ or Fortran-90 program to use the library with only minor changes to the source code. In most cases only a few type statements and (for Fortran-90 programs) read/write statements need to be changed. PSLQ and numerical quadrature programs are included.
Modified Example Code
#include <qd/qd_real.h>
#include <qd/fpu.h>
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
typedef qd_real T;
int
main()
{
unsigned int old_cw;
fpu_fix_start(&old_cw);
typedef GeMatrix<FullStorage<T, ColMajor> > Matrix;
typedef DenseVector<Array<T> > Vector;
const int n = 5;
Matrix A(n, n), VL(n, n), VR(n, n);
Vector wr(n), wi(n);
A = 2, 3, -1, 0, 2,
-6, -5, 0, 2, -6,
2, -5, 6, -6, 2,
2, 3, -1, 0, 8,
-6, -5, 10, 2, -6;
cerr << "A = " << A << endl;
Vector work;
lapack::ev(true, true, A, wr, wi, VL, VR, work);
cerr << "wr = " << wr << endl;
cerr << "wi = " << wi << endl;
cerr << "VL = " << VL << endl;
cerr << "VR = " << VR << endl;
fpu_fix_end(&old_cw);
}
Comments on Modifications
Only a few modifications were made which will be pointed out in the following:
Include the qd-headers
#include <qd/fpu.h>
For using quad-double precision in this example change the typedef to
For understanding the next code lines we first take a look at the QD Library documentation:
The algorithms in the QD library assume IEEE double precision floating point arithmetic. Since Intel x86 processors have extended (80-bit) floating point registers, the round-to-double flag must be enabled in the control word of the FPU for this library to function properly under x86 processors. The function fpu_fix_start turns on the round-to-double bit in the FPU control word, while fpu_fix_end will restore the original state.
So the first thing we do in main is turning on the correct rounding ...
main()
{
unsigned int old_cw;
fpu_fix_start(&old_cw);
... and at the end restore FPU rounding behavior as mentioned above.
}
Compile and Run
$shell> cd flens/examples $shell> clang++ -Wall -std=c++0x -I../.. -I/opt/local/include/ -c qd-lapack-geev.cc $shell> clang++ -o qd-lapack-getrf qd-lapack-geev.o /opt/local/lib/libqd.a $shell> ./qd-lapack-geev A = 2.000000e+00 3.000000e+00 -1.000000e+00 0.000000e+00 2.000000e+00 -6.000000e+00 -5.000000e+00 0.000000e+00 2.000000e+00 -6.000000e+00 2.000000e+00 -5.000000e+00 6.000000e+00 -6.000000e+00 2.000000e+00 2.000000e+00 3.000000e+00 -1.000000e+00 0.000000e+00 8.000000e+00 -6.000000e+00 -5.000000e+00 1.000000e+01 2.000000e+00 -6.000000e+00 wr = -1.132127e+01 5.956297e+00 5.956297e+00 -6.761839e-01 -2.915142e+00 wi = 0.000000e+00 2.746028e+00 -2.746028e+00 0.000000e+00 0.000000e+00 VL = -4.911328e-01 -4.148256e-01 -1.144990e-01 -8.524688e-01 -4.048219e-01 -1.671857e-01 2.705139e-01 -1.445942e-01 -4.852388e-01 -7.652287e-01 3.833882e-01 -4.898286e-01 2.748545e-01 -1.672458e-01 -3.290069e-01 3.524856e-01 6.171098e-01 0.000000e+00 -9.794343e-02 -3.047207e-01 -6.779417e-01 9.784495e-02 1.218572e-01 1.661517e-02 2.223600e-01 VR = 2.178595e-01 -1.289807e-02 -2.607406e-02 -5.892439e-01 4.229696e-01 -5.383990e-01 -2.257901e-01 2.433586e-02 5.366796e-01 -7.582926e-01 7.596711e-02 4.970590e-01 2.115828e-01 3.618556e-02 -2.178560e-01 5.387600e-01 4.834818e-01 -2.549198e-01 -6.028721e-01 4.455418e-01 -6.055002e-01 5.975004e-01 0.000000e+00 1.535864e-03 -1.096688e-02
Using mpfr::real (C++ Interface for mpfr)
The following description is taken from the mpfr::real page:
The class mpfr::real is a high-level C++ wrapper for the GNU MPFR library, a C library for multiple-precision floating-point computations with correct rounding.
The objects of mpfr::real can (almost) be used like doubles or other fundamental floating-point types, thus avoiding the use of C functions to manipulate MPFR's low-level data type directly. In addition to all arithmetic and comparison operators available for fundamental floating-point types, mpfr::real also supports the set of mathematical functions for double from math.h/cmath. Finally, std::istream operator >> and std::ostream operator << are implemented for mpfr::reals. This allows to substitute double with mpfr::real with no further modifications of the code in many cases.
Modified Example Code
#define REAL_ENABLE_CONVERSION_OPERATORS
#include <external/real.hpp>
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
typedef mpfr::real<53> T;
int
main()
{
typedef GeMatrix<FullStorage<T, ColMajor> > Matrix;
typedef DenseVector<Array<T> > Vector;
const int n = 5;
Matrix A(n, n), VL(n, n), VR(n, n);
Vector wr(n), wi(n);
A = 2, 3, -1, 0, 2,
-6, -5, 0, 2, -6,
2, -5, 6, -6, 2,
2, 3, -1, 0, 8,
-6, -5, 10, 2, -6;
cerr << "A = " << A << endl;
Vector work;
// int optSize = ev_wsq(true, true, A);
// Vector work(optSize);
lapack::ev(true, true, A, wr, wi, VL, VR, work);
cerr << "wr = " << wr << endl;
cerr << "wi = " << wi << endl;
cerr << "VL = " << VL << endl;
cerr << "VR = " << VR << endl;
}
Comments on Modifications
Include header file for mpfr::real before flens.cxx and make sure that conversion operators are enabled
#include <external/real.hpp>
#include <flens/flens.cxx>
Make typedef for using mpfr::real
Vector for workspace. If this vector has zero length then lapack::ev will do a worksize query and also resize work.
You also could do a worksize query manually
// Vector work(optSize);
Call lapack::ev to compute eigenvalues \(w = w_r+i w_i\), left eigenvectors \(V_L\) and right eigenvectors \(V_R\).
Compile and Run
$shell> cd flens/examples $shell> clang++ -Wall -std=c++0x -I../.. -I/opt/local/include -o mpfr-real-lapack-geev mpfr-real-lapack-geev.cc -L /opt/local/lib -lmpfr -lgmp $shell> ./mpfr-real-lapack-geev A = 2.000000e+00 3.000000e+00 -1.000000e+00 0.000000e+00 2.000000e+00 -6.000000e+00 -5.000000e+00 0.000000e+00 2.000000e+00 -6.000000e+00 2.000000e+00 -5.000000e+00 6.000000e+00 -6.000000e+00 2.000000e+00 2.000000e+00 3.000000e+00 -1.000000e+00 0.000000e+00 8.000000e+00 -6.000000e+00 -5.000000e+00 1.000000e+01 2.000000e+00 -6.000000e+00 wr = -1.132127e+01 5.956297e+00 5.956297e+00 -6.761839e-01 -2.915142e+00 wi = 0.000000e+00 2.746028e+00 -2.746028e+00 0.000000e+00 0.000000e+00 VL = -4.911328e-01 -4.148256e-01 -1.144990e-01 8.524688e-01 -4.048219e-01 -1.671857e-01 2.705139e-01 -1.445942e-01 4.852388e-01 -7.652287e-01 3.833882e-01 -4.898286e-01 2.748545e-01 1.672458e-01 -3.290069e-01 3.524856e-01 6.171098e-01 0.000000e+00 9.794343e-02 -3.047207e-01 -6.779417e-01 9.784495e-02 1.218572e-01 -1.661517e-02 2.223600e-01 VR = 2.178595e-01 1.289807e-02 2.607406e-02 5.892439e-01 4.229696e-01 -5.383990e-01 2.257901e-01 -2.433586e-02 -5.366796e-01 -7.582926e-01 7.596711e-02 -4.970590e-01 -2.115828e-01 -3.618556e-02 -2.178560e-01 5.387600e-01 -4.834818e-01 2.549198e-01 6.028721e-01 4.455418e-01 -6.055002e-01 -5.975004e-01 0.000000e+00 -1.535864e-03 -1.096688e-02