//////////////////////////////////////////////////////////////////////////////// // // 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 . // //////////////////////////////////////////////////////////////////////////////// // // testjade.C // // use: testjade nprow npcol m mb // //////////////////////////////////////////////////////////////////////////////// #include #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" #include "jade.h" const bool print = false; enum MatrixType { FRANK, SCALED_FRANK, LAPLACIAN, SMALL }; const MatrixType mtype = SCALED_FRANK; 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; } double frank(int n, int i, int j) { return n - max(i,j); } double aijf(int n, int k, int i, int j) { switch ( mtype ) { case FRANK : if ( k == 0 ) return frank(n,i,j); else return frank(n,n-i,n-j); case SCALED_FRANK : if ( k == 0 ) return frank(n,i,j)/(n*n); else return frank(n,n-i,n-j)/(n*n); case LAPLACIAN : if ( i == j ) { return 2.0 + k; } else if ( (i-j)*(i-j) == 1 ) { return 1.0; } case SMALL : { if ( k == 0 ) { if ( i == j && i == 0 ) return 1.0; else if ( i == j && i > 0 ) return 4.0; else return 0.0; } else { if ( i == j && i == 0 ) return 1.0; else if ( i == j && i > 0 ) return 1.0; else return 0.5; } } } return 0.0; } int main(int argc, char **argv) { // use: testjade nprow npcol n nb const int maxsweep = 30; const double tol = 1.e-8; const int nmat = 2; int mype; int npes; #ifdef USE_MPI MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD, &npes); MPI_Comm_rank(MPI_COMM_WORLD, &mype); #else npes=1; mype=0; #endif Timer tm; assert(argc==5); int nprow=atoi(argv[1]); int npcol=atoi(argv[2]); int m_a=atoi(argv[3]); int mb_a=atoi(argv[4]); int n_a = m_a; int nb_a = mb_a; if(mype == 0) { //infile >> nprow >> npcol; cout << "nprow=" << nprow << ", npcol=" << npcol << endl; //infile >> m_a >> n_a >> mb_a >> nb_a >> ta; cout << "nmat=" << nmat << " m_a=" << m_a << ", n_a=" << n_a << endl; cout << "tol=" << tol << endl; } #ifdef USE_MPI MPI_Bcast(&nprow, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&npcol, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&m_a, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&n_a, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&mb_a, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(&nb_a, 1, MPI_INT, 0, MPI_COMM_WORLD); #endif { Context ctxt(MPI_COMM_WORLD,nprow,npcol); if ( mype == 0 ) { cout << " Context " << ctxt.ictxt() << ": " << ctxt.nprow() << "x" << ctxt.npcol() << endl; } vector a(nmat); vector > adiag(nmat); for ( int k = 0; k < nmat; k++ ) { a[k] = new DoubleMatrix(ctxt,m_a,n_a,mb_a,nb_a); adiag[k].resize(n_a); } DoubleMatrix u(ctxt,m_a,n_a,mb_a,nb_a); if ( mype == 0 ) { cout << " m_a x n_a / mb_a x nb_a / ta = " << a[0]->m() << "x" << a[0]->n() << " / " << a[0]->mb() << "x" << a[0]->nb() << endl; } for ( int k = 0; k < nmat; k++ ) { for ( int m = 0; m < a[k]->nblocks(); m++ ) for ( int l = 0; l < a[k]->mblocks(); l++ ) for ( int y = 0; y < a[k]->nbs(m); y++ ) for ( int x = 0; x < a[k]->mbs(l); x++ ) { int i = a[k]->i(l,x); int j = a[k]->j(m,y); //double aij = a.i(l,x) * 10 + a.j(m,y); double aij = aijf(a[k]->n(),k,i,j); int iii = x + l*a[k]->mb(); int jjj = y + m*a[k]->nb(); int ival = iii + jjj * a[k]->mloc(); (*a[k])[ival] = aij; } //cout << " a[" << k << "]=" << endl; //cout << (*a[k]); } tm.start(); int nsweep = jade(maxsweep,tol,a,u,adiag); tm.stop(); if ( ctxt.onpe0() ) { cout << " m=" << m_a << " mb=" << mb_a << " ctxt: " << ctxt.nprow() << "x" << ctxt.npcol() << " nsweep=" << nsweep << " time: " << tm.real() << endl; } for ( int k = 0; k < nmat; k++ ) sort(adiag[k].begin(),adiag[k].end()); if ( nmat == 1 && (mtype == FRANK || mtype == SCALED_FRANK) ) { vector e_exact(adiag[0].size()); for ( int i = 0; i < e_exact.size(); i++ ) { e_exact[i] = 1.0 / ( 2.0 * ( 1.0 - cos( ((2*i+1)*M_PI)/(2*n_a+1)) ) ); if ( mtype == SCALED_FRANK ) { int n = a[0]->n(); e_exact[i] /= (n*n); } } sort(e_exact.begin(),e_exact.end()); if ( mype == 0 ) { for ( int k = 0; k < nmat; k++ ) { double asum = 0.0; for ( int i = 0; i < e_exact.size(); i++ ) { //cout << ctxt.mype() << ": eig[" << k << "][" << i << "]= " // << adiag[k][i] // << " " << e_exact[i] << endl; asum += fabs(adiag[k][i]-e_exact[i]); } cout << "a[" << k << "] sum of abs eigenvalue errors: " << asum << endl; } } } if ( print ) { // a[k] contains AU. // compute the product C = U^T A U DoubleMatrix c(*a[0]); for ( int k = 0; k < nmat; k++ ) { c.gemm('t','n',1.0,u,(*a[k]),0.0); if ( mype == 0 ) { cout << " a[" << k << "]=" << endl; cout << c; } } cout << " u=" << endl; cout << u; } } #ifdef USE_MPI MPI_Finalize(); #endif }