GopherProxy

	tmain.cpp - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	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> > (…
	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.eigenvec…
	50 check(fabs(checksum - (Floattype)msize) < epsmsizemsize);
	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…
	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 < epsmsizemsize*2);
	72 checksum = arma::sum(arma::sum(e_a - Diag.eigenvalues()));
	73 cout << "Eigenvalues identical: \t";
	74 check(checksum < epsmsize2);
	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 }