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 <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;
}

Comments on Example Code

#include <iostream>

With header flens.cxx all of FLENS gets included.

#include <flens/flens.cxx>

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 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 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);

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

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);
}

Comments on Modifications

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

Include the qd-headers

#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;
}

Comments on Modifications

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