tmain.cpp - numeric - C++ library with numerical algorithms
(HTM) git clone git://src.adamsgaard.dk/numeric
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tmain.cpp (2762B)
---
1 #include <iostream>
2 #include <armadillo>
3 #include "header.h"
4 #include "jacobi.h"
5
6 int main(int argc, char* argv[])
7 {
8 // Namespace declarations
9 using std::cout;
10
11 // Define matrix size
12 Lengthtype msize = 6;
13 if (argc == 2) {
14 if ((strcmp(argv[1], "-h") == 0) || (strcmp(argv[1], "--help") == 0)) {
15 cout << "Usage: " << argv[0] << " [matrix size]\n"
16 << "If matrix size is not specified, "
17 << "the matrix width and length will be "
18 << msize << ".\n";
19 return 1;
20 }
21 msize = atoi(argv[1]); // Use specified matrix size
22 }
23
24 // Calculate machine precision
25 Floattype eps = 1.0f;
26 while (1.0f + eps != 1.0f)
27 eps /= 2.0f;
28 //cout << "Machine precision of '" << typeid(eps).name()
29 // << "' type is: eps = " << eps << '\n';
30
31 Floattype checksum;
32 // Generate input matrix A, which is symmetric
33 cout << "\n\033[1;33m--- Input data check ---\033[0m\n";
34 arma::Mat<Floattype> A = symmatu(arma::randu< arma::Mat<Floattype> > (msize, msize));
35 checksum = arma::sum(arma::sum(A - A.t()));
36 cout << "Symmetry check: A = A^T: ";
37 check(checksum < eps);
38 if (verbose == true) {
39 A.print("Original matrix:");
40 }
41
42 // Perform Jacobi diagonalization of matrix A
43 Jacobi Diag = Jacobi(A);
44
45 cout << "\n\033[1;33m--- Diagonalization check ---\033[0m\n";
46 if (verbose == true)
47 Diag.trans().print("Transformed matrix (At):");
48 cout << "Check: V V^T = 1: \t";
49 checksum = arma::sum(arma::sum(Diag.eigenvectors().t() * Diag.eigenvectors()));
50 check(fabs(checksum - (Floattype)msize) < eps*msize*msize);
51 if (verbose == true) {
52 Diag.eigenvalues().print("Eigenvalues (e):");
53 Diag.eigenvectors().print("Eigenvectors (V):");
54 Diag.eigenvectors().t().print("V^T");
55 (Diag.eigenvectors() * Diag.eigenvectors().t()).print("V V^T");
56 (Diag.eigenvectors().t() * A * Diag.eigenvectors()).print("V^T A V");
57 (Diag.eigenvectors().t() * Diag.trans() * Diag.eigenvectors()).print("V^T At V");
58 }
59
60 // Armadillo implementation
61 arma::Mat<Floattype> V_a (msize, msize);
62 arma::Col<Floattype> e_a (msize);
63 arma::eig_sym(e_a, V_a, A);
64 if (verbose == true) {
65 e_a.print("Armadillo eigenvalues:");
66 V_a.print("Armadillo eigenvectors:");
67 }
68 cout << "\n\033[1;33m--- Armadillo comparison ---\033[0m\n";
69 checksum = arma::sum(arma::sum(V_a - Diag.eigenvectors()));
70 cout << "Eigenvectors identical:\t";
71 check(checksum < eps*msize*msize*2);
72 checksum = arma::sum(arma::sum(e_a - Diag.eigenvalues()));
73 cout << "Eigenvalues identical: \t";
74 check(checksum < eps*msize*2);
75
76
77 // Return successfully
78 return 0;
79 }
80
81 void check(const bool statement)
82 {
83 using std::cout;
84 if (statement == true)
85 cout << "\t\033[0;32mPassed\033[0m\n";
86 else
87 cout << "\t\033[1;31mFail!!\033[0m\n";
88 }