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 <iostream>
#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(truetrue, A, wr, wi, VL, VR, work);

cerr << "wr = " << wr << endl;
cerr << "wi = " << wi << endl;
cerr << "VL = " << VL << endl;
cerr << "VR = " << VR << endl;
}

With header flens.cxx all of FLENS gets included.

#include <flens/flens.cxx>

Define some convenient typedefs for the matrix/vector types of our system of linear equations.

typedef GeMatrix<FullStorage<T, ColMajor> >   Matrix;
typedef DenseVector<Array<T> >                Vector;

Vector for workspace. If this vector has zero length then lapack::ev will do a worksize query and also resize work.

Vector   work;

You also could do a worksize query manually

// int     optSize = ev_wsq(true, true, A);
// Vector  work(optSize);

Call lapack::ev to compute eigenvalues $$w = w_r+i w_i$$, left eigenvectors $$V_L$$ and right eigenvectors $$V_R$$.

lapack::ev(truetrue, A, wr, wi, VL, VR, work);

Compile

$shell> cd flens/examples$shell> g++ -std=c++11 -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 <iostream>

#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(truetrue, 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);
}

Only a few modifications were made which will be pointed out in the following:

#include <qd/qd_real.h>
#include <qd/fpu.h>

For using quad-double precision in this example change the typedef to

typedef qd_real   T;

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 ...

int
main()
{
unsigned int old_cw;
fpu_fix_start(&old_cw);

... and at the end restore FPU rounding behavior as mentioned above.

fpu_fix_end(&old_cw);
}

Compile and Run

$shell> cd flens/examples$shell> g++ -Wall -std=c++11 -I../.. -I/opt/local/include/ -c qd-lapack-geev.cc
$shell> g++ -o qd-lapack-geev 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

#include <iostream>

#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(truetrue, A, wr, wi, VL, VR, work);

cerr << "wr = " << wr << endl;
cerr << "wi = " << wi << endl;
cerr << "VL = " << VL << endl;
cerr << "VR = " << VR << endl;
}

Include header file for mpfr::real before flens.cxx and make sure that conversion operators are enabled

#define REAL_ENABLE_CONVERSION_OPERATORS
#include <external/real.hpp>
#include <flens/flens.cxx>

Make typedef for using mpfr::real

typedef mpfr::real<53>   T;

Vector for workspace. If this vector has zero length then lapack::ev will do a worksize query and also resize work.

Vector   work;

You also could do a worksize query manually

// int     optSize = ev_wsq(true, true, A);
// Vector  work(optSize);

Call lapack::ev to compute eigenvalues $$w = w_r+i w_i$$, left eigenvectors $$V_L$$ and right eigenvectors $$V_R$$.

lapack::ev(truetrue, A, wr, wi, VL, VR, work);

Compile and Run

$shell> cd flens/examples$shell> g++ -Wall -std=c++11 -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