GopherProxy

	tmain.cpp - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	LICENSE
	---
	tmain.cpp (4726B)
	---
	1 #include <iostream>
	2 #include <fstream>
	3 #include <ctime>
	4 #include <armadillo>
	5 #include "header.h"
	6 #include "qrfunc.h"
	7
	8 int main(int argc, char* argv[])
	9 {
	10 // Namespace declarations
	11 using std::cout;
	12
	13 // Timer variables
	14 double tic1, toc1, tic2, toc2;
	15
	16 // Create 2D matrices from Armadillo templates
	17 Lengthtype width;
	18 if (argc == 1) // If no command line argument is given...
	19 width = 4; // Matrices are 4x4.
	20 else
	21 width = atoi(argv[1]); // Else the specified size
	22
	23 const Lengthtype height = width;
	24 cout << "\nInitializing " << width << " by " << height << " matrices.\…
	25 arma::mat A = arma::randu<arma::mat>(width, height);
	26
	27 // QR decomposition is performed upon initialization of QR class
	28 tic1 = clock(); // Start clock1
	29 QR qr(A);
	30 toc1 = clock(); // Stop clock1
	31
	32 //// QR decomposition check
	33 cout << "\n\033[1;33m--- QR decomposition check ---\033[0m\n";
	34 if (verbose == true) {
	35 // Display values to stdout
	36 qr.A.print("Original A:");
	37 }
	38
	39 // Check QR decomposition
	40 if (verbose == true) {
	41 qr.Q.print("Q, after QR dec:");
	42 qr.Q.t().print("Q^T:");
	43 qr.R.print("R, after QR dec:");
	44 cout << '\n';
	45 }
	46
	47 // Check that matrix is orthogonal
	48 arma::mat QTQ = qr.Q.t()*qr.Q;
	49 Floattype checksum = arma::sum(arma::sum(QTQ));
	50 if (verbose == true) {
	51 QTQ.print("Q^T Q:");
	52 cout << "sum = " << checksum << '\n';
	53 }
	54 cout << "Check: Q^T Q = 1: \t";
	55 check(checksum-(Floattype)height < 1e-12f);
	56
	57 cout << "Check: QR = A by \|\|QR-A\|\| = 0: ";
	58 checksum = sum(sum((qr.Q*qr.R)-qr.A));
	59 check(checksum < 1e-12f);
	60
	61
	62 //// Solving linear equations
	63 cout << "\n\033[1;33m--- Solving linear equations: Ax=b ---\033[0m\n";
	64 cout << "Solving QRx=b.\n";
	65 arma::vec b = arma::randu<arma::vec>(qr.size());
	66
	67 // Perform back-substitution of QR system
	68 if (verbose == true) {
	69 b.print("Vector b:");
	70 }
	71
	72 arma::vec x = qr.backsub(b);
	73
	74 if (verbose == true) {
	75 x.print("Solution, x:");
	76 }
	77
	78 cout << "Check: Ax = b by \|Ax-b\| = 0: ";
	79 checksum = arma::sum(arma::sum(qr.A*x - b));
	80 check(checksum < 1e-12f);
	81
	82 //// Calculating the determinant
	83 cout << "\n\033[1;33m--- Determinant of A ---\033[0m\n";
	84 cout << "\|det(A)\| = " << qr.det() << '\n';
	85
	86 //// Calculating the inverse
	87 cout << "\n\033[1;33m--- Inverse of A ---\033[0m\n";
	88 arma::mat Ainv = qr.inverse();
	89 if (verbose == true)
	90 Ainv.print("A^(-1):");
	91 cout << "Check: (A^(-1))^(-1) = A: ";
	92 QR qrinv(Ainv);
	93 arma::mat Ainvinv = qrinv.inverse();
	94 bool equal = true; // Elementwise comparison of values
	95 for (Lengthtype i=0; i<width; ++i) {
	96 for (Lengthtype j=0; j<height; ++j) {
	97 if (fabs(A(i,j)-Ainvinv(i,j)) > 1e12f) {
	98 equal = false;
	99 cout << "At (" << i << "," << j << ") = "
	100 << "A = " << A(i,j)
	101 << ", (A^(-1))^(-1) = " << Ainvinv(i,j) << '\n';
	102 }
	103 }
	104 }
	105 check(equal);
	106
	107 //// Use the Armadillo built-in QR decomposition
	108 tic2 = clock(); // Start clock2
	109 arma::mat Q, R;
	110 arma::qr(Q, R, A);
	111 toc2 = clock(); // Stop clock2
	112
	113 //// Statistics
	114 // Print resulting time of homegrown function and Armadillo function
	115 cout << "\n\033[1;33m--- Performance comparison ---\033[0m\n";
	116 double t1 = (toc1 - tic1)/(CLOCKS_PER_SEC);
	117 double t2 = (toc2 - tic2)/(CLOCKS_PER_SEC);
	118 cout << "Homegrown implementation: \t" << t1 << " s \n"
	119 << "Armadillo implementation: \t" << t2 << " s \n";
	120
	121 cout << "Benchmarking performance across a range of sizes...\n";
	122
	123 // Write results to file
	124 std::ofstream outstream;
	125 // Open outfile for write
	126 outstream.open("performance.dat");
	127 double homegrown_time, armadillo_time;
	128
	129 // Define sizes
	130 Lengthtype dims[] = {32, 64, 128, 256, 512, 1024, 2048};
	131 Lengthtype ndims = sizeof(dims)/sizeof(int); // Number of entries in d…
	132
	133 // Loop through sizes and record performance
	134 for (Lengthtype i=0; i<ndims; ++i) {
	135 cout << " " << dims[i] << std::flush;
	136 // Generate input matrix
	137 arma::mat An = arma::randu<arma::mat>(dims[i],dims[i]);
	138
	139 // Homegrown implementation
	140 tic1 = clock();
	141 QR qrn = QR(An);
	142 toc1 = clock();
	143
	144 // Armadillo implementation
	145 tic2 = clock();
	146 arma::mat Qn, Rn;
	147 arma::qr(Qn, Rn, An);
	148 toc2 = clock();
	149
	150 // Record time spent
	151 homegrown_time = (toc1 - tic1)/(CLOCKS_PER_SEC);
	152 armadillo_time = (toc2 - tic2)/(CLOCKS_PER_SEC);
	153
	154 // Write time to file in three columns
	155 outstream << dims[i] << '\t' << homegrown_time << '\t' << armadillo_…
	156 }
	157
	158 cout << '\n';
	159 // Close file
	160 outstream.close();
	161
	162 // Return successfully
	163 return 0;
	164 }
	165
	166 void check(const bool statement)
	167 {
	168 using std::cout;
	169 if (statement == true)
	170 cout << "\t\033[0;32mPassed\033[0m\n";
	171 else
	172 cout << "\t\033[1;31mFail!!\033[0m\n";
	173 }