//////////////////////////////////////////////////////////////////////////////// // // Copyright (c) 2008 The Regents of the University of California // // This file is part of Qbox // // Qbox is distributed under the terms of the GNU General Public License // as published by the Free Software Foundation, either version 2 of // the License, or (at your option) any later version. // See the file COPYING in the root directory of this distribution // or . // //////////////////////////////////////////////////////////////////////////////// // // testMatrix.C // //////////////////////////////////////////////////////////////////////////////// // // multiply a matrix a(m,k) by b(k,n) to get c(m,n) // using blocks of size (mb,nb) on a process grid (nprow,npcol) // // use: testMatrix input_file [-check] // input_file: // nprow npcol // m_a n_a mb_a nb_a transa // m_b n_b mb_b nb_b transb // m_c n_c mb_c nb_c // #include #include #include #include #include #include #include #include using namespace std; #include "Timer.h" #ifdef USE_MPI #include #endif #include "Context.h" #include "Matrix.h" double aa(int i, int j) { return 1.0/(i+1)+2.0*i/(j+1); } double bb(int i, int j) { return i-j-3; } int main(int argc, char **argv) { int mype; int npes; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD, &npes); MPI_Comm_rank(MPI_COMM_WORLD, &mype); char* infilename = argv[1]; ifstream infile(infilename); assert(argc == 2 || argc == 3); if ( !(argc==2 || argc==3) ) { cout << "use: testMatrix inputfile [-check]" << endl; return 1; } bool tcheck = false; if ( argc == 3 ) { if ( !strcmp(argv[2],"-check") ) tcheck = true; else { cerr << " invalid argv[2]" << endl; MPI_Abort(MPI_COMM_WORLD,2); } } Timer tm; int nprow, npcol; int m_a, n_a, mb_a, nb_a; int m_b, n_b, mb_b, nb_b; int m_c, n_c, mb_c, nb_c; char ta, tb; if(mype == 0) { infile >> nprow >> npcol; cout<<"nprow="< 1.e-8 ) { cout << " error at element (" << i << "," << j << ") " << c[ival] << " " << sum << endl; exit(1); } } cout << " results checked" << endl; } cout << " CPU/Real: " << setw(8) << tm.cpu() << " / " << setw(8) << tm.real(); if ( tm.real() > 0.0 ) { int kmax = ( ta == 'n' ) ? a.n() : a.m(); cout << " MFlops: " << (2.0e-6*m_c*n_c*kmax) / tm.real() << endl; } #if 1 double norma=a.nrm2(); if(mype == 0)cout<<"Norm(a)="< w(c.m()); c.syevd('l',w,z); //c.syevx('l',w,z,1.e-5); if (mype == 0) cout << " done" << endl; tm.stop(); if (mype == 0) cout << "Eigenproblem time: " << tm.real() << endl; // complex eigenvalue problem ComplexMatrix cc(c.context(),m_c,n_c,mb_c,nb_c); for ( int m = 0; m < cc.nblocks(); m++ ) for ( int l = 0; l < cc.mblocks(); l++ ) for ( int y = 0; y < cc.nbs(m); y++ ) for ( int x = 0; x < cc.mbs(l); x++ ) { int i = cc.i(l,x); int j = cc.j(m,y); int iii = x + l*cc.mb(); int jjj = y + m*cc.nb(); int ival = iii + jjj * cc.mloc(); if ( i == j ) cc[ival] = i + 1.e-5*drand48(); else cc[ival] = complex(1.e-5*drand48(), 1.e-5*drand48()); } tm.reset(); tm.start(); if (mype == 0) cout << "Complex Eigenproblem... "; ComplexMatrix zz(cc.context(),cc.n(),cc.n(),cc.nb(),cc.nb()); valarray ww(cc.m()); cc.heev('l',ww,zz); //c.syevx('l',w,z,1.e-5); if (mype == 0) cout << " done" << endl; tm.stop(); if (mype == 0) cout << "Complex Eigenproblem time: " << tm.real() << endl; } // Gram-Schmidt orthogonalization of matrix a for ( int m = 0; m < a.nblocks(); m++ ) for ( int l = 0; l < a.mblocks(); l++ ) for ( int y = 0; y < a.nbs(m); y++ ) for ( int x = 0; x < a.mbs(l); x++ ) { int i = a.i(l,x); int j = a.j(m,y); //double aij = aa(i,j); int iii = x + l*a.mb(); int jjj = y + m*a.nb(); int ival = iii + jjj * a.mloc(); if ( i == j ) a[ival] = i + 1.e-6*drand48(); else a[ival] = 1.e-6*drand48(); } DoubleMatrix s(a.context(),a.n(),a.n(),a.nb(),a.nb()); tm.reset(); tm.start(); s.syrk('l','t',2.0,a,0.0); if (mype == 0) cout << "Gram syrk time: " << tm.real() << endl; tm.reset(); tm.start(); s.syr('l',-1.0,a,0,'r'); tm.stop(); if (mype == 0) cout << "Gram syr time: " << tm.real() << endl; // Cholesky decomposition tm.reset(); tm.start(); s.potrf('l'); // Cholesky decomposition: S = L * L^T tm.stop(); if (mype == 0) cout << "Gram Cholesky time: " << tm.real() << endl; // Triangular solve tm.reset(); tm.start(); // solve triangular system X * L^T = C a.trsm('r','l','t','n',1.0,s); tm.stop(); if (mype == 0) cout << "Gram triangular solve time: " << tm.real() << endl; #endif } MPI_Finalize(); }