1 #include <Eigen/Core>
2 #include <iostream>
3 #include "timer.h"
4 #include "ArpackFun.h"
5
6 using Eigen::MatrixXd;
7 using Eigen::VectorXd;
8 using Eigen::MatrixXcd;
9 using Eigen::VectorXcd;
10 using Eigen::Lower;
11 typedef Eigen::Map<VectorXd> MapVec;
12
eigs_sym_F77(MatrixXd & M,VectorXd & init_resid,int k,int m,double & time_used,double & prec_err,int & nops)13 void eigs_sym_F77(MatrixXd &M, VectorXd &init_resid, int k, int m,
14 double &time_used, double &prec_err, int &nops)
15 {
16 double start, end;
17 prec_err = -1.0;
18 start = get_wall_time();
19
20 // Begin ARPACK
21 //
22 // Initial value of ido
23 int ido = 0;
24 // 'I' means standard eigen value problem, A * x = lambda * x
25 char bmat = 'I';
26 // dimension of A (n by n)
27 int n = M.rows();
28 // Specify selection criteria
29 // "LM": largest magnitude
30 char which[3] = {'L', 'M', '\0'};
31 // Number of eigenvalues requested
32 int nev = k;
33 // Precision
34 double tol = 1e-10;
35 // Residual vector
36 double *resid = new double[n]();
37 std::copy(init_resid.data(), init_resid.data() + n, resid);
38 // Number of Ritz values used
39 int ncv = m;
40 // Vector of eigenvalues
41 VectorXd evals(nev);
42 // Matrix of eigenvectors
43 MatrixXd evecs(n, ncv);
44
45 // Store final results of eigenvectors
46 // double *V = new double[n * ncv]();
47 double *V = evecs.data();
48 // Leading dimension of V, required by FORTRAN
49 int ldv = n;
50 // Control parameters
51 int *iparam = new int[11]();
52 iparam[1 - 1] = 1; // ishfts
53 iparam[3 - 1] = 1000; // maxitr
54 iparam[7 - 1] = 1; // mode
55 // Some pointers
56 int *ipntr = new int[11]();
57 /* workd has 3 columns.
58 * ipntr[2] - 1 ==> first column to store B * X,
59 * ipntr[1] - 1 ==> second to store Y,
60 * ipntr[0] - 1 ==> third to store X. */
61 double *workd = new double[3 * n]();
62 int lworkl = ncv * (ncv + 8);
63 double *workl = new double[lworkl]();
64 // Error flag. 0 means random initialization,
65 // otherwise using resid as initial value
66 int info = 1;
67
68 saupd(ido, bmat, n, which,
69 nev, tol, resid,
70 ncv, V, ldv,
71 iparam, ipntr, workd,
72 workl, lworkl, info);
73 // ido == -1 or ido == 1 means more iterations needed
74 while (ido == -1 || ido == 1)
75 {
76 MapVec vec_in(&workd[ipntr[0] - 1], n);
77 MapVec vec_out(&workd[ipntr[1] - 1], n);
78 vec_out.noalias() = M.selfadjointView<Lower>() * vec_in;
79
80 saupd(ido, bmat, n, which,
81 nev, tol, resid,
82 ncv, V, ldv,
83 iparam, ipntr, workd,
84 workl, lworkl, info);
85 }
86
87 // info > 0 means warning, < 0 means error
88 if (info > 0)
89 std::cout << "warnings occured" << std::endl;
90 if (info < 0)
91 {
92 delete[] workl;
93 delete[] workd;
94 delete[] ipntr;
95 delete[] iparam;
96 delete[] resid;
97
98 std::cout << "errors occured" << std::endl;
99 end = get_wall_time();
100 time_used = (end - start) * 1000;
101
102 return;
103 }
104
105 // Retrieve results
106 //
107 // Whether to calculate eigenvectors or not.
108 bool rvec = true;
109 // 'A' means to calculate Ritz vectors
110 // 'P' to calculate Schur vectors
111 char howmny = 'A';
112 // Vector of eigenvalues
113 double *d = evals.data();
114 // Used to store results, will use V instead.
115 double *Z = V;
116 // Leading dimension of Z, required by FORTRAN
117 int ldz = n;
118 // Shift
119 double sigma = 0;
120 // Error information
121 int ierr = 0;
122
123 // Number of converged eigenvalues
124 int nconv = 0;
125 // Number of iterations
126 int niter = 0;
127
128 // Use seupd() to retrieve results
129 seupd(rvec, howmny, d,
130 Z, ldz, sigma, bmat,
131 n, which, nev, tol,
132 resid, ncv, V, ldv,
133 iparam, ipntr, workd, workl,
134 lworkl, ierr);
135
136 // Obtain 'nconv' converged eigenvalues
137 nconv = iparam[5 - 1];
138 // 'niter' number of iterations
139 niter = iparam[9 - 1];
140
141 // Free memory of temp arrays
142 delete[] workl;
143 delete[] workd;
144 delete[] ipntr;
145 delete[] iparam;
146 delete[] resid;
147
148 // ierr < 0 means error
149 if (ierr < 0)
150 {
151 std::cout << "errors occured" << std::endl;
152 end = get_wall_time();
153 time_used = (end - start) * 1000;
154
155 return;
156 }
157
158 /* std::cout << "computed eigenvalues D = \n" << evals.transpose() << std::endl;
159 std::cout << "first 5 rows of computed eigenvectors U = \n" <<
160 evecs.topLeftCorner(5, nconv) << std::endl;
161 std::cout << "nconv = " << nconv << std::endl;
162 std::cout << "nops = " << niter << std::endl; */
163
164 end = get_wall_time();
165 time_used = (end - start) * 1000;
166 MatrixXd err = M * evecs.leftCols(nev) - evecs.leftCols(nev) * evals.asDiagonal();
167 prec_err = err.cwiseAbs().maxCoeff();
168 nops = niter;
169 }
170
eigs_gen_F77(MatrixXd & M,VectorXd & init_resid,int k,int m,double & time_used,double & prec_err,int & nops)171 void eigs_gen_F77(MatrixXd &M, VectorXd &init_resid, int k, int m,
172 double &time_used, double &prec_err, int &nops)
173 {
174 double start, end;
175 prec_err = -1.0;
176 start = get_wall_time();
177
178 // Begin ARPACK
179 //
180 // Initial value of ido
181 int ido = 0;
182 // 'I' means standard eigen value problem, A * x = lambda * x
183 char bmat = 'I';
184 // dimension of A (n by n)
185 int n = M.rows();
186 // Specify selection criteria
187 // "LM": largest magnitude
188 char which[3] = {'L', 'M', '\0'};
189 // Number of eigenvalues requested
190 int nev = k;
191 // Precision
192 double tol = 1e-10;
193 // Residual vector
194 double *resid = new double[n]();
195 std::copy(init_resid.data(), init_resid.data() + n, resid);
196 // Number of Ritz values used
197 int ncv = m;
198 // Vector of eigenvalues
199 VectorXd evals_re(nev + 1);
200 VectorXd evals_im(nev + 1);
201 // Matrix of eigenvectors
202 MatrixXd evecs(n, ncv);
203
204 // Store final results of eigenvectors
205 // double *V = new double[n * ncv]();
206 double *V = evecs.data();
207 // Leading dimension of V, required by FORTRAN
208 int ldv = n;
209 // Control parameters
210 int *iparam = new int[11]();
211 iparam[1 - 1] = 1; // ishfts
212 iparam[3 - 1] = 1000; // maxitr
213 iparam[7 - 1] = 1; // mode
214 // Some pointers
215 int *ipntr = new int[14]();
216 /* workd has 3 columns.
217 * ipntr[2] - 1 ==> first column to store B * X,
218 * ipntr[1] - 1 ==> second to store Y,
219 * ipntr[0] - 1 ==> third to store X. */
220 double *workd = new double[3 * n]();
221 int lworkl = 3 * ncv * ncv + 6 * ncv;
222 double *workl = new double[lworkl]();
223 // Error flag. 0 means random initialization,
224 // otherwise using resid as initial value
225 int info = 1;
226
227 naupd(ido, bmat, n, which,
228 nev, tol, resid,
229 ncv, V, ldv,
230 iparam, ipntr, workd,
231 workl, lworkl, info);
232 // ido == -1 or ido == 1 means more iterations needed
233 while (ido == -1 || ido == 1)
234 {
235 MapVec vec_in(&workd[ipntr[0] - 1], n);
236 MapVec vec_out(&workd[ipntr[1] - 1], n);
237 vec_out.noalias() = M * vec_in;
238
239 naupd(ido, bmat, n, which,
240 nev, tol, resid,
241 ncv, V, ldv,
242 iparam, ipntr, workd,
243 workl, lworkl, info);
244 }
245
246 // info > 0 means warning, < 0 means error
247 if (info > 0)
248 std::cout << "warnings occured" << std::endl;
249 if (info < 0)
250 {
251 delete[] workl;
252 delete[] workd;
253 delete[] ipntr;
254 delete[] iparam;
255 delete[] resid;
256
257 std::cout << "errors occured" << std::endl;
258 end = get_wall_time();
259 time_used = (end - start) * 1000;
260
261 return;
262 }
263
264 // Retrieve results
265 //
266 // Whether to calculate eigenvectors or not.
267 bool rvec = true;
268 // 'A' means to calculate Ritz vectors
269 // 'P' to calculate Schur vectors
270 char howmny = 'A';
271 // Vector of eigenvalues
272 double *dr = evals_re.data();
273 double *di = evals_im.data();
274 // Used to store results, will use V instead.
275 double *Z = V;
276 // Leading dimension of Z, required by FORTRAN
277 int ldz = n;
278 // Shift
279 double sigmar = 0;
280 double sigmai = 0;
281 // Working space
282 double *workv = new double[3 * ncv]();
283 // Error information
284 int ierr = 0;
285
286 // Number of converged eigenvalues
287 int nconv = 0;
288 // Number of iterations
289 int niter = 0;
290
291 // Use seupd() to retrieve results
292 neupd(rvec, howmny, dr, di,
293 Z, ldz, sigmar, sigmai, workv,
294 bmat, n, which, nev, tol,
295 resid, ncv, V, ldv, iparam,
296 ipntr, workd, workl, lworkl, ierr);
297
298 // Obtain 'nconv' converged eigenvalues
299 nconv = iparam[5 - 1];
300 // 'niter' number of iterations
301 niter = iparam[9 - 1];
302
303 // Free memory of temp arrays
304 delete[] workv;
305 delete[] workl;
306 delete[] workd;
307 delete[] ipntr;
308 delete[] iparam;
309 delete[] resid;
310
311 // ierr < 0 means error
312 if (ierr < 0)
313 {
314 std::cout << "errors occured" << std::endl;
315 end = get_wall_time();
316 time_used = (end - start) * 1000;
317
318 return;
319 }
320
321 /* VectorXcd evals(evals_re.size());
322 evals.real() = evals_re;
323 evals.imag() = evals_im;
324 std::cout << "computed eigenvalues D = \n" << evals << std::endl;
325 std::cout << "first 5 rows of computed eigenvectors U = \n" <<
326 evecs.topLeftCorner(5, nconv) << std::endl;
327 std::cout << "nconv = " << nconv << std::endl;
328 std::cout << "nops = " << niter << std::endl; */
329
330 end = get_wall_time();
331 time_used = (end - start) * 1000;
332 nops = niter;
333 }
334