1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // testjade.cpp
16 //
17 // use: testjade nprow npcol m mb
18 //
19 ////////////////////////////////////////////////////////////////////////////////
20 
21 #include <cassert>
22 #include <cstdlib>
23 #include <cmath>
24 #include <iostream>
25 #include <iomanip>
26 #include <string>
27 #include <fstream>
28 #include <valarray>
29 #include <algorithm>
30 using namespace std;
31 
32 #include "Timer.h"
33 
34 #ifdef USE_MPI
35 #include <mpi.h>
36 #endif
37 
38 #include "Context.h"
39 #include "Matrix.h"
40 #include "jade.h"
41 
42 const bool print = false;
43 
44 enum MatrixType { FRANK, SCALED_FRANK, LAPLACIAN, SMALL };
45 const MatrixType mtype = SCALED_FRANK;
46 
aa(int i,int j)47 double aa(int i, int j) { return 1.0/(i+1)+2.0*i/(j+1); }
bb(int i,int j)48 double bb(int i, int j) { return i-j-3; }
49 
frank(int n,int i,int j)50 double frank(int n, int i, int j) { return n - max(i,j); }
aijf(int n,int k,int i,int j)51 double aijf(int n, int k, int i, int j)
52 {
53   switch ( mtype )
54   {
55     case FRANK :
56       if ( k == 0 )
57         return frank(n,i,j);
58       else
59         return frank(n,n-i,n-j);
60 
61     case SCALED_FRANK :
62       if ( k == 0 )
63         return frank(n,i,j)/(n*n);
64       else
65         return frank(n,n-i,n-j)/(n*n);
66 
67     case LAPLACIAN :
68       if ( i == j )
69       {
70         return 2.0 + k;
71       }
72       else if ( (i-j)*(i-j) == 1 )
73       {
74         return 1.0;
75       }
76 
77     case SMALL :
78     {
79       if ( k == 0 )
80       {
81         if ( i == j && i == 0 )
82           return 1.0;
83         else if ( i == j && i > 0 )
84           return 4.0;
85         else
86           return 0.0;
87       }
88       else
89       {
90         if ( i == j && i == 0 )
91           return 1.0;
92         else if ( i == j && i > 0 )
93           return 1.0;
94         else
95           return 0.5;
96       }
97     }
98   }
99   return 0.0;
100 }
101 
main(int argc,char ** argv)102 int main(int argc, char **argv)
103 {
104   // use: testjade nprow npcol n nb
105   const int maxsweep = 30;
106   const double tol = 1.e-8;
107   const int nmat = 2;
108 
109   int mype;
110   int npes;
111 #ifdef USE_MPI
112   MPI_Init(&argc,&argv);
113   MPI_Comm_size(MPI_COMM_WORLD, &npes);
114   MPI_Comm_rank(MPI_COMM_WORLD, &mype);
115 #else
116   npes=1;
117   mype=0;
118 #endif
119 
120   Timer tm;
121   if ( argc != 5 )
122   {
123     cout << "use: testjade nprow npcol n nb" << endl;
124     return 1;
125   }
126   int nprow=atoi(argv[1]);
127   int npcol=atoi(argv[2]);
128   int m_a=atoi(argv[3]);
129   int mb_a=atoi(argv[4]);
130   int n_a = m_a;
131   int nb_a = mb_a;
132   if(mype == 0)
133   {
134     //infile >> nprow >> npcol;
135     cout << "nprow=" << nprow << ", npcol=" << npcol << endl;
136     //infile >> m_a >> n_a >> mb_a >> nb_a >> ta;
137     cout << "nmat=" << nmat << " m_a=" << m_a << ", n_a=" << n_a << endl;
138     cout << "tol=" << tol << endl;
139   }
140 #ifdef USE_MPI
141   MPI_Bcast(&nprow, 1, MPI_INT, 0, MPI_COMM_WORLD);
142   MPI_Bcast(&npcol, 1, MPI_INT, 0, MPI_COMM_WORLD);
143   MPI_Bcast(&m_a, 1, MPI_INT, 0, MPI_COMM_WORLD);
144   MPI_Bcast(&n_a, 1, MPI_INT, 0, MPI_COMM_WORLD);
145   MPI_Bcast(&mb_a, 1, MPI_INT, 0, MPI_COMM_WORLD);
146   MPI_Bcast(&nb_a, 1, MPI_INT, 0, MPI_COMM_WORLD);
147 #endif
148   {
149     Context ctxt(MPI_COMM_WORLD,nprow,npcol);
150 
151     if ( mype == 0 )
152     {
153       cout << " Context " << ctxt.ictxt()
154            << ": " << ctxt.nprow() << "x" << ctxt.npcol() << endl;
155     }
156 
157     vector<DoubleMatrix*> a(nmat);
158     vector<vector<double> > adiag(nmat);
159     for ( int k = 0; k < nmat; k++ )
160     {
161       a[k] = new DoubleMatrix(ctxt,m_a,n_a,mb_a,nb_a);
162       adiag[k].resize(n_a);
163     }
164     DoubleMatrix u(ctxt,m_a,n_a,mb_a,nb_a);
165 
166     if ( mype == 0 )
167     {
168       cout << " m_a x n_a / mb_a x nb_a / ta = "
169            << a[0]->m() << "x" << a[0]->n() << " / "
170            << a[0]->mb() << "x" << a[0]->nb() << endl;
171     }
172 
173     for ( int k = 0; k < nmat; k++ )
174     {
175       for ( int m = 0; m < a[k]->nblocks(); m++ )
176         for ( int l = 0; l < a[k]->mblocks(); l++ )
177           for ( int y = 0; y < a[k]->nbs(m); y++ )
178             for ( int x = 0; x < a[k]->mbs(l); x++ )
179             {
180               int i = a[k]->i(l,x);
181               int j = a[k]->j(m,y);
182               //double aij = a.i(l,x) * 10 + a.j(m,y);
183               double aij = aijf(a[k]->n(),k,i,j);
184               int iii = x + l*a[k]->mb();
185               int jjj = y + m*a[k]->nb();
186               int ival = iii + jjj * a[k]->mloc();
187               (*a[k])[ival] = aij;
188             }
189       //cout << " a[" << k << "]=" << endl;
190       //cout << (*a[k]);
191     }
192 
193     tm.start();
194     int nsweep = jade(maxsweep,tol,a,u,adiag);
195     tm.stop();
196     if ( ctxt.onpe0() )
197     {
198       cout << " m=" << m_a << " mb=" << mb_a
199            << " ctxt: " << ctxt.nprow() << "x" << ctxt.npcol()
200            << " nsweep=" << nsweep << " time: " << tm.real() << endl;
201     }
202 
203     for ( int k = 0; k < nmat; k++ )
204       sort(adiag[k].begin(),adiag[k].end());
205 
206     if ( nmat == 1 && (mtype == FRANK || mtype == SCALED_FRANK) )
207     {
208       vector<double> e_exact(adiag[0].size());
209       for ( int i = 0; i < e_exact.size(); i++ )
210       {
211         e_exact[i] = 1.0 / ( 2.0 * ( 1.0 - cos( ((2*i+1)*M_PI)/(2*n_a+1)) ) );
212         if ( mtype == SCALED_FRANK )
213         {
214           int n = a[0]->n();
215           e_exact[i] /= (n*n);
216         }
217       }
218       sort(e_exact.begin(),e_exact.end());
219 
220       if ( mype == 0 )
221       {
222         for ( int k = 0; k < nmat; k++ )
223         {
224           double asum = 0.0;
225           for ( int i = 0; i < e_exact.size(); i++ )
226           {
227             //cout << ctxt.mype() << ": eig[" << k << "][" << i << "]= "
228             //     << adiag[k][i]
229             //     << "  " << e_exact[i] << endl;
230             asum += fabs(adiag[k][i]-e_exact[i]);
231           }
232           cout << "a[" << k << "] sum of abs eigenvalue errors: "
233                << asum << endl;
234         }
235       }
236 
237     }
238 
239     if ( print )
240     {
241       // a[k] contains AU.
242       // compute the product C = U^T A U
243       DoubleMatrix c(*a[0]);
244       for ( int k = 0; k < nmat; k++ )
245       {
246         c.gemm('t','n',1.0,u,(*a[k]),0.0);
247         if ( mype == 0 )
248         {
249           cout << " a[" << k << "]=" << endl;
250           cout << c;
251         }
252       }
253       cout << " u=" << endl;
254       cout << u;
255     }
256   }
257 
258 #ifdef USE_MPI
259   MPI_Finalize();
260 #endif
261 }
262