Feast.cxx
1 #ifndef SELDON_FILE_FEAST_CXX
2 #define SELDON_FILE_FEAST_CXX
3 
4 namespace Seldon
5 {
6 
7  // main function to find eigenvalues and eigenvectors with Feast (MKL implementation)
8 #ifdef SELDON_WITH_VIRTUAL
9  template<class T>
10  void FindEigenvaluesFeast(EigenProblem_Base<T>& var,
11  Vector<T>& eigen_values,
12  Vector<T>& lambda_imag,
13  Matrix<T, General, ColMajor>& eigen_vectors)
14 #else
15  template<class EigenProblem, class T, class Allocator1,
16  class Allocator2, class Allocator3>
17  void FindEigenvaluesFeast(EigenProblem& var,
18  Vector<T, VectFull, Allocator1>& eigen_values,
19  Vector<T, VectFull, Allocator2>& lambda_imag,
20  Matrix<T, General, ColMajor, Allocator3>& eigen_vectors)
21 #endif
22  {
23  // initialization of feast parameters
24  IVect fpm(128);
25  fpm.Fill(0);
26  feastinit_(fpm.GetData());
27 
28  fpm(0) = 0;
29  if (var.GetPrintLevel() > 0)
30  fpm(0) = 1;
31 
32  double tol = var.GetStoppingCriterion();
33  fpm(2) = -log10(tol);
34 
35  for (int i = 64; i < fpm.GetM(); i++)
36  fpm(i) = 0;
37 
38  fpm(3) = var.GetNbMaximumIterations();
39 
40  FeastParam& param = var.GetFeastParameters();
41  if (param.EstimateNumberEigenval())
42  fpm(13) = 2;
43 
44  // regular mode only
45  // no shift is applied, the eigenvalues are directly searched into a circle
46  var.SetComputationalMode(var.REGULAR_MODE);
47 
48 #ifdef SELDON_WITH_MPI
49  MPI_Comm solver_comm; // duplication needed for feast
50  MPI_Comm_dup(var.GetCommunicator(), &solver_comm);
51 
52  fpm(8) = MPI_Comm_c2f(var.GetRciCommunicator());
53  fpm(48) = MPI_Comm_c2f(solver_comm);
54  fpm(58) = MPI_Comm_c2f(var.GetGlobalCommunicator());
55 #endif
56 
57  // initialization of computation
58  int n = var.GetM();
59 
60  typedef typename ClassComplexType<T>::Tcplx Tcplx;
61  typedef typename ClassComplexType<T>::Treal Treal;
62 
63  int m0 = var.GetNbAskedEigenvalues()+1;
64 
65  Treal emin = param.GetLowerBoundInterval();
66  Treal emax = param.GetUpperBoundInterval();
67  Tcplx Emid = param.GetCircleCenterSpectrum();
68  Treal r = param.GetCircleRadiusSpectrum();
69  fpm(17) = param.GetRatioEllipseSpectrum();
70  fpm(18) = param.GetAngleEllipseSpectrum();
71 
72  T zero; Tcplx one;
73  SetComplexZero(zero);
74  SetComplexOne(one);
75 
76  bool herm = var.IsHermitianProblem();
77  bool sym = var.IsSymmetricProblem();
78  if (!herm)
79  {
80  if (var.UseCholeskyFactoForMass())
81  {
82  cout << "Cholesky case not implemented in Feast interface" << endl;
83  abort();
84  }
85  }
86 
87  if (herm || sym)
88  fpm(14) = 2;
89  else
90  fpm(14) = 1;
91 
92  if (herm)
93  {
94  if (param.GetTypeIntegration() >= 0)
95  fpm(15) = param.GetTypeIntegration();
96 
97  if (param.GetNumOfQuadraturePoints() > 0)
98  fpm(1) = param.GetNumOfQuadraturePoints();
99  else
100  fpm(1) = 8;
101  }
102  else
103  {
104  if (param.GetTypeIntegration() >= 0)
105  fpm(15) = param.GetTypeIntegration();
106 
107  if (param.GetNumOfQuadraturePoints() > 0)
108  fpm(7) = param.GetNumOfQuadraturePoints();
109  else
110  fpm(7) = 16;
111  }
112 
113  // work arrays
114  Matrix<T, General, ColMajor> work;
115  Matrix<Tcplx, General, ColMajor> zwork;
116  Matrix<Tcplx, General, ColMajor> workc;
117  Vector<Tcplx> xc(n), yc(n);
118  Vector<T> aq, bq;
119  Vector<Tcplx> zaq, zbq;
120  Vector<T> x(n), y(n);
121 
122  workc.Reallocate(n, m0);
123  if (herm)
124  {
125  work.Reallocate(n, m0);
126  aq.Reallocate(m0*m0);
127  bq.Reallocate(m0*m0);
128  }
129  else if (sym)
130  {
131  zwork.Reallocate(n, m0);
132  zaq.Reallocate(m0*m0);
133  zbq.Reallocate(m0*m0);
134  }
135  else
136  {
137  zwork.Reallocate(n, m0);
138  zaq.Reallocate(m0*m0);
139  zbq.Reallocate(m0*m0);
140  }
141 
142  work.Zero(); zwork.Zero();
143  workc.Zero(); xc.Zero(); yc.Zero();
144  x.Zero(); y.Zero();
145  aq.Zero(); bq.Zero(); zaq.Zero(); zbq.Zero();
146 
147  // output arrays
148  Treal epsout(0); Tcplx ze;
149  Vector<Treal> res(m0);
150 
151  int loop = 0, m = 0;
152  Matrix<Tcplx, General, ColMajor> eigen_vec_cplx;
153  Vector<Treal> lambda; Vector<Tcplx> lambda_cplx;
154 
155  if (herm)
156  {
157  eigen_vectors.Reallocate(n, m0);
158  lambda.Reallocate(m0);
159  }
160  else if (sym)
161  {
162  eigen_vec_cplx.Reallocate(n, m0);
163  lambda_cplx.Reallocate(m0);
164  }
165  else
166  {
167  eigen_vec_cplx.Reallocate(n, 2*m0);
168  lambda_cplx.Reallocate(m0);
169  res.Reallocate(2*m0);
170  }
171 
172  res.Zero();
173  eigen_vectors.Zero(); lambda.Zero();
174  lambda_cplx.Zero(); eigen_vec_cplx.Zero();
175 
176  if (var.DiagonalMass())
177  {
178  // if the mass matrix is diagonal :
179  // the diagonal is computed
180  var.ComputeDiagonalMass();
181 
182  // and M^{-1/2} is computed
183  var.FactorizeDiagonalMass();
184  }
185  else if (var.UseCholeskyFactoForMass())
186  {
187  // if the user wants to use the Cholesky factorization of the mass matrix:
188 
189  // computing mass matrix in a convenient form
190  var.ComputeMassForCholesky();
191 
192  // performing Cholesky factorization
193  var.FactorizeCholeskyMass();
194  }
195  else
196  var.ComputeMassMatrix();
197 
198  // evaluating stiffness matrix K
199  var.ComputeStiffnessMatrix();
200 
201  // main loop (reverse communication interface)
202  int ijob = -1, info = 0;
203  while (ijob != 0)
204  {
205  // feast is called
206  CallFeast(ijob, n, herm, sym, ze, work, workc, aq, bq, zwork, zaq, zbq,
207  fpm, epsout, loop, emin, emax, Emid, r, m0, lambda, eigen_vectors,
208  lambda_cplx, eigen_vec_cplx, m, res, info);
209 
210  if (ijob == 10)
211  {
212  // we have to factorize ze M - K
213  var.ComputeAndFactorizeStiffnessMatrix(ze, -one);
214  }
215  else if (ijob == 11)
216  {
217  // solves (ze K - M) y = workc(n, m0)
218  for (int k = 0; k < fpm(22); k++)
219  {
220  for (int i = 0; i < n; i++)
221  xc(i) = workc(i, k);
222 
223  if (var.DiagonalMass())
224  var.MltSqrtDiagonalMass(xc);
225  else if (var.UseCholeskyFactoForMass())
226  var.MltCholeskyMass(SeldonNoTrans, xc);
227 
228  var.ComputeSolution(xc, yc);
229  var.IncrementLinearSolves();
230 
231  if (var.DiagonalMass())
232  var.MltSqrtDiagonalMass(yc);
233  else if (var.UseCholeskyFactoForMass())
234  {
235  var.MltCholeskyMass(SeldonTrans, yc);
236  var.IncrementLinearSolves();
237  }
238 
239  for (int i = 0; i < n; i++)
240  workc(i, k) = yc(i);
241  }
242  }
243  else if (ijob == 20)
244  {
245  // factorize (ze K - M)^H
246  // already done in ijob = 10
247  }
248  else if (ijob == 21)
249  {
250  // solves (ze K - M)^H y = workc(n, m0)
251  for (int k = 0; k < fpm(22); k++)
252  {
253  for (int i = 0; i < n; i++)
254  xc(i) = workc(i, k);
255 
256  if (var.DiagonalMass())
257  var.MltSqrtDiagonalMass(xc);
258  else if (var.UseCholeskyFactoForMass())
259  var.MltCholeskyMass(SeldonNoTrans, xc);
260 
261  Conjugate(xc);
262  var.ComputeSolution(SeldonTrans, xc, yc);
263  var.IncrementLinearSolves();
264  Conjugate(yc);
265 
266  if (var.DiagonalMass())
267  var.MltSqrtDiagonalMass(yc);
268  else if (var.UseCholeskyFactoForMass())
269  {
270  var.MltCholeskyMass(SeldonTrans, yc);
271  var.IncrementLinearSolves();
272  }
273 
274  for (int i = 0; i < n; i++)
275  workc(i, k) = yc(i);
276  }
277  }
278  else if (ijob == 30)
279  {
280  // multiplication by matrix A
281  int i = fpm(23), j = fpm(23) + fpm(24)-1;
282  if (herm)
283  for (int k = i; k <= j; k++)
284  {
285  GetCol(eigen_vectors, k-1, x);
286 
287  if (var.DiagonalMass())
288  var.MltInvSqrtDiagonalMass(x);
289  else if (var.UseCholeskyFactoForMass())
290  var.SolveCholeskyMass(SeldonTrans, x);
291 
292  var.MltStiffness(x, y);
293  var.IncrementProdMatVect();
294 
295  if (var.DiagonalMass())
296  var.MltInvSqrtDiagonalMass(y);
297  else if (var.UseCholeskyFactoForMass())
298  {
299  var.SolveCholeskyMass(SeldonNoTrans, y);
300  var.IncrementLinearSolves();
301  }
302 
303  for (int p = 0; p < n; p++)
304  work(p, k-1) = y(p);
305  }
306  else
307  for (int k = i; k <= j; k++)
308  {
309  GetCol(eigen_vec_cplx, k-1, xc);
310 
311  if (var.DiagonalMass())
312  var.MltInvSqrtDiagonalMass(xc);
313 
314  var.MltStiffness(xc, yc);
315  var.IncrementProdMatVect();
316 
317  if (var.DiagonalMass())
318  var.MltInvSqrtDiagonalMass(yc);
319 
320  SetCol(yc, k-1, zwork);
321  }
322  }
323  else if (ijob == 31)
324  {
325  // multiplication by matrix A^H
326  int i = fpm(33), j = fpm(33) + fpm(34)-1;
327  if (herm)
328  {
329  cout << "impossible" << endl;
330  abort();
331  }
332  else
333  for (int k = i; k <= j; k++)
334  {
335  GetCol(eigen_vec_cplx, k-1, xc);
336 
337  if (var.DiagonalMass())
338  var.MltInvSqrtDiagonalMass(xc);
339 
340  var.MltStiffness(SeldonConjTrans, xc, yc);
341  var.IncrementProdMatVect();
342 
343  if (var.DiagonalMass())
344  var.MltInvSqrtDiagonalMass(yc);
345 
346  SetCol(yc, k-1, zwork);
347  }
348  }
349  else if (ijob == 40)
350  {
351  // multiplication by matrix B
352  int i = fpm(23), j = fpm(23) + fpm(24)-1;
353  if (herm)
354  for (int k = i; k <= j; k++)
355  {
356  for (int p = 0; p < n; p++)
357  x(p) = eigen_vectors(p, k-1);
358 
359  if (var.DiagonalMass() || var.UseCholeskyFactoForMass())
360  Copy(x, y);
361  else
362  {
363  var.MltMass(x, y);
364  var.IncrementProdMatVect();
365  }
366 
367  for (int p = 0; p < n; p++)
368  work(p, k-1) = y(p);
369  }
370  else
371  for (int k = i; k <= j; k++)
372  {
373  GetCol(eigen_vec_cplx, k-1, xc);
374 
375  if (var.DiagonalMass())
376  Copy(xc, yc);
377  else
378  {
379  var.MltMass(xc, yc);
380  var.IncrementProdMatVect();
381  }
382 
383  SetCol(yc, k-1, zwork);
384  }
385  }
386  else if (ijob == 41)
387  {
388  // multiplication by matrix B^H
389  int i = fpm(33), j = fpm(33) + fpm(34)-1;
390  if (herm)
391  for (int k = i; k <= j; k++)
392  {
393  for (int p = 0; p < n; p++)
394  x(p) = eigen_vectors(p, k-1);
395 
396  if (var.DiagonalMass() || var.UseCholeskyFactoForMass())
397  Copy(x, y);
398  else
399  {
400  var.MltMass(SeldonConjTrans, x, y);
401  var.IncrementProdMatVect();
402  }
403 
404  for (int p = 0; p < n; p++)
405  work(p, k-1) = y(p);
406  }
407  else
408  for (int k = i; k <= j; k++)
409  {
410  GetCol(eigen_vec_cplx, k-1, xc);
411 
412  if (var.DiagonalMass())
413  Copy(xc, yc);
414  else
415  {
416  var.MltMass(SeldonConjTrans, xc, yc);
417  var.IncrementProdMatVect();
418  }
419 
420  SetCol(yc, k-1, zwork);
421  }
422  }
423  }
424 
425  if (info != 0)
426  {
427  if (info == 3)
428  {
429  cout << "Feast returns error code 3 " << endl;
430  cout << "It usually means that there are more than "
431  << m0-1 << " eigenvalues in the ";
432 
433  if (sym)
434  cout << "interval [" << emin << ", " << emax << "] " << endl;
435  else
436  cout << "disk of center " << Emid << " and radius" << r << endl;
437 
438  cout << "You can try to change the ";
439  if (sym)
440  cout << "interval ";
441  else
442  cout << "disk ";
443 
444  cout << "or increase the number of asked eigenvalues" << endl;
445 
446  abort();
447  }
448  else if (info == 1)
449  {
450  cout << "No eigenvalues have been found in the ";
451  if (sym)
452  cout << "interval [" << emin << ", " << emax << "] " << endl;
453  else
454  cout << "disk of center " << Emid << " and radius" << r << endl;
455 
456  cout << "You can try to change the ";
457  if (sym)
458  cout << "interval. " << endl;
459  else
460  cout << "disk. " << endl;
461 
462  abort();
463  }
464  else if (info == 2)
465  {
466  cout << "Feast did not converge, loop = " << loop << endl;
467  cout << "Maximum number of iterations = " << var.GetNbMaximumIterations() << endl;
468  cout << "Residual = " << res << endl;
469  }
470  else if (info == 5)
471  {
472  if (var.GetGlobalRankCommunicator() == 0)
473  cout << "Estimated number of eigenvalues = " << m << endl;
474 
475  return;
476  }
477  else if (info == 6)
478  {
479  if (var.GetGlobalRankCommunicator() == 0)
480  cout << "Warning : subspace is not biorthogonal" << endl;
481 
482  /*if (sym)
483  for (int i = 0; i < m; i++)
484  cout << "Errors on lambda " << i << " = " << res(i) << endl;
485  else
486  for (int i = 0; i < m; i++)
487  cout << "Errors on lambda " << i << " = " << res(i) << " " << res(i+m0) << endl;
488  */
489  }
490  else if (info > 6)
491  {
492  cout << "An error occurred during feast eigenvalue solving " << endl;
493  cout << "info = " << info << endl;
494  abort();
495  }
496  }
497 
498  if ((var.GetPrintLevel() > 0) && (var.GetGlobalRankCommunicator() == 0))
499  cout << "Feast converged in " << loop << " cycles" << endl;
500 
501  work.Clear(); workc.Clear(); zwork.Clear();
502  aq.Clear(); bq.Clear(); zaq.Clear(); zbq.Clear();
503 
504  eigen_values.Reallocate(m);
505  lambda_imag.Reallocate(m);
506  lambda_imag.Zero();
507  if (herm)
508  {
509  for (int i = 0; i < m; i++)
510  eigen_values(i) = lambda(i);
511 
512  eigen_vectors.Resize(n, m);
513  }
514  else
515  var.FillComplexEigenvectors(m, Emid, tol, lambda_cplx, eigen_vec_cplx,
516  eigen_values, lambda_imag, eigen_vectors);
517 
518  lambda.Clear();
519  lambda_cplx.Clear(); eigen_vec_cplx.Clear();
520 
521  T shiftr = var.GetShiftValue(), shifti = var.GetImagShiftValue();
522  ApplyScalingEigenvec(var, eigen_values, lambda_imag, eigen_vectors,
523  shiftr, shifti);
524 
525  }
526 
527 
528  void CallFeast(int& ijob, int n, bool herm, bool sym, complex<double>& ze,
529  Matrix<double, General, ColMajor>& work,
530  Matrix<complex<double>, General, ColMajor>& workc,
531  Vector<double>& aq, Vector<double>& bq,
532  Matrix<complex<double>, General, ColMajor>& zwork,
533  Vector<complex<double> >& zaq, Vector<complex<double> >& zbq,
534  IVect& fpm, double& epsout, int& loop,
535  double emin, double emax, complex<double> Emid, double r,
536  int m0, Vector<double>& lambda,
537  Matrix<double, General, ColMajor>& eigen_vectors,
538  Vector<complex<double> >& lambda_cplx,
539  Matrix<complex<double>, General, ColMajor>& eigen_vec_cplx,
540  int& m, Vector<double>& res, int& info)
541  {
542  if (sym)
543  dfeast_srci_(&ijob, &n, reinterpret_cast<void*>(&ze), work.GetData(),
544  workc.GetDataVoid(), aq.GetData(), bq.GetData(),
545  fpm.GetData(), &epsout, &loop, &emin, &emax,
546  &m0, lambda.GetData(), eigen_vectors.GetData(),
547  &m, res.GetData(), &info);
548  else
549  dfeast_grci_(&ijob, &n, reinterpret_cast<void*>(&ze), zwork.GetDataVoid(),
550  workc.GetDataVoid(), zaq.GetDataVoid(), zbq.GetDataVoid(),
551  fpm.GetData(), &epsout, &loop, &Emid, &r,
552  &m0, lambda_cplx.GetDataVoid(), eigen_vec_cplx.GetDataVoid(),
553  &m, res.GetData(), &info);
554  }
555 
556  void CallFeast(int& ijob, int n, bool herm, bool sym, complex<double>& ze,
557  Matrix<complex<double>, General, ColMajor>& work,
558  Matrix<complex<double>, General, ColMajor>& workc,
559  Vector<complex<double> >& aq, Vector<complex<double> >& bq,
560  Matrix<complex<double>, General, ColMajor>& zwork,
561  Vector<complex<double> >& zaq, Vector<complex<double> >& zbq,
562  IVect& fpm, double& epsout, int& loop,
563  double emin, double emax, complex<double> Emid, double r,
564  int m0, Vector<double>& lambda,
565  Matrix<complex<double>, General, ColMajor>& eigen_vectors,
566  Vector<complex<double> >& lambda_cplx,
567  Matrix<complex<double>, General, ColMajor>& eigen_vec_cplx,
568  int& m, Vector<double>& res, int& info)
569  {
570  if (herm)
571  zfeast_hrci_(&ijob, &n, reinterpret_cast<void*>(&ze), work.GetDataVoid(),
572  workc.GetDataVoid(), aq.GetDataVoid(), bq.GetDataVoid(),
573  fpm.GetData(), &epsout, &loop, &emin, &emax, &m0,
574  lambda.GetData(), eigen_vectors.GetDataVoid(),
575  &m, res.GetData(), &info);
576  else if (sym)
577  zfeast_srci_(&ijob, &n, reinterpret_cast<void*>(&ze), zwork.GetDataVoid(),
578  workc.GetDataVoid(), zaq.GetDataVoid(), zbq.GetDataVoid(),
579  fpm.GetData(), &epsout, &loop, &Emid, &r, &m0,
580  lambda_cplx.GetData(), eigen_vec_cplx.GetDataVoid(),
581  &m, res.GetData(), &info);
582  else
583  zfeast_grci_(&ijob, &n, reinterpret_cast<void*>(&ze), zwork.GetDataVoid(),
584  workc.GetDataVoid(), zaq.GetDataVoid(), zbq.GetDataVoid(),
585  fpm.GetData(), &epsout, &loop, &Emid, &r, &m0,
586  lambda_cplx.GetData(), eigen_vec_cplx.GetDataVoid(),
587  &m, res.GetData(), &info);
588  }
589 
590  void CallFeastP(int& ijob, int p, int n, bool sym, complex<double>& ze,
591  Matrix<complex<double>, General, ColMajor>& zwork,
592  Matrix<complex<double>, General, ColMajor>& workc,
593  Vector<complex<double> >& zaq, Vector<complex<double> >& zbq,
594  IVect& fpm, double& epsout, int& loop, complex<double> Emid, double r,
595  int m0, Vector<complex<double> >& lambda_cplx,
596  Matrix<complex<double>, General, ColMajor>& eigen_vec_cplx,
597  int& m, Vector<double>& res, int& info)
598  {
599  if (sym)
600  zfeast_srcipev_(&ijob, &p, &n, reinterpret_cast<void*>(&ze),
601  zwork.GetDataVoid(), workc.GetDataVoid(),
602  zaq.GetDataVoid(), zbq.GetDataVoid(),
603  fpm.GetData(), &epsout, &loop, &Emid, &r, &m0,
604  lambda_cplx.GetData(), eigen_vec_cplx.GetDataVoid(),
605  &m, res.GetData(), &info);
606  else
607  zfeast_grcipev_(&ijob, &p, &n, reinterpret_cast<void*>(&ze),
608  zwork.GetDataVoid(), workc.GetDataVoid(),
609  zaq.GetDataVoid(), zbq.GetDataVoid(),
610  fpm.GetData(), &epsout, &loop, &Emid, &r, &m0,
611  lambda_cplx.GetData(), eigen_vec_cplx.GetDataVoid(),
612  &m, res.GetData(), &info);
613  }
614 
615 #ifdef SELDON_WITH_VIRTUAL
616  template<class T>
617  void FindEigenvaluesFeast(PolynomialEigenProblem_Base<T>& var,
618  Vector<T>& eigen_values,
619  Vector<T>& lambda_imag,
620  Matrix<T, General, ColMajor>& eigen_vectors)
621  {
622  // initialization of feast parameters
623  IVect fpm(128);
624  fpm.Fill(0);
625  feastinit_(fpm.GetData());
626 
627  fpm(0) = 1;
628  if (var.GetPrintLevel() > 0)
629  fpm(0) = 1;
630 
631  double tol = var.GetStoppingCriterion();
632  fpm(2) = -log10(tol);
633 
634  for (int i = 64; i < fpm.GetM(); i++)
635  fpm(i) = 0;
636 
637  fpm(3) = var.GetNbMaximumIterations();
638 
639  FeastParam& param = var.GetFeastParameters();
640  if (param.EstimateNumberEigenval())
641  fpm(13) = 2;
642 
643 #ifdef SELDON_WITH_MPI
644  MPI_Comm solver_comm; // duplication needed for feast
645  MPI_Comm_dup(var.GetCommunicator(), &solver_comm);
646 
647  fpm(8) = MPI_Comm_c2f(var.GetRciCommunicator());
648  fpm(48) = MPI_Comm_c2f(solver_comm);
649  fpm(58) = MPI_Comm_c2f(var.GetGlobalCommunicator());
650 #endif
651 
652  // initialization of computation
653  int n = var.GetM();
654  int p = var.GetPolynomialDegree();
655 
656  typedef typename ClassComplexType<T>::Tcplx Tcplx;
657  typedef typename ClassComplexType<T>::Treal Treal;
658 
659  int m0 = var.GetNbAskedEigenvalues()+1;
660 
661  Tcplx Emid = param.GetCircleCenterSpectrum();
662  Treal r = param.GetCircleRadiusSpectrum();
663  fpm(17) = param.GetRatioEllipseSpectrum();
664  fpm(18) = param.GetAngleEllipseSpectrum();
665 
666  bool sym = var.IsSymmetricProblem();
667  if (sym)
668  fpm(14) = 2;
669  else
670  fpm(14) = 1;
671 
672  if (param.GetTypeIntegration() >= 0)
673  fpm(15) = param.GetTypeIntegration();
674 
675  if (param.GetNumOfQuadraturePoints() > 0)
676  fpm(7) = param.GetNumOfQuadraturePoints();
677  else
678  fpm(7) = 16;
679 
680  // work arrays
681  Matrix<Tcplx, General, ColMajor> zwork;
682  Matrix<Tcplx, General, ColMajor> workc;
683  Vector<Tcplx> xc(n), yc(n);
684  Vector<Tcplx> zaq, zbq;
685 
686  workc.Reallocate(n, m0);
687  zwork.Reallocate(n, m0);
688  zaq.Reallocate(p*p*m0*m0);
689  zbq.Reallocate(p*p*m0*m0);
690 
691  zwork.Zero();
692  workc.Zero(); xc.Zero(); yc.Zero();
693  zaq.Zero(); zbq.Zero();
694 
695  // output arrays
696  Treal epsout(0); Tcplx ze;
697  Vector<Treal> res(m0);
698 
699  int loop = 0, m = 0;
700  Matrix<Tcplx, General, ColMajor> eigen_vec_cplx;
701  Vector<Treal> lambda; Vector<Tcplx> lambda_cplx;
702 
703  eigen_vec_cplx.Reallocate(n, m0);
704  lambda_cplx.Reallocate(m0);
705 
706  res.Zero();
707  eigen_vectors.Zero(); lambda.Zero();
708  lambda_cplx.Zero(); eigen_vec_cplx.Zero();
709 
710  Vector<Tcplx> coef(p+1); coef.Zero();
711 
712  // main loop (reverse communication interface)
713  int ijob = -1, info = 0;
714  while (ijob != 0)
715  {
716  // feast is called
717  CallFeastP(ijob, p, n, sym, ze, zwork, workc, zaq, zbq,
718  fpm, epsout, loop, Emid, r, m0, lambda_cplx, eigen_vec_cplx,
719  m, res, info);
720 
721  if (ijob == 10)
722  {
723  // we have to factorize P(Ze)
724  SetComplexOne(coef(0));
725  for (int i = 0; i < p; i++)
726  coef(i+1) = coef(i)*ze;
727 
728  var.FactorizeOperator(coef);
729  }
730  else if (ijob == 11)
731  {
732  // solves P(ze) y = workc(n, m0)
733  for (int k = 0; k < fpm(22); k++)
734  {
735  for (int i = 0; i < n; i++)
736  xc(i) = workc(i, k);
737 
738  var.SolveOperator(SeldonNoTrans, xc, yc);
739  var.IncrementLinearSolves();
740 
741  for (int i = 0; i < n; i++)
742  workc(i, k) = yc(i);
743  }
744  }
745  else if (ijob == 20)
746  {
747  // factorize P(Ze)^H
748  // already done in ijob = 10
749  }
750  else if (ijob == 21)
751  {
752  // solves P(ze)^H y = workc(n, m0)
753  for (int k = 0; k < fpm(22); k++)
754  {
755  for (int i = 0; i < n; i++)
756  xc(i) = workc(i, k);
757 
758  Conjugate(xc);
759  var.SolveOperator(SeldonTrans, xc, yc);
760  var.IncrementLinearSolves();
761  Conjugate(yc);
762 
763  for (int i = 0; i < n; i++)
764  workc(i, k) = yc(i);
765  }
766  }
767  else if (ijob == 30)
768  {
769  // multiplication by matrix A[fpm(56)]
770  int i = fpm(23), j = fpm(23) + fpm(24)-1;
771  for (int k = i; k <= j; k++)
772  {
773  GetCol(eigen_vec_cplx, k-1, xc);
774 
775  var.MltOperator(fpm(56)-1, SeldonNoTrans, xc, yc);
776  var.IncrementProdMatVect();
777 
778  SetCol(yc, k-1, zwork);
779  }
780  }
781  else if (ijob == 31)
782  {
783  // multiplication by matrix A^H
784  int i = fpm(33), j = fpm(33) + fpm(34)-1;
785  for (int k = i; k <= j; k++)
786  {
787  GetCol(eigen_vec_cplx, k-1, xc);
788 
789  var.MltOperator(fpm(56)-1, SeldonConjTrans, xc, yc);
790  var.IncrementProdMatVect();
791 
792  SetCol(yc, k-1, zwork);
793  }
794  }
795  }
796  zwork.Clear(); workc.Clear(); zaq.Clear(); zbq.Clear();
797 
798  var.FillComplexEigenvectors(m, Emid, tol, lambda_cplx, eigen_vec_cplx,
799  eigen_values, lambda_imag, eigen_vectors);
800 
801  eigen_vec_cplx.Clear(); lambda_cplx.Clear();
802 
803  var.DistributeEigenvectors(eigen_vectors);
804  }
805 #endif
806 
807 }
808 
809 #endif
Seldon::Conjugate
void Conjugate(Matrix< T, Prop, Storage, Allocator > &A)
A is replaced by its conjugate.
Definition: Functions_Matrix.cxx:2915
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::ApplyScalingEigenvec
void ApplyScalingEigenvec(EigenPb &var, Vector1 &eigen_values, Vector1 &lambda_imag, Matrix1 &eigen_vectors, const T0 &shiftr, const T0 &shifti)
modification of eigenvectors to take into account scaling by mass matrix
Definition: EigenvalueSolver.cxx:846