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