#include #include #include "vector.h" #include "matrix.h" #include "timer.h" // another frend operator to compute eigenvalues template inline int Eigenvalues(Matrix& H, Vector& En) { static Vector work(4*H.size0()); if (H.size0()!=H.size1()) {std::cerr<<"Can not diagonalize non-square matrix!"< Vector operator*(const Matrix& A, const Vector& b) { assert(A.size1()==b.size()); // sizes match Vector c(A.size0()); for (int i=0; i void TestMatrixVector(int N, int M) { cout<<"*************************************************"< A(N,M); A.random(); cout<<"Matrix A is "< b(M); b.random(); cout<<"Vector b is "< c = A*b; cout<<"Product A*b is "< void TestMatrixMatrix(int N0) { cout<<"*************************************************"< A(N0,N0), B(N0,N0); A.random(); B.random(); cout<<"Matrix A is "< C; // C.Product(A,B); C = A * B; cout<<"Matrix A*B is "< En(A.size0()); Eigenvalues(A, En); cout<<"Eigevalues of A are: "< detA = 1; for (int i=0; i detB = 1; for (int i=0; i detC = 1; for (int i=0; i double fpo_matrix(int N0, T) { return 2.*N0*N0*N0;} template<> double fpo_matrix(int N0, double) { return 2.*N0*N0*N0;} template<> double fpo_matrix(int N0, complex) { return 2.*4*N0*N0*N0;} template void TestSpeed() { cout<<"*************************************************"<(N0*1.5); Matrix A(N0,N0), B(N0,N0); Matrix C; A.random(); B.random(); int Nmax = static_cast(1e8/(N0*N0*N0)); if (Nmax>1e6) Nmax = static_cast(1e6); if (Nmax<1) Nmax = 1; Timer t; for (int i=0; i(N0, T(0.0)); double Mflops = 1e-6/time_per_FPO; cout<<"It takes "<(Mflops)<<" Mflops"< test_type; // typedef double test_type; TestMatrixVector(5, 6); TestMatrixMatrix(5); TestSpeed(); return 0; }