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 (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.\n";
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 dims
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_time << '\n';
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 }