OS: Ubuntu 22.04.2 LTS.
Download eigen from website: https://eigen.tuxfamily.org/index.php?title=Main_Page
Small Eigen Project
Write a file main.cpp which use eigen library.
#include <iostream>
#include <Eigen/Dense>
using Eigen::MatrixXd;
int main()
{
MatrixXd m(2,2);
m(0,0) = 3;
m(1,0) = 2.5;
m(0,1) = -1;
m(1,1) = m(1,0) + m(0,1);
std::cout << m << std::endl;
return 0;
}Build and run it.
➜ useEigen g++ -I /home/stephen/Downloads/eigen-3.4.0/ main.cpp -o main
➜ useEigen ./main
3 -1
2.5 1.5Use Lapack In Eigen In C++ Project
Write test code to use lapack, build and run the project.
#include <iostream>
#include <fstream>
#include "Eigen/src/misc/lapacke.h"
using namespace std;
int main(int argc, char** argv)
{
// check for an argument
if (argc<2){
cout << "Usage: " << argv[0] << " " << " filename" << endl;
return -1;
}
int n,m;
double *data;
// read in a text file that contains a real matrix stored in column major format
// but read it into row major format
ifstream fin(argv[1]);
if (!fin.is_open()){
cout << "Failed to open " << argv[1] << endl;
return -1;
}
fin >> n >> m; // n is the number of rows, m the number of columns
data = new double[n*m];
for (int i=0;i<n;i++){
for (int j=0;j<m;j++){
fin >> data[j*n+i];
}
}
if (fin.fail() || fin.eof()){
cout << "Error while reading " << argv[1] << endl;
return -1;
}
fin.close();
// check that matrix is square
if (n != m){
cout << "Matrix is not square" <<endl;
return -1;
}
// allocate data
char Nchar='N';
double *eigReal=new double[n];
double *eigImag=new double[n];
double *vl,*vr;
int one=1;
int lwork=6*n;
double *work=new double[lwork];
int info;
// calculate eigenvalues using the DGEEV subroutine
dgeev_(&Nchar,&Nchar,&n,data,&n,eigReal,eigImag,
vl,&one,vr,&one,
work,&lwork,&info);
// check for errors
if (info!=0){
cout << "Error: dgeev returned error code " << info << endl;
return -1;
}
// output eigenvalues to stdout
cout << "--- Eigenvalues ---" << endl;
for (int i=0;i<n;i++){
cout << "( " << eigReal[i] << " , " << eigImag[i] << " )\n";
}
cout << endl;
// deallocate
delete [] data;
delete [] eigReal;
delete [] eigImag;
delete [] work;
return 0;
}Create a file matrix.txt
3 3
-1.0 -8.0 0.0
-1.0 1.0 -5.0
3.0 0.0 2.0Run the executable file by command line
g++ -I /home/stephen/Downloads/eigen-3.4.0/ test.cpp -o main -D EIGEN_USE_LAPACKE=1 -llapack
./main matrix.txt
Compute SVD By Lapack In Eigen In C++ Project
Here is a project to compute the SVD for a 3×3 matrix. The interface introduction about dgesdd can be found at https://netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_gad8e0f1c83a78d3d4858eaaa88a1c5ab1.html. You can also find dgesdd function declaration in Eigen/src/misc/lapacke.h in Eigen project.
void LAPACK_dgesdd( char* jobz, lapack_int* m, lapack_int* n, double* a,
lapack_int* lda, double* s, double* u, lapack_int* ldu,
double* vt, lapack_int* ldvt, double* work,
lapack_int* lwork, lapack_int* iwork, lapack_int *info );Here is a cpp file using dgesdd to compute svd for original 3×3 matrix.
#include <stdlib.h>
#include <stdio.h>
#include "Eigen/src/misc/lapacke.h"
/* DGESDD prototype */
// extern void dgesdd( char* jobz, int* m, int* n, double* a,
// int* lda, double* s, double* u, int* ldu, double* vt, int* ldvt,
// double* work, int* lwork, int* iwork, int* info );
/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, int m, int n, double* a, int lda );
/* Parameters */
#define M 3
#define N 3
#define LDA M
#define LDU M
#define LDVT N
/* Main program */
int main() {
/* Locals */
int m = M, n = N, lda = LDA, ldu = LDU, ldvt = LDVT, info, lwork;
double wkopt;
double* work;
/* Local arrays */
/* iwork dimension should be at least 8*min(m,n) */
int iwork[8*N];
double s[N], u[LDU*M], vt[LDVT*N];
double a[LDA*N] = {
8915.40487296, -603.33651003, -107.43002595,
-438.80704991, -6950.77014448, -579.47443937,
11.03958149,-144.92407957,1730.06287165
};
/* Executable statements */
printf( " DGESDD Example Program Results\n" );
/* Query and allocate the optimal workspace */
lwork = -1;
// dgesdd_( "Singular vectors", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt,
// &lwork, iwork, &info );
dgesdd_( "S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt,
&lwork, iwork, &info );
lwork = (int)wkopt;
work = (double*)malloc( lwork*sizeof(double) );
/* Compute SVD */
// dgesdd_( "Singular vectors", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work,
// &lwork, iwork, &info );
dgesdd_( "S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work,
&lwork, iwork, &info );
/* Check for convergence */
if( info > 0 ) {
printf( "The algorithm computing SVD failed to converge.\n" );
exit( 1 );
}
/* Print singular values */
print_matrix( "Singular values", 1, n, s, 1 );
/* Print left singular vectors */
print_matrix( "Left singular vectors (stored columnwise)", m, n, u, ldu );
/* Print right singular vectors */
print_matrix( "Right singular vectors (stored rowwise)", n, n, vt, ldvt );
/* Free workspace */
free( (void*)work );
return 0;
} /* End of DGESDD Example */
/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, double* a, int lda ) {
int i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) printf( " %.8f", a[i+j*lda] );
printf( "\n" );
}
}Build and run it.
g++ -I /home/stephen/Downloads/eigen-3.4.0/ final2.cpp -o svd -D EIGEN_USE_LAPACKE=1 -llapack
./svd
DGESDD Example Program Results
Singular values
8936.65538773 6988.40028468 1736.15736644
Left singular vectors (stored columnwise)
-0.99701678 -0.07692282 0.00635843
0.07612665 -0.99359764 -0.08347732
0.01273903 -0.08274424 0.99648939
Right singular vectors (stored rowwise)
-0.99993861 -0.01108052 0.00000001
-0.01108052 0.99993861 -0.00000082
0.00000000 0.00000082 1.00000000