Slepc.cxx
1 #ifndef SELDON_FILE_SLEPC_CXX
2 
3 #include <slepceps.h>
4 #include <slepcpep.h>
5 #include <slepcnep.h>
6 #include <slepc/private/stimpl.h>
7 #include "petsc/src/mat/impls/aij/seq/aij.h"
8 
9 typedef struct {
10  PetscInt nmat,maxnmat;
11  PetscScalar *coeff;
12  Mat *A;
13  Vec t;
14 } NEP_NLEIGS_MATSHELL;
15 
16 
17 namespace Seldon
18 {
19 
20  const char* SlepcParam::GetEigensolverChar() const
21  {
22  switch (type_solver)
23  {
24  case POWER : return EPSPOWER;
25  case SUBSPACE : return EPSSUBSPACE;
26  case ARNOLDI : return EPSARNOLDI;
27  case LANCZOS : return EPSLANCZOS;
28  case KRYLOVSCHUR : return EPSKRYLOVSCHUR;
29  case GD : return EPSGD;
30  case JD : return EPSJD;
31  case RQCG : return EPSRQCG;
32  case LOBPCG : return EPSLOBPCG;
33  case CISS : return EPSCISS;
34  case LAPACK : return EPSLAPACK;
35  case ARPACK : return EPSARPACK;
36  case BLZPACK : return EPSBLZPACK;
37  case TRLAN : return EPSTRLAN;
38  case BLOPEX : return EPSBLOPEX;
39  case PRIMME : return EPSPRIMME;
40  //case FEAST : return EPSFEAST;
41  }
42 
43  return "";
44  }
45 
46 
47  void SetParametersSlepc(const SlepcParam& param, EPS& eps)
48  {
49  int ierr = EPSSetType(eps, param.GetEigensolverChar());
50  if (ierr != 0)
51  {
52  cout << "Chosen type = " << param.GetEigensolverChar() << " Not enabled in Slepc ?" << endl;
53  abort();
54  }
55 
56  int solver = param.GetEigensolverType();
57  if (solver == param.BLOPEX)
58  {
59 #ifdef SLEPC_HAVE_BLOPEX
60  PetscInt bs = param.GetBlockSize();
61  if (bs > 0)
62  EPSBLOPEXSetBlockSize(eps, bs);
63 #else
64  cout << "Slepc not compiled with BLOPEX" << endl;
65  abort();
66 #endif
67  }
68  else if (solver == param.BLZPACK)
69  {
70 #ifdef SLEPC_HAVE_BLZPACK
71  PetscInt bs = param.GetBlockSize();
72  if (bs > 0)
73  EPSBlzpackSetBlockSize(eps, bs);
74 
75  bs = param.GetNumberOfSteps();
76  if (bs > 0)
77  EPSBlzpackSetNSteps(eps, bs);
78 #else
79  cout << "Slepc not compiled with BLZPACK" << endl;
80  abort();
81 #endif
82  }
83  else if (solver == param.CISS)
84  {
85  int type_extraction = param.GetExtractionType();
86  if (type_extraction >= 0)
87  {
88  if (type_extraction == param.EXTRACT_RITZ)
89  EPSCISSSetExtraction(eps, EPS_CISS_EXTRACTION_RITZ);
90  else
91  EPSCISSSetExtraction(eps, EPS_CISS_EXTRACTION_HANKEL);
92  }
93 
94  int type_quad = param.GetQuadratureRuleType();
95  if (type_quad >= 0)
96  {
97  if (type_quad == param.QUADRULE_TRAPEZE)
98  EPSCISSSetQuadRule(eps, EPS_CISS_QUADRULE_TRAPEZOIDAL);
99  else
100  EPSCISSSetQuadRule(eps, EPS_CISS_QUADRULE_CHEBYSHEV);
101  }
102 
103  PetscScalar a; PetscBool complex_number = PETSC_FALSE;
104  if (IsComplexNumber(a))
105  complex_number = PETSC_TRUE;
106 
107  if (param.GetInnerSteps() > 0)
108  EPSCISSSetRefinement(eps, param.GetInnerSteps(), param.GetOuterSteps());
109 
110  if (param.GetNumberIntegrationPoints() > 0)
111  EPSCISSSetSizes(eps, param.GetNumberIntegrationPoints(), param.GetBlockSize(),
112  param.GetMomentSize(), param.GetNumberPartitions(),
113  param.GetMaximumBlockSize(), complex_number);
114 
115  if (param.GetThresholdRank() > 0)
116  EPSCISSSetThreshold(eps, param.GetThresholdRank(), param.GetThresholdSpurious());
117  }
118  else if (solver == param.FEAST)
119  {
120 #ifdef SLEPC_HAVE_FEAST
121  if (param.GetNumberIntegrationPoints() > 0)
122  EPSFEASTSetNumPoints(eps, param.GetNumberIntegrationPoints());
123 #else
124  cout << "Slepc not compiled with FEAST" << endl;
125  abort();
126 #endif
127  }
128  else if (solver == param.GD)
129  {
130  if (param.GetBorthogonalization() >= 0)
131  {
132  PetscBool borth = PETSC_FALSE;
133  if (param.GetBorthogonalization() >= 1)
134  borth = PETSC_TRUE;
135 
136  EPSGDSetBOrth(eps, borth);
137  }
138 
139  PetscInt bs = param.GetBlockSize();
140  if (bs > 0)
141  EPSGDSetBlockSize(eps, bs);
142 
143  if (param.GetDoubleExpansion() >= 0)
144  {
145  PetscBool exp = PETSC_FALSE;
146  if (param.GetDoubleExpansion() >= 1)
147  exp = PETSC_TRUE;
148 
149  EPSGDSetDoubleExpansion(eps, exp);
150  }
151 
152  if (param.GetInitialSize() > 0)
153  EPSGDSetInitialSize(eps, param.GetInitialSize());
154 
155  if (param.GetKrylovRestart() >= 0)
156  {
157  PetscBool restart = PETSC_FALSE;
158  if (param.GetKrylovRestart() >= 1)
159  restart = PETSC_TRUE;
160 
161  EPSGDSetKrylovStart(eps, restart);
162  }
163 
164  if (param.GetRestartNumber() > 0)
165  EPSGDSetRestart(eps, param.GetRestartNumber(), param.GetRestartNumberAdd());
166 
167  // deprecated function : SetWindowSizes
168  //if (param.GetNumberConvergedVectors() > 0)
169  // EPSGDSetWindowSizes(eps, param.GetNumberConvergedVectors(),
170  // param.GetNumberConvergedVectorsProjected());
171  }
172  else if (solver == param.JD)
173  {
174  if (param.GetBorthogonalization() >= 0)
175  {
176  PetscBool borth = PETSC_FALSE;
177  if (param.GetBorthogonalization() >= 1)
178  borth = PETSC_TRUE;
179 
180  EPSJDSetBOrth(eps, borth);
181  }
182 
183  PetscInt bs = param.GetBlockSize();
184  if (bs > 0)
185  EPSJDSetBlockSize(eps, bs);
186 
187  if (param.GetInitialSize() > 0)
188  EPSJDSetInitialSize(eps, param.GetInitialSize());
189 
190  if (param.GetKrylovRestart() >= 0)
191  {
192  PetscBool restart = PETSC_FALSE;
193  if (param.GetKrylovRestart() >= 1)
194  restart = PETSC_TRUE;
195 
196  EPSJDSetKrylovStart(eps, restart);
197  }
198 
199  if (param.GetRestartNumber() > 0)
200  EPSJDSetRestart(eps, param.GetRestartNumber(), param.GetRestartNumberAdd());
201 
202  // deprecated function : SetWindowSizes
203  //if (param.GetNumberConvergedVectors() > 0)
204  // EPSJDSetWindowSizes(eps, param.GetNumberConvergedVectors(),
205  // param.GetNumberConvergedVectorsProjected());
206  }
207  else if (solver == param.KRYLOVSCHUR)
208  {
209  if (param.UseNonLockingVariant())
210  EPSKrylovSchurSetLocking(eps, PETSC_FALSE);
211  else
212  EPSKrylovSchurSetLocking(eps, PETSC_TRUE);
213 
214  if (param.GetRestartRatio() > 0)
215  EPSKrylovSchurSetRestart(eps, param.GetRestartRatio());
216  }
217  else if (solver == param.LOBPCG)
218  {
219  if (param.GetBlockSize() > 0)
220  EPSLOBPCGSetBlockSize(eps, param.GetBlockSize());
221 
222  if (param.UseNonLockingVariant())
223  EPSLOBPCGSetLocking(eps, PETSC_FALSE);
224  else
225  EPSLOBPCGSetLocking(eps, PETSC_TRUE);
226 
227  if (param.GetRestartRatio() > 0)
228  EPSLOBPCGSetRestart(eps, param.GetRestartRatio());
229  }
230  else if (solver == param.PRIMME)
231  {
232 #ifdef SLEPC_HAVE_PRIMME
233  if (param.GetBlockSize() > 0)
234  EPSPRIMMESetBlockSize(eps, param.GetBlockSize());
235 
236  if (param.GetMethod().size() > 1)
237  {
238  if (param.GetMethod() == "DYNAMIC")
239  EPSPRIMMESetMethod(eps, EPS_PRIMME_DYNAMIC);
240  else if (param.GetMethod() == "DEFAULT_MIN_TIME")
241  EPSPRIMMESetMethod(eps, EPS_PRIMME_DEFAULT_MIN_TIME);
242  else if (param.GetMethod() == "DEFAULT_MIN_MATVECS")
243  EPSPRIMMESetMethod(eps, EPS_PRIMME_DEFAULT_MIN_MATVECS);
244  else if (param.GetMethod() == "ARNOLDI")
245  EPSPRIMMESetMethod(eps, EPS_PRIMME_ARNOLDI);
246  else if (param.GetMethod() == "GD")
247  EPSPRIMMESetMethod(eps, EPS_PRIMME_GD);
248  else if (param.GetMethod() == "GD_PLUSK")
249  EPSPRIMMESetMethod(eps, EPS_PRIMME_GD_PLUSK);
250  else if (param.GetMethod() == "GD_OLSEN_PLUSK")
251  EPSPRIMMESetMethod(eps, EPS_PRIMME_GD_OLSEN_PLUSK);
252  else if (param.GetMethod() == "JD_OLSEN_PLUSK")
253  EPSPRIMMESetMethod(eps, EPS_PRIMME_JD_OLSEN_PLUSK);
254  else if (param.GetMethod() == "RQI")
255  EPSPRIMMESetMethod(eps, EPS_PRIMME_RQI);
256  else if (param.GetMethod() == "JDQR")
257  EPSPRIMMESetMethod(eps, EPS_PRIMME_JDQR);
258  else if (param.GetMethod() == "JDQMR")
259  EPSPRIMMESetMethod(eps, EPS_PRIMME_JDQMR);
260  else if (param.GetMethod() == "JDQMR_ETOL")
261  EPSPRIMMESetMethod(eps, EPS_PRIMME_JDQMR_ETOL);
262  else if (param.GetMethod() == "SUBSPACE_ITERATION")
263  EPSPRIMMESetMethod(eps, EPS_PRIMME_SUBSPACE_ITERATION);
264  else if (param.GetMethod() == "LOBPCG_ORTHOBASIS")
265  EPSPRIMMESetMethod(eps, EPS_PRIMME_LOBPCG_ORTHOBASIS);
266  else if (param.GetMethod() == "LOBPCG_ORTHOBASISW")
267  EPSPRIMMESetMethod(eps, EPS_PRIMME_LOBPCG_ORTHOBASISW);
268  }
269 #else
270  cout << "Slepc not compiled with PRIMME" << endl;
271  abort();
272 #endif
273  }
274  else if (solver == param.POWER)
275  {
276  if (param.GetShiftType() >= 0)
277  {
278  if (param.GetShiftType() == param.SHIFT_CONSTANT)
279  EPSPowerSetShiftType(eps, EPS_POWER_SHIFT_CONSTANT);
280  else if (param.GetShiftType() == param.SHIFT_RAYLEIGH)
281  EPSPowerSetShiftType(eps, EPS_POWER_SHIFT_RAYLEIGH);
282  else if (param.GetShiftType() == param.SHIFT_WILKINSON)
283  EPSPowerSetShiftType(eps, EPS_POWER_SHIFT_WILKINSON);
284  }
285  }
286  else if (solver == param.RQCG)
287  {
288  if (param.GetNumberOfSteps() > 0)
289  EPSRQCGSetReset(eps, param.GetNumberOfSteps());
290  }
291  }
292 
293 
295  void CopyPointerPetsc(const Vec& x, Vector<PetscScalar>& y)
296  {
297  // it is assumed that the vector x is stored in a contiguous array
298  PetscInt n;
299  VecGetLocalSize(x, &n);
300 
301  PetscScalar* x_array;
302  VecGetArrayRead(x, (const PetscScalar**)&x_array);
303 
304  y.SetData(n, x_array);
305  }
306 
307 
308  void AllocatePetscVector(const MPI_Comm& comm, Vec& Vr, int n, int nglob,
309  Vector<PetscScalar>& Vr_vec)
310  {
311  VecCreate(comm, &Vr);
312  VecSetSizes(Vr, n, nglob);
313  VecSetFromOptions(Vr);
314 
315  CopyPointerPetsc(Vr, Vr_vec);
316  }
317 
318 
320  PetscErrorCode MatMult_Matrix(Mat mat, Vec x, Vec y)
321  {
322  Vector<PetscScalar> xvec0, yvec;
323  CopyPointerPetsc(x, xvec0);
324  CopyPointerPetsc(y, yvec);
325 
326  void* ctx;
327  MatShellGetContext(mat, &ctx);
328 
329  EigenProblem_Base<PetscScalar>& var_eig
330  = *reinterpret_cast<EigenProblem_Base<PetscScalar>* >(ctx);
331 
332  var_eig.IncrementProdMatVect();
333  Vector<PetscScalar> xvec(xvec0);
334  if (var_eig.DiagonalMass() || (var_eig.UseCholeskyFactoForMass()))
335  {
336  // standard eigenvalue problem
337  if (var_eig.GetComputationalMode() == var_eig.REGULAR_MODE)
338  {
339  if (var_eig.DiagonalMass())
340  var_eig.MltInvSqrtDiagonalMass(xvec);
341  else
342  var_eig.SolveCholeskyMass(Seldon::SeldonTrans, xvec);
343 
344  var_eig.MltStiffness(xvec, yvec);
345  var_eig.IncrementProdMatVect();
346 
347  if (var_eig.DiagonalMass())
348  var_eig.MltInvSqrtDiagonalMass(yvec);
349  else
350  {
351  var_eig.SolveCholeskyMass(Seldon::SeldonNoTrans, yvec);
352  var_eig.IncrementLinearSolves();
353  }
354  }
355  else
356  {
357  if (var_eig.DiagonalMass())
358  var_eig.MltSqrtDiagonalMass(xvec);
359  else
360  var_eig.MltCholeskyMass(Seldon::SeldonNoTrans, xvec);
361 
362  var_eig.ComputeSolution(xvec, yvec);
363  var_eig.IncrementLinearSolves();
364 
365  if (var_eig.DiagonalMass())
366  var_eig.MltSqrtDiagonalMass(yvec);
367  else
368  {
369  var_eig.MltCholeskyMass(Seldon::SeldonTrans, yvec);
370  var_eig.IncrementLinearSolves();
371  }
372  }
373  }
374  else
375  {
376  if (var_eig.GetComputationalMode() == var_eig.INVERT_MODE)
377  {
378  if (var_eig.GetTypeSpectrum() != var_eig.CENTERED_EIGENVALUES)
379  {
380  var_eig.MltStiffness(xvec0, xvec);
381  var_eig.ComputeSolution(xvec, yvec);
382  }
383  else
384  {
385  var_eig.MltMass(xvec0, xvec);
386  var_eig.ComputeSolution(xvec, yvec);
387  }
388 
389  var_eig.IncrementLinearSolves();
390  var_eig.IncrementProdMatVect();
391  }
392  else if (var_eig.GetComputationalMode() == var_eig.REGULAR_MODE)
393  {
394  var_eig.MltStiffness(xvec, yvec);
395  var_eig.IncrementProdMatVect();
396  }
397  else
398  {
399  cout << "Not implemented" << endl;
400  abort();
401  }
402  }
403 
404  xvec0.Nullify(); yvec.Nullify();
405  return 0;
406  }
407 
408 
410  PetscErrorCode MatMult_MassMatrix(Mat mat, Vec x, Vec y)
411  {
412  Vector<PetscScalar> xvec, yvec;
413  CopyPointerPetsc(x, xvec);
414  CopyPointerPetsc(y, yvec);
415 
416  void* ctx;
417  MatShellGetContext(mat, &ctx);
418 
419  EigenProblem_Base<PetscScalar>& var_eig
420  = *reinterpret_cast<EigenProblem_Base<PetscScalar>* >(ctx);
421 
422  var_eig.MltMass(xvec, yvec);
423  var_eig.IncrementProdMatVect();
424 
425  xvec.Nullify(); yvec.Nullify();
426  return 0;
427  }
428 
429 
431  PetscErrorCode MatMultTranspose_Matrix(Mat A, Vec x, Vec y)
432  {
433  throw Undefined("Function MatMultTranspose_Matrix not implemented");
434  return 0;
435  }
436 
437 
439  PetscErrorCode MatGetDiagonal_Matrix(Mat A, Vec d)
440  {
441  throw Undefined("Function MatGetDiagonal_Matrix not implemented");
442  return 0;
443  }
444 
445 
446  template<class T>
447  bool PutEigenpairLapackForm(int num, int nev, T& Lr, T& Li, Vector<T>& Vr, Vector<T>& Vi,
448  T& Lr_next, T& Li_next, Vector<T>& Vr_next, Vector<T>& Vi_next,
449  Vector<T>& eigen_values, Vector<T>& lambda_imag,
450  Matrix<T, General, ColMajor>& eigen_vectors)
451  {
452  bool eigen_pair = false;
453  if ((Li != T(0)) && (num < nev-1))
454  eigen_pair = true;
455 
456  int n = Vr.GetM();
457  if (eigen_pair)
458  {
459  eigen_values(num) = Lr;
460  lambda_imag(num) = Li;
461  eigen_values(num+1) = Lr_next;
462  lambda_imag(num+1) = Li_next;
463  for (int j = 0; j < n; j++)
464  {
465  eigen_vectors(j, num) = Vr(j);
466  eigen_vectors(j, num+1) = Vi(j);
467  }
468  }
469  else
470  {
471  eigen_values(num) = Lr;
472  lambda_imag(num) = Li;
473  for (int j = 0; j < n; j++)
474  eigen_vectors(j, num) = Vr(j);
475  }
476 
477  return eigen_pair;
478  }
479 
480 
481  template<class T>
482  bool PutEigenpairLapackForm(int num, int nev, complex<T>& Lr, complex<T>& Li,
483  Vector<complex<T> >& Vr, Vector<complex<T> >& Vi,
484  complex<T>& Lr_next, complex<T>& Li_next,
485  Vector<complex<T> >& Vr_next, Vector<complex<T> >& Vi_next,
486  Vector<complex<T> >& eigen_values, Vector<complex<T> >& lambda_imag,
487  Matrix<complex<T>, General, ColMajor>& eigen_vectors)
488  {
489  int n = Vr.GetM();
490  eigen_values(num) = Lr;
491  lambda_imag(num) = Li;
492  for (int j = 0; j < n; j++)
493  eigen_vectors(j, num) = Vr(j);
494 
495  return false;
496  }
497 
498 
499  void FindEigenvaluesSlepc_(EigenProblem_Base<PetscScalar>& var,
500  Vector<PetscScalar>& eigen_values,
501  Vector<PetscScalar>& lambda_imag,
502  Matrix<PetscScalar, General, ColMajor>& eigen_vectors)
503  {
504  // initializing of computation
505  PetscScalar shiftr = var.GetShiftValue(), shifti = var.GetImagShiftValue();
506  PetscScalar zero, one;
507  SetComplexZero(zero);
508  SetComplexOne(one);
509 
510  int print_level = var.GetPrintLevel();
511  SlepcParam& param = var.GetSlepcParameters();
512  int rank_proc; MPI_Comm_rank(var.GetCommunicator(), &rank_proc);
513 
514  Mat stiff, mass;
515  int nev = var.GetNbAskedEigenvalues();
516  int ncv = var.GetNbArnoldiVectors();
517  int n = var.GetM();
518  int nglob = var.GetGlobalM();
519 
520  // creation of a matrix-free Petsc structure
521  MatCreateShell(var.GetCommunicator(), n, n, nglob, nglob,
522  reinterpret_cast<void*>(&var), &stiff);
523 
524  MatCreateShell(var.GetCommunicator(), n, n, nglob, nglob,
525  reinterpret_cast<void*>(&var), &mass);
526 
527  //MatSetFromOptions(stiff); MatSetFromOptions(mass);
528  MatShellSetOperation(stiff, MATOP_MULT, (void(*)())MatMult_Matrix);
529  MatShellSetOperation(stiff, MATOP_MULT_TRANSPOSE, (void(*)())MatMultTranspose_Matrix);
530  //MatShellSetOperation(stiff, MATOP_GET_DIAGONAL, (void(*)())MatGetDiagonal_Matrix);
531 
532  MatShellSetOperation(mass, MATOP_MULT, (void(*)())MatMult_MassMatrix);
533 
534  // creation of the eigensolver
535  EPS eps;
536  EPSCreate(var.GetCommunicator(), &eps);
537  SetParametersSlepc(param, eps);
538 
539  // type of eigenproblem (hermitian/generalized)
540  bool generalized = true;
541  bool isherm = var.IsHermitianProblem();
542  if (var.DiagonalMass() || var.UseCholeskyFactoForMass())
543  generalized = false;
544  else if (var.GetComputationalMode() == var.INVERT_MODE)
545  {
546  generalized = false;
547  isherm = false;
548  }
549 
550  if (generalized)
551  EPSSetOperators(eps, stiff, mass);
552  else
553  EPSSetOperators(eps, stiff, NULL);
554 
555  if (isherm)
556  {
557  if (generalized)
558  EPSSetProblemType(eps, EPS_GHEP);
559  else
560  EPSSetProblemType(eps, EPS_HEP);
561  }
562  else
563  {
564  if (generalized)
565  EPSSetProblemType(eps, EPS_PGNHEP);
566  else
567  EPSSetProblemType(eps, EPS_NHEP);
568  }
569 
570  EPSWhich which(EPS_LARGEST_MAGNITUDE);
571  switch (var.GetTypeSorting())
572  {
573  case EigenProblem_Base<PetscScalar>::SORTED_REAL : which = EPS_LARGEST_REAL; break;
574  case EigenProblem_Base<PetscScalar>::SORTED_IMAG : which = EPS_LARGEST_IMAGINARY; break;
575  case EigenProblem_Base<PetscScalar>::SORTED_MODULUS : which = EPS_LARGEST_MAGNITUDE; break;
576  case EigenProblem_Base<PetscScalar>::SORTED_USER : which = EPS_WHICH_USER; break;
577  }
578 
579  if (var.GetTypeSpectrum() == var.SMALL_EIGENVALUES)
580  {
581  switch (var.GetTypeSorting())
582  {
583  case EigenProblem_Base<PetscScalar>::SORTED_REAL : which = EPS_SMALLEST_REAL; break;
584  case EigenProblem_Base<PetscScalar>::SORTED_IMAG : which = EPS_SMALLEST_IMAGINARY; break;
585  case EigenProblem_Base<PetscScalar>::SORTED_MODULUS : which = EPS_SMALLEST_MAGNITUDE; break;
586  }
587  }
588 
589  if ((var.GetComputationalMode() == var.REGULAR_MODE) &&
590  (var.GetTypeSpectrum() == var.CENTERED_EIGENVALUES))
591  {
592  PetscScalar target = shiftr;
593  switch (var.GetTypeSorting())
594  {
595  case EigenProblem_Base<PetscScalar>::SORTED_REAL : which = EPS_TARGET_REAL; break;
596  case EigenProblem_Base<PetscScalar>::SORTED_IMAG : which = EPS_TARGET_IMAGINARY; target = shifti; break;
597  case EigenProblem_Base<PetscScalar>::SORTED_MODULUS : which = EPS_TARGET_MAGNITUDE; break;
598  }
599 
600  EPSSetTarget(eps, target);
601  }
602 
603  EPSSetWhichEigenpairs(eps, which);
604  if (which == EPS_WHICH_USER)
605  EPSSetEigenvalueComparison(eps, &EigenProblem_Base<PetscScalar>::GetComparisonEigenvalueSlepc, &var);
606 
607  double tol = var.GetStoppingCriterion();
608  int nb_max_iter = var.GetNbMaximumIterations();
609  EPSSetTolerances(eps, tol, nb_max_iter);
610  EPSSetDimensions(eps, nev, ncv, PETSC_DEFAULT);
611  EPSSetFromOptions(eps);
612 
613  // computing needed operators
614  if (var.DiagonalMass() || var.UseCholeskyFactoForMass())
615  {
616  if (var.DiagonalMass())
617  {
618  // computation of M
619  var.ComputeDiagonalMass();
620 
621  // computation of M^{-1/2}
622  var.FactorizeDiagonalMass();
623  }
624  else
625  {
626  // computation of M for Cholesky factorisation
627  var.ComputeMassForCholesky();
628 
629  // computation of Cholesky factorisation M = L L^T
630  var.FactorizeCholeskyMass();
631  }
632 
633  if (var.GetComputationalMode() != var.REGULAR_MODE)
634  var.ComputeAndFactorizeStiffnessMatrix(-shiftr, one);
635  }
636  else
637  {
638  if (var.GetComputationalMode() == var.INVERT_MODE)
639  {
640  if (var.GetTypeSpectrum() != var.CENTERED_EIGENVALUES)
641  {
642  // computation and factorisation of mass matrix
643  var.ComputeAndFactorizeStiffnessMatrix(one, zero);
644 
645  // computation of stiffness matrix
646  var.ComputeStiffnessMatrix();
647  }
648  else
649  {
650  // computation and factorization of (K - sigma M)
651  var.ComputeAndFactorizeStiffnessMatrix(-shiftr, one);
652 
653  // computation of M
654  var.ComputeMassMatrix();
655  }
656  }
657  else if (var.GetComputationalMode() == var.REGULAR_MODE)
658  {
659  // factorization of the mass matrix
660  var.ComputeAndFactorizeStiffnessMatrix(one, zero);
661 
662  // computation of stiffness and mass matrix
663  var.ComputeStiffnessMatrix();
664  var.ComputeMassMatrix();
665 
666  cout << "Case not implemented" << endl;
667  abort();
668  }
669  else
670  {
671  cout << "Case not implemented" << endl;
672  abort();
673  }
674  }
675 
676  // the eigenvalue problem is solved
677  int ierr = EPSSolve(eps);
678  if (ierr != 0)
679  {
680  cout << "Error during solution of eigensystem = " << ierr << endl;
681  abort();
682  }
683 
684  if (print_level >= 2)
685  {
686  PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
687  EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD);
688  PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
689  }
690 
691  EPSConvergedReason reason;
692  EPSGetConvergedReason(eps, &reason);
693  if (reason < 0)
694  {
695  if (rank_proc == 0)
696  cout << "The solver did not converge" << endl;
697 
698  if (reason == EPS_DIVERGED_ITS)
699  throw SolverMaximumIterationError("FindEigenvaluesSlepc", "Maximum number of iterations reached");
700  else
701  throw SolverDivergenceError("FindEigenvaluesSlepc", "The solver diverged");
702  }
703 
704  // eigenvalues and eigenvectors are extracted
705  Vec Vr, Vi, Vr_next, Vi_next;
706  Vector<PetscScalar> Vr_vec, Vi_vec, Vr_vec_next, Vi_vec_next;
707  AllocatePetscVector(var.GetCommunicator(), Vr, n, nglob, Vr_vec);
708  AllocatePetscVector(var.GetCommunicator(), Vi, n, nglob, Vi_vec);
709  AllocatePetscVector(var.GetCommunicator(), Vr_next, n, nglob, Vr_vec_next);
710  AllocatePetscVector(var.GetCommunicator(), Vi_next, n, nglob, Vi_vec_next);
711 
712  eigen_values.Reallocate(nev);
713  lambda_imag.Reallocate(nev);
714  eigen_vectors.Reallocate(n, nev);
715  int num = 0;
716  PetscScalar Lr, Li, Lr_next, Li_next;
717  bool eigen_pair = true;
718  while (num < nev)
719  {
720  if (eigen_pair)
721  EPSGetEigenpair(eps, num, &Lr, &Li, Vr, Vi);
722 
723  if (num < nev-1)
724  EPSGetEigenpair(eps, num+1, &Lr_next, &Li_next, Vr_next, Vi_next);
725 
726  eigen_pair = PutEigenpairLapackForm(num, nev, Lr, Li, Vr_vec, Vi_vec,
727  Lr_next, Li_next, Vr_vec_next, Vi_vec_next,
728  eigen_values, lambda_imag, eigen_vectors);
729 
730  if (eigen_pair)
731  num += 2;
732  else
733  {
734  Lr = Lr_next; Li = Li_next;
735  Vr_vec = Vr_vec_next; Vi_vec = Vi_vec_next;
736  num++;
737  }
738  }
739 
740  Vr_vec.Nullify();
741  Vi_vec.Nullify();
742  Vr_vec_next.Nullify();
743  Vi_vec_next.Nullify();
744 
745  // modifies eigenvalues and eigenvectors if needed
746  ApplyScalingEigenvec(var, eigen_values, lambda_imag, eigen_vectors,
747  shiftr, shifti);
748 
749  // temporary objects are destroyed
750  VecDestroy(&Vr);
751  VecDestroy(&Vi);
752  EPSDestroy(&eps);
753  MatDestroy(&stiff);
754  MatDestroy(&mass);
755 
756  // clears eigenproblem
757  var.Clear();
758 
759  }
760 
761  /*****************************
762  * Interface with PEP solver *
763  *****************************/
764 
765 
766  void ApplySpectralTransform(double shift, Vector<double>& eigen_values, Vector<double>& lambda_imag)
767  {
768  for (int i = 0; i < eigen_values.GetM(); i++)
769  {
770  if (lambda_imag(i) == 0.0)
771  eigen_values(i) = shift + 1.0/eigen_values(i);
772  else
773  {
774  complex<double> val(eigen_values(i), lambda_imag(i));
775  complex<double> z = shift + 1.0/val;
776  eigen_values(i) = real(z); lambda_imag(i) = imag(z);
777  }
778  }
779  }
780 
781  void ApplySpectralTransform(complex<double> shift, Vector<complex<double> >& eigen_values, Vector<complex<double> >& lambda_imag)
782  {
783  for (int i = 0; i < eigen_values.GetM(); i++)
784  eigen_values(i) = shift + 1.0/eigen_values(i);
785  }
786 
787 
788  void SetParametersSlepc(const SlepcParamPep& param, PEP& pep)
789  {
790  switch(param.GetEigensolverType())
791  {
792  case SlepcParamPep::TOAR : PEPSetType(pep, PEPTOAR); break;
793  case SlepcParamPep::STOAR : PEPSetType(pep, PEPSTOAR); break;
794  case SlepcParamPep::QARNOLDI : PEPSetType(pep, PEPQARNOLDI); break;
795  case SlepcParamPep::LINEAR : PEPSetType(pep, PEPLINEAR); break;
796  case SlepcParamPep::JD : PEPSetType(pep, PEPJD); break;
797  }
798  }
799 
800  struct MySlepcOperator
801  {
802  PolynomialEigenProblem_Base<PetscScalar>* var;
803  int num_op;
804  };
805 
807  PetscErrorCode MatMult_PepOperator(Mat mat, Vec x, Vec y)
808  {
809  Vector<PetscScalar> xvec, yvec;
810  CopyPointerPetsc(x, xvec);
811  CopyPointerPetsc(y, yvec);
812 
813  void* ctx;
814  MatShellGetContext(mat, &ctx);
815 
816  MySlepcOperator& op
817  = *reinterpret_cast<MySlepcOperator*>(ctx);
818 
819  op.var->MltOperator(op.num_op, SeldonNoTrans, xvec, yvec);
820 
821  xvec.Nullify(); yvec.Nullify();
822  return 0;
823  }
824 
826  PetscErrorCode MatMult_PepOperatorTrans(Mat mat, Vec x, Vec y)
827  {
828  Vector<PetscScalar> xvec, yvec;
829  CopyPointerPetsc(x, xvec);
830  CopyPointerPetsc(y, yvec);
831 
832  void* ctx;
833  MatShellGetContext(mat, &ctx);
834 
835  MySlepcOperator& op
836  = *reinterpret_cast<MySlepcOperator*>(ctx);
837 
838  op.var->MltOperator(op.num_op, SeldonTrans, xvec, yvec);
839 
840  xvec.Nullify(); yvec.Nullify();
841  return 0;
842  }
843 
845  PetscErrorCode MatSolve_PepOperator(Mat mat, Vec x, Vec y)
846  {
847  Vector<PetscScalar> xvec, yvec;
848  CopyPointerPetsc(x, xvec);
849  CopyPointerPetsc(y, yvec);
850 
851  void* ctx;
852  MatShellGetContext(mat, &ctx);
853 
854  MySlepcOperator& op
855  = *reinterpret_cast<MySlepcOperator*>(ctx);
856 
857  if (op.var->UseSpectralTransformation())
858  op.var->SolveOperator(SeldonNoTrans, xvec, yvec);
859  else
860  op.var->SolveMass(SeldonNoTrans, xvec, yvec);
861 
862  op.var->IncrementLinearSolves();
863 
864  xvec.Nullify(); yvec.Nullify();
865  return 0;
866  }
867 
869  PetscErrorCode MatSolveTrans_PepOperator(Mat mat, Vec x, Vec y)
870  {
871  cout << "Not implemented" << endl;
872  abort();
873  }
874 
876  void FindEigenvaluesSlepc_(PolynomialEigenProblem_Base<PetscScalar>& var,
877  Vector<PetscScalar>& eigen_values,
878  Vector<PetscScalar>& lambda_imag,
879  Matrix<PetscScalar, General, ColMajor>& eigen_vectors)
880  {
881  // initializing of computation
882  PetscScalar shift = var.GetShiftValue();
883  PetscScalar zero, one;
884  SetComplexZero(zero);
885  SetComplexOne(one);
886 
887  // degree of polynom
888  int deg_pol = var.GetPolynomialDegree();
889  int nev = var.GetNbAskedEigenvalues();
890  int n = var.GetM();
891  int nglob = var.GetGlobalM();
892 
893  // creation of shell matrices
894  Vector<Mat> op(deg_pol+1);
895  Vector<MySlepcOperator> my_op(deg_pol+1);
896  for (int k = 0; k <= deg_pol; k++)
897  {
898  my_op(k).var = &var;
899  my_op(k).num_op = k;
900  MatCreateShell(var.GetCommunicator(), n, n, nglob, nglob,
901  reinterpret_cast<void*>(&my_op(k)), &op(k));
902 
903  MatShellSetOperation(op(k), MATOP_MULT, (void(*)())MatMult_PepOperator);
904  MatShellSetOperation(op(k), MATOP_MULT_TRANSPOSE, (void(*)())MatMult_PepOperatorTrans);
905  }
906 
907  // creation of the eigensolver
908  PEP solver; ST st;
909  PEPCreate(var.GetCommunicator(), &solver);
910 
911  PEPSetOperators(solver, deg_pol+1, op.GetData());
912  PEPGetST(solver, &st);
913  st->matsolve = (PetscErrorCode(*)(Mat, Vec, Vec))MatSolve_PepOperator;
914  st->matsolve_trans = (PetscErrorCode(*)(Mat, Vec, Vec))MatSolveTrans_PepOperator;
915  //st->D = NULL;
916  STSetMatMode(st, ST_MATMODE_SHELL);
917  //STSetType(st, STSINVERT);
918 
919  // sorting and selection of spectrum
920  PEPWhich which(PEP_LARGEST_MAGNITUDE);
921  switch (var.GetTypeSorting())
922  {
923  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_REAL : which = PEP_LARGEST_REAL; break;
924  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_IMAG : which = PEP_LARGEST_IMAGINARY; break;
925  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_MODULUS : which = PEP_LARGEST_MAGNITUDE; break;
926  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_USER : which = PEP_WHICH_USER; break;
927  }
928 
929  if (var.GetTypeSpectrum() == var.SMALL_EIGENVALUES)
930  {
931  switch (var.GetTypeSorting())
932  {
933  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_REAL : which = PEP_SMALLEST_REAL; break;
934  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_IMAG : which = PEP_SMALLEST_IMAGINARY; break;
935  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_MODULUS : which = PEP_SMALLEST_MAGNITUDE; break;
936  }
937  }
938 
939  if ((!var.UseSpectralTransformation()) &&
940  (var.GetTypeSpectrum() == var.CENTERED_EIGENVALUES))
941  {
942  switch (var.GetTypeSorting())
943  {
944  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_REAL : which = PEP_TARGET_REAL; break;
945  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_IMAG : which = PEP_TARGET_IMAGINARY; break;
946  case PolynomialEigenProblem_Base<PetscScalar>::SORTED_MODULUS : which = PEP_TARGET_MAGNITUDE; break;
947  }
948 
949  PEPSetTarget(solver, shift);
950  }
951 
952  PEPSetWhichEigenpairs(solver, which);
953  if (which == PEP_WHICH_USER)
954  PEPSetEigenvalueComparison(solver, &PolynomialEigenProblem_Base<PetscScalar>::GetComparisonEigenvalueSlepc, &var);
955 
956  // tolerance and number of iterations
957  double tol = var.GetStoppingCriterion();
958  int nb_max_iter = var.GetNbMaximumIterations();
959  PEPSetTolerances(solver, tol, nb_max_iter);
960  PEPSetDimensions(solver, nev, PETSC_DECIDE, PETSC_DECIDE);
961 
962  if (!var.UseSpectralTransformation())
963  var.FactorizeMass();
964  else
965  {
966  // binomial coefficients are computed
967  Matrix<int> binom_coef(deg_pol+1, deg_pol+1);
968  binom_coef.Zero();
969  binom_coef(0, 0) = 1;
970  binom_coef(1, 0) = 1; binom_coef(1, 1) = 1;
971  for (int n = 2; n <= deg_pol; n++)
972  {
973  binom_coef(n, 0) = 1; binom_coef(n, n) = 1;
974  for (int k = 1; k < n; k++)
975  binom_coef(n, k) = binom_coef(n-1, k-1) + binom_coef(n-1, k);
976  }
977 
978  // setting operators to compute
979  Vector<PetscScalar> coef(deg_pol+1);
980  for (int k = 0; k <= deg_pol; k++)
981  {
982  coef.Zero(); PetscScalar pow_shift = 1.0;
983  for (int j = 0; j <= deg_pol-k; j++)
984  {
985  coef(j+k) = double(binom_coef(j+k, k))*pow_shift;
986  pow_shift *= shift;
987  }
988 
989  var.ComputeOperator(deg_pol-k, coef);
990  if (k == 0)
991  var.FactorizeOperator(coef);
992  }
993  }
994 
995  PEPSetScale(solver, PEP_SCALE_SCALAR, PETSC_DECIDE, PETSC_NULL, PETSC_NULL,
996  PETSC_DECIDE, PETSC_DECIDE);
997 
998  // other parameters of PEP
999  SetParametersSlepc(var.GetSlepcParameters(), solver);
1000 
1001  // the eigenvalue problem is solved
1002  int ierr = PEPSolve(solver);
1003  if (ierr != 0)
1004  {
1005  cout << "Error during solution of eigensystem = " << ierr << endl;
1006  abort();
1007  }
1008 
1009  PEPConvergedReason reason;
1010  PEPGetConvergedReason(solver, &reason);
1011  if (reason < 0)
1012  {
1013  cout << "Failed to converged " << reason << endl;
1014  abort();
1015  }
1016 
1017  int print_level = var.GetPrintLevel();
1018  if (print_level >= 4)
1019  {
1020  PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
1021  PEPErrorView(solver, PEP_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD);
1022  PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
1023  }
1024 
1025  int rank_proc = 0;
1026 #ifdef SELDON_WITH_MPI
1027  MPI_Comm_rank(var.GetCommunicator(), &rank_proc);
1028 #endif
1029 
1030  if ((print_level >= 1) && (rank_proc == 0))
1031  cout << "Number of linear solves = " << var.GetNbLinearSolves() << endl;
1032 
1033  // eigenvalues and eigenvectors are extracted
1034  Vec Vr, Vi, Vr_next, Vi_next;
1035  Vector<PetscScalar> Vr_vec, Vi_vec, Vr_vec_next, Vi_vec_next;
1036  AllocatePetscVector(var.GetCommunicator(), Vr, n, nglob, Vr_vec);
1037  AllocatePetscVector(var.GetCommunicator(), Vi, n, nglob, Vi_vec);
1038  AllocatePetscVector(var.GetCommunicator(), Vr_next, n, nglob, Vr_vec_next);
1039  AllocatePetscVector(var.GetCommunicator(), Vi_next, n, nglob, Vi_vec_next);
1040 
1041  eigen_values.Reallocate(nev);
1042  lambda_imag.Reallocate(nev);
1043  eigen_vectors.Reallocate(n, nev);
1044  int num = 0;
1045  PetscScalar Lr, Li, Lr_next, Li_next;
1046  bool eigen_pair = true;
1047  while (num < nev)
1048  {
1049  if (eigen_pair)
1050  PEPGetEigenpair(solver, num, &Lr, &Li, Vr, Vi);
1051 
1052  if (num < nev-1)
1053  PEPGetEigenpair(solver, num+1, &Lr_next, &Li_next, Vr_next, Vi_next);
1054 
1055  eigen_pair = PutEigenpairLapackForm(num, nev, Lr, Li, Vr_vec, Vi_vec,
1056  Lr_next, Li_next, Vr_vec_next, Vi_vec_next,
1057  eigen_values, lambda_imag, eigen_vectors);
1058 
1059  if (eigen_pair)
1060  num += 2;
1061  else
1062  {
1063  Lr = Lr_next; Li = Li_next;
1064  Vr_vec = Vr_vec_next; Vi_vec = Vi_vec_next;
1065  num++;
1066  }
1067  }
1068 
1069  Vr_vec.Nullify();
1070  Vi_vec.Nullify();
1071  Vr_vec_next.Nullify();
1072  Vi_vec_next.Nullify();
1073 
1074  if (var.UseSpectralTransformation())
1075  ApplySpectralTransform(shift, eigen_values, lambda_imag);
1076 
1077  // temporary objects are destroyed
1078  VecDestroy(&Vr);
1079  VecDestroy(&Vi);
1080  PEPDestroy(&solver);
1081  for (int k = 0; k <= deg_pol; k++)
1082  MatDestroy(&op(k));
1083 
1084  var.DistributeEigenvectors(eigen_vectors);
1085  }
1086 
1087 
1088  /*****************************
1089  * Interface with NEP solver *
1090  *****************************/
1091 
1093  struct NepSlepcOperator
1094  {
1095  // if true, the jacobian T'(L) is required
1096  bool jacobian;
1097  PetscScalar L; // value of L
1098  int num_operator; // for split formulation (operator number)
1099  NonLinearEigenProblem_Base<PetscScalar>* var; // pointer to the non-linear problem
1100 
1101  NepSlepcOperator()
1102  { jacobian = false; num_operator = -1; var = NULL; SetComplexZero(L); }
1103 
1104  };
1105 
1107  struct NepSlepcPreconditioning
1108  {
1109  PetscScalar L; // value of L
1110  NonLinearEigenProblem_Base<PetscScalar>* var; // pointer to the non-linear problem
1111  int type_solver; // eigensolver to be used
1112  Vector<PetscScalar> Ltab, coef_tab; // list of values for L and related coefficients
1113  Vector<int> numL_tab; // operator numbers
1114 
1115  NepSlepcPreconditioning()
1116  { SetComplexZero(L); type_solver = 0; }
1117  };
1118 
1120  void UpdatePrecond_Nep(PC pc, bool first_construct)
1121  {
1122  NepSlepcPreconditioning* shell;
1123 
1124  PCShellGetContext(pc, (void**)&shell);
1125 
1126  Mat Amat, Pmat;
1127  PCGetOperators(pc, &Amat, &Pmat);
1128 
1129  MatType type;
1130  MatGetType(Amat, &type);
1131 
1132  if (strcmp(type,MATSHELL) == 0)
1133  {
1134  // shell-matrix
1135  void* ctxF;
1136  MatShellGetContext(Amat, &ctxF);
1137 
1138  PetscScalar L;
1139  if (shell->type_solver == SlepcParamNep::NLEIGS)
1140  {
1141  NEP_NLEIGS_MATSHELL* ctx
1142  = reinterpret_cast<NEP_NLEIGS_MATSHELL*>(ctxF);
1143 
1144  // nleigs solver => we retrieve the values of L and related coefficients
1145  Vector<PetscScalar> L(ctx->nmat), coef(ctx->nmat);
1146  Vector<int> numL(ctx->nmat);
1147  for (int i = 0; i < ctx->nmat; i++)
1148  {
1149  MatShellGetContext(ctx->A[i], &ctxF);
1150 
1151  NepSlepcOperator& op
1152  = *reinterpret_cast<NepSlepcOperator*>(ctxF);
1153 
1154  L(i) = op.L; coef(i) = ctx->coeff[i];
1155  numL(i) = op.num_operator;
1156  shell->var->CheckValueL(L(i));
1157  }
1158 
1159  // new_coef will be true if the preconditioning changed
1160  bool new_coef = false;
1161  if (ctx->nmat != shell->Ltab.GetM())
1162  new_coef = true;
1163  else
1164  {
1165  for (int i = 0; i < ctx->nmat; i++)
1166  if ((abs(shell->Ltab(i) - L(i)) > 1e-12) || (abs(shell->coef_tab(i) - coef(i)) > 1e-12) || (shell->numL_tab(i) != numL(i)))
1167  new_coef = true;
1168  }
1169 
1170  // preconditioning is recomputed if needed
1171  if (new_coef || first_construct)
1172  {
1173  //cout << " Value of L for preconditioning = " << L << endl;
1174  shell->Ltab = L; shell->coef_tab = coef; shell->numL_tab = numL;
1175  if (shell->var->UseSplitMatrices())
1176  shell->var->ComputeSplitPreconditioning(numL, coef);
1177  else
1178  shell->var->ComputePreconditioning(L, coef);
1179  }
1180  }
1181  else
1182  {
1183  // other solvers => a single value of L
1184  NepSlepcOperator& op
1185  = *reinterpret_cast<NepSlepcOperator*>(ctxF);
1186 
1187  L = op.L;
1188  shell->var->CheckValueL(L);
1189 
1190  //cout << " Value of L for preconditioning = " << L << endl;
1191  // preconditioning is recomputed only if coefficients changed
1192  if ((abs(L - shell->L) >= 1e-12) || first_construct)
1193  {
1194  shell->L = L;
1195  shell->var->ComputePreconditioning(L);
1196  }
1197  }
1198  }
1199  else
1200  {
1201  // case where the matrix is stored
1202  PetscInt m, n;
1203  MatGetSize(Amat, &m, &n);
1204 
1205  Mat_SeqAIJ* A = (Mat_SeqAIJ*) Amat->data;
1206 
1207  // the matrix is duplicated
1208  DistributedMatrix<PetscScalar, General, ArrayRowSparse> Ad;
1209  Ad.Reallocate(m, n);
1210 
1211  for (int i = 0; i < m; i++)
1212  {
1213  int nz = A->i[i+1] - A->i[i];
1214  Ad.ReallocateRow(i, nz);
1215  for (int j = 0; j < nz; j++)
1216  {
1217  Ad.Index(i, j) = A->j[A->i[i]+j];
1218  Ad.Value(i, j) = A->a[A->i[i]+j];
1219  }
1220  }
1221 
1222  // and preconditioning computed
1223  shell->var->ComputeExplicitPreconditioning(Ad);
1224  }
1225  }
1226 
1227  // first construction of preconditioning
1228  PetscErrorCode ConstructPreconditioning_NepOperator(PC pc)
1229  {
1230  UpdatePrecond_Nep(pc, true);
1231 
1232  return 0;
1233  }
1234 
1235  // preconditioning may be updated (if L is different)
1236  PetscErrorCode UpdatePreconditioning_NepOperator(PC pc, KSP ksp, Vec b, Vec x)
1237  {
1238  UpdatePrecond_Nep(pc, false);
1239  return 0;
1240  }
1241 
1242  // applies the preconditioning
1243  PetscErrorCode ApplyPreconditioning_NepOperator(PC pc, Vec x, Vec y)
1244  {
1245  Vector<PetscScalar> xvec, yvec;
1246  CopyPointerPetsc(x, xvec);
1247  CopyPointerPetsc(y, yvec);
1248 
1249  NepSlepcPreconditioning* shell;
1250 
1251  PCShellGetContext(pc, (void**)&shell);
1252  shell->var->IncrementLinearSolves();
1253  shell->var->ApplyPreconditioning(SeldonNoTrans, xvec, yvec);
1254 
1255  xvec.Nullify(); yvec.Nullify();
1256  return 0;
1257  }
1258 
1259  // provides the singular points of T(L)
1260  PetscErrorCode SetSingularities_NepOperator(NEP nep, PetscInt *maxnp,
1261  PetscScalar *xi, void *ctx)
1262  {
1263  NonLinearEigenProblem_Base<PetscScalar>& var
1264  = *reinterpret_cast<NonLinearEigenProblem_Base<PetscScalar>* >(ctx);
1265 
1266  const Vector<PetscScalar>& s = var.GetSingularities();
1267  *maxnp = min(*maxnp, PetscInt(s.GetM()));
1268  for (int k = 0; k < *maxnp; k++)
1269  xi[k] = s(k);
1270 
1271  return 0;
1272  }
1273 
1274  // sets other parameters of slepc for nep
1275  void SetParametersSlepc(const SlepcParamNep& param, NEP& nep,
1276  NonLinearEigenProblem_Base<PetscScalar>& var,
1277  NepSlepcPreconditioning& prec)
1278  {
1279  switch(param.GetEigensolverType())
1280  {
1281  case SlepcParamNep::RII : NEPSetType(nep, NEPRII); break;
1282  case SlepcParamNep::SLP : NEPSetType(nep, NEPSLP); break;
1283  case SlepcParamNep::NARNOLDI : NEPSetType(nep, NEPNARNOLDI); break;
1284  case SlepcParamNep::CISS : NEPSetType(nep, NEPCISS); break;
1285  case SlepcParamNep::INTERPOL : NEPSetType(nep, NEPINTERPOL); break;
1286  case SlepcParamNep::NLEIGS : NEPSetType(nep, NEPNLEIGS); break;
1287  }
1288 
1289  KSP ksp; PC pc; bool precond = false;
1290  if (param.GetEigensolverType() == SlepcParamNep::RII)
1291  {
1292  precond = true;
1293  NEPRIISetLagPreconditioner(nep, 0);
1294  NEPRIIGetKSP(nep, &ksp);
1295  if (var.IsHermitianProblem())
1296  NEPRIISetHermitian(nep, PETSC_TRUE);
1297  }
1298  else if (param.GetEigensolverType() == SlepcParamNep::SLP)
1299  {
1300  precond = true;
1301  NEPSLPGetKSP(nep, &ksp);
1302  }
1303  else if (param.GetEigensolverType() == SlepcParamNep::NARNOLDI)
1304  {
1305  precond = true;
1306  NEPNArnoldiSetLagPreconditioner(nep, 0);
1307  NEPNArnoldiGetKSP(nep, &ksp);
1308  NEPNArnoldiSetLagPreconditioner(nep, 0);
1309  }
1310  else if (param.GetEigensolverType() == SlepcParamNep::NLEIGS)
1311  {
1312  precond = true;
1313  if (!var.UseSplitMatrices())
1314  NEPNLEIGSSetSingularitiesFunction(nep, SetSingularities_NepOperator, reinterpret_cast<void*>(&var));
1315 
1316  if (param.GetRKShifts().GetM() > 0)
1317  {
1318  int ns = param.GetRKShifts().GetM();
1319  PetscScalar shifts[ns];
1320  for (int k = 0; k < ns; k++)
1321  shifts[k] = param.GetRKShifts()(k);
1322 
1323  NEPNLEIGSSetRKShifts(nep, ns, shifts);
1324  }
1325 
1326  int nsolve; KSP* vec_ksp;
1327  NEPNLEIGSGetKSPs(nep, &nsolve, &vec_ksp);
1328  ksp = vec_ksp[0];
1329  if (nsolve > 1)
1330  {
1331  cout << "Too many solves required in nleigs = " << nsolve << endl;
1332  cout << "Not implemented" << endl;
1333  abort();
1334  }
1335 
1336  RG rg;
1337  NEPGetRG(nep, &rg);
1338  RGSetType(rg, RGINTERVAL);
1339  RGIntervalSetEndpoints(rg, param.GetLrMin(), param.GetLrMax(),
1340  param.GetLiMin(), param.GetLiMax());
1341 
1342  if (!param.InsideRegion(var.GetShiftValue()))
1343  {
1344  cout << "Target " << var.GetShiftValue() << " not inside the region"
1345  << " [ " << param.GetLrMin() << ", " << param.GetLrMax() << "] x ["
1346  << param.GetLiMin() << ", " << param.GetLiMax() << "]" << endl;
1347 
1348  abort();
1349  }
1350 
1351  if (param.FullBasis())
1352  NEPNLEIGSSetFullBasis(nep, PETSC_TRUE);
1353  else
1354  NEPNLEIGSSetFullBasis(nep, PETSC_FALSE);
1355 
1356  if (param.LockingVariant())
1357  NEPNLEIGSSetLocking(nep, PETSC_TRUE);
1358  else
1359  NEPNLEIGSSetLocking(nep, PETSC_FALSE);
1360 
1361  NEPNLEIGSSetRestart(nep, param.GetRestartNleigs());
1362 
1363  if ((param.GetInterpolationDegree() >= 0) && (param.GetInterpolationTolerance() > 0))
1364  NEPNLEIGSSetInterpolation(nep, param.GetInterpolationTolerance(),
1365  param.GetInterpolationDegree());
1366  }
1367  else
1368  {
1369  cout << "Solver " << param.GetEigensolverType() << " not interfaced in nep" << endl;
1370  abort();
1371  }
1372 
1373  if (var.ExplicitMatrix() && param.UseDefaultPetscSolver())
1374  precond = false;
1375 
1376  if (precond)
1377  {
1378  prec.type_solver = param.GetEigensolverType();
1379  prec.var = &var;
1380  if (var.ExactPreconditioning())
1381  KSPSetType(ksp, KSPPREONLY); // (dans le cas ou le preconditionneur est exact)
1382  else
1383  KSPSetType(ksp, KSPBCGS);
1384 
1385  KSPGetPC(ksp, &pc);
1386  PCSetType(pc, PCSHELL);
1387  PCShellSetContext(pc, reinterpret_cast<void*>(&prec));
1388  PCShellSetApply(pc, ApplyPreconditioning_NepOperator);
1389  //PCShellSetSetUp(pc, ConstructPreconditioning_NepOperator);
1390  PCShellSetPreSolve(pc, UpdatePreconditioning_NepOperator);
1391  }
1392 
1393  // options from the command line
1394  if (param.UseCommandLineOptions())
1395  NEPSetFromOptions(nep);
1396  }
1397 
1398  // computation of operator T(lambda)
1399  PetscErrorCode FormFunctionNEP(NEP nep, PetscScalar lambda, Mat A, Mat B, void* ctx)
1400  {
1401  void* ctxF;
1402  MatShellGetContext(A, &ctxF);
1403 
1404  NepSlepcOperator& op
1405  = *reinterpret_cast<NepSlepcOperator*>(ctxF);
1406 
1407  op.L = lambda;
1408  op.var->CheckValueL(lambda);
1409  op.var->ComputeOperator(op.L);
1410 
1411  return 0;
1412  }
1413 
1414  // computation of operator T'(lambda)
1415  PetscErrorCode FormJacobianNEP(NEP nep, PetscScalar lambda, Mat A, void* ctx)
1416  {
1417  void* ctxF;
1418  MatShellGetContext(A, &ctxF);
1419 
1420  NepSlepcOperator& op
1421  = *reinterpret_cast<NepSlepcOperator*>(ctxF);
1422 
1423  op.L = lambda;
1424  op.var->CheckValueL(lambda);
1425  op.var->ComputeJacobian(op.L);
1426 
1427  return 0;
1428  }
1429 
1430  // conversion from a Seldon matrix to a PetscMatrix
1431  void ConvertToPetsc(DistributedMatrix<PetscScalar, General, ArrayRowSparse>& A, Mat Ap)
1432  {
1433  Vector<int> nnz(A.GetM());
1434  for (int i = 0; i < A.GetM(); i++)
1435  nnz(i) = A.GetRowSize(i);
1436 
1437  MatSeqAIJSetPreallocation(Ap, PETSC_DEFAULT, nnz.GetData());
1438  MatSetUp(Ap);
1439 
1440  for (int i = 0; i < A.GetM(); i++)
1441  MatSetValues(Ap, 1, &i, A.GetRowSize(i), A.GetIndex(i), A.GetData(i), INSERT_VALUES);
1442 
1443  MatAssemblyBegin(Ap, MAT_FINAL_ASSEMBLY);
1444  MatAssemblyEnd(Ap, MAT_FINAL_ASSEMBLY);
1445  }
1446 
1447  // computes T(L) explicitely
1448  PetscErrorCode FormFunctionExpNEP(NEP nep, PetscScalar lambda, Mat Ap, Mat Bp, void* ctx)
1449  {
1450  NepSlepcOperator& op
1451  = *reinterpret_cast<NepSlepcOperator*>(ctx);
1452 
1453  DistributedMatrix<PetscScalar, General, ArrayRowSparse> A;
1454  op.var->CheckValueL(lambda);
1455  op.var->ComputeOperatorExplicit(lambda, A);
1456 
1457  ConvertToPetsc(A, Ap);
1458 
1459  return 0;
1460  }
1461 
1462  // computes T'(L) explicitely
1463  PetscErrorCode FormJacobianExpNEP(NEP nep, PetscScalar lambda, Mat Ap, void* ctx)
1464  {
1465  NepSlepcOperator& op
1466  = *reinterpret_cast<NepSlepcOperator*>(ctx);
1467 
1468  DistributedMatrix<PetscScalar, General, ArrayRowSparse> A;
1469  op.var->CheckValueL(lambda);
1470  op.var->ComputeJacobianExplicit(lambda, A);
1471 
1472  ConvertToPetsc(A, Ap);
1473 
1474  return 0;
1475  }
1476 
1478  PetscErrorCode MatMult_NepOperator(Mat mat, Vec x, Vec y)
1479  {
1480  Vector<PetscScalar> xvec, yvec;
1481  CopyPointerPetsc(x, xvec);
1482  CopyPointerPetsc(y, yvec);
1483 
1484  void* ctx;
1485  MatShellGetContext(mat, &ctx);
1486 
1487  NepSlepcOperator& op
1488  = *reinterpret_cast<NepSlepcOperator*>(ctx);
1489 
1490  if (op.num_operator == -1)
1491  {
1492  op.var->CheckValueL(op.L);
1493  op.var->IncrementProdMatVect();
1494  if (op.jacobian)
1495  op.var->MltJacobian(op.L, SeldonNoTrans, xvec, yvec);
1496  else
1497  op.var->MltOperator(op.L, SeldonNoTrans, xvec, yvec);
1498  }
1499  else
1500  {
1501  op.var->IncrementProdMatVect();
1502  op.var->MltOperatorSplit(op.num_operator, SeldonNoTrans, xvec, yvec);
1503  }
1504 
1505  xvec.Nullify(); yvec.Nullify();
1506  return 0;
1507  }
1508 
1510  PetscErrorCode MatMultTrans_NepOperator(Mat mat, Vec x, Vec y)
1511  {
1512  Vector<PetscScalar> xvec, yvec;
1513  CopyPointerPetsc(x, xvec);
1514  CopyPointerPetsc(y, yvec);
1515 
1516  void* ctx;
1517  MatShellGetContext(mat, &ctx);
1518 
1519  NepSlepcOperator& op
1520  = *reinterpret_cast<NepSlepcOperator*>(ctx);
1521 
1522  if (op.num_operator == -1)
1523  {
1524  op.var->IncrementProdMatVect();
1525  op.var->CheckValueL(op.L);
1526  if (op.jacobian)
1527  op.var->MltJacobian(op.L, SeldonTrans, xvec, yvec);
1528  else
1529  op.var->MltOperator(op.L, SeldonTrans, xvec, yvec);
1530  }
1531  else
1532  {
1533  op.var->MltOperatorSplit(op.num_operator, SeldonTrans, xvec, yvec);
1534  op.var->IncrementProdMatVect();
1535  }
1536 
1537  xvec.Nullify(); yvec.Nullify();
1538  return 0;
1539  }
1540 
1541  // operator is duplicated
1542  PetscErrorCode MatDuplicate_NepOperator(Mat A, MatDuplicateOption op, Mat* B)
1543  {
1544  NepSlepcOperator *actx,*bctx;
1545  PetscInt n, nloc;
1546  MPI_Comm comm;
1547 
1548  MatShellGetContext(A, (void**)&actx);
1549  MatGetLocalSize(A, &nloc, NULL);
1550  MatGetSize(A, &n, NULL);
1551 
1552  bctx = new NepSlepcOperator;
1553  bctx->jacobian = actx->jacobian;
1554  bctx->L = actx->L;
1555  bctx->var = actx->var;
1556  bctx->num_operator = actx->num_operator;
1557 
1558  PetscObjectGetComm((PetscObject)A, &comm);
1559  MatCreateShell(comm, nloc, nloc, n, n,(void*)bctx, B);
1560  MatShellSetOperation(*B, MATOP_MULT, (void(*)(void))MatMult_NepOperator);
1561  MatShellSetOperation(*B, MATOP_MULT_TRANSPOSE, (void(*)(void))MatMultTrans_NepOperator);
1562  MatShellSetOperation(*B, MATOP_DUPLICATE, (void(*)(void))MatDuplicate_NepOperator);
1563 
1564  return 0;
1565  }
1566 
1568  void FindEigenvaluesSlepc_(NonLinearEigenProblem_Base<PetscScalar>& var,
1569  Vector<PetscScalar>& eigen_values,
1570  Vector<PetscScalar>& lambda_imag,
1571  Matrix<PetscScalar, General, ColMajor>& eigen_vectors)
1572  {
1573  // initializing of computation
1574  PetscScalar shift = var.GetShiftValue();
1575  PetscScalar zero, one;
1576  SetComplexZero(zero);
1577  SetComplexOne(one);
1578 
1579  // dimensions
1580  int nev = var.GetNbAskedEigenvalues();
1581  int n = var.GetM();
1582  int nglob = var.GetGlobalM();
1583 
1584  // creation of the eigensolver
1585  NEP solver;
1586  NEPCreate(var.GetCommunicator(), &solver);
1587 
1588  // creation of shell matrices
1589  Mat EvalF, EvalJacob;
1590  NepSlepcOperator operT;
1591  operT.jacobian = false; operT.var = &var;
1592  NepSlepcOperator operTp;
1593  operTp.jacobian = true; operTp.var = &var;
1594 
1595 
1596  // matrices for the split formulation
1597  Vector<Mat> EvalF_Split(var.GetNbSplitMatrices());
1598  Vector<NepSlepcOperator> operSplit(var.GetNbSplitMatrices());
1599  Vector<FN> FunctionF_Split(var.GetNbSplitMatrices());
1600 
1601  if (var.UseSplitMatrices())
1602  {
1603  if (var.ExplicitMatrix())
1604  {
1605  for (int i = 0; i < var.GetNbSplitMatrices(); i++)
1606  {
1607  MatCreate(var.GetCommunicator(), &EvalF_Split(i));
1608  MatSetSizes(EvalF_Split(i), n, n, nglob, nglob);
1609  MatSetType(EvalF_Split(i), MATSEQAIJ);
1610 
1611  DistributedMatrix<PetscScalar, General, ArrayRowSparse> A;
1612  var.ComputeOperatorSplitExplicit(i, A);
1613 
1614  ConvertToPetsc(A, EvalF_Split(i));
1615  }
1616  }
1617  else
1618  {
1619  for (int i = 0; i < var.GetNbSplitMatrices(); i++)
1620  {
1621  operSplit(i).jacobian = false; operSplit(i).var = &var;
1622  operSplit(i).num_operator = i;
1623  MatCreateShell(var.GetCommunicator(), n, n, nglob, nglob,
1624  reinterpret_cast<void*>(&operSplit(i)), &EvalF_Split(i));
1625 
1626  MatShellSetOperation(EvalF_Split(i), MATOP_MULT, (void(*)(void))MatMult_NepOperator);
1627  MatShellSetOperation(EvalF_Split(i), MATOP_MULT_TRANSPOSE, (void(*)(void))MatMultTrans_NepOperator);
1628  MatShellSetOperation(EvalF_Split(i), MATOP_DUPLICATE, (void(*)(void))MatDuplicate_NepOperator);
1629  }
1630  }
1631 
1632  for (int i = 0; i < var.GetNbSplitMatrices(); i++)
1633  {
1634  const Vector<PetscScalar>& P = var.GetNumeratorSplitFct(i);
1635  const Vector<PetscScalar>& Q = var.GetDenominatorSplitFct(i);
1636  FNCreate(var.GetCommunicator(), &FunctionF_Split(i));
1637  FNSetType(FunctionF_Split(i), FNRATIONAL);
1638  FNRationalSetNumerator(FunctionF_Split(i), P.GetM(), P.GetData());
1639  FNRationalSetDenominator(FunctionF_Split(i), Q.GetM(), Q.GetData());
1640  }
1641 
1642  NEPSetProblemType(solver, NEP_RATIONAL);
1643  NEPSetSplitOperator(solver, var.GetNbSplitMatrices(), EvalF_Split.GetData(),
1644  FunctionF_Split.GetData(), DIFFERENT_NONZERO_PATTERN);
1645  }
1646  else
1647  {
1648  if (var.ExplicitMatrix())
1649  {
1650  MatCreate(var.GetCommunicator(), &EvalF);
1651  MatSetSizes(EvalF, n, n, nglob, nglob);
1652  MatSetType(EvalF, MATSEQAIJ);
1653  MatSetUp(EvalF);
1654  }
1655  else
1656  {
1657  MatCreateShell(var.GetCommunicator(), n, n, nglob, nglob,
1658  reinterpret_cast<void*>(&operT), &EvalF);
1659 
1660  MatShellSetOperation(EvalF, MATOP_MULT, (void(*)(void))MatMult_NepOperator);
1661  MatShellSetOperation(EvalF, MATOP_MULT_TRANSPOSE, (void(*)(void))MatMultTrans_NepOperator);
1662  MatShellSetOperation(EvalF, MATOP_DUPLICATE, (void(*)(void))MatDuplicate_NepOperator);
1663  }
1664 
1665  if (var.ExplicitMatrix())
1666  {
1667  MatCreate(var.GetCommunicator(), &EvalJacob);
1668  MatSetSizes(EvalJacob, n, n, nglob, nglob);
1669  MatSetType(EvalJacob, MATSEQAIJ);
1670  MatSetUp(EvalJacob);
1671  }
1672  else
1673  {
1674  MatCreateShell(var.GetCommunicator(), n, n, nglob, nglob,
1675  reinterpret_cast<void*>(&operTp), &EvalJacob);
1676 
1677  MatShellSetOperation(EvalJacob, MATOP_MULT, (void(*)(void))MatMult_NepOperator);
1678  MatShellSetOperation(EvalJacob, MATOP_MULT_TRANSPOSE, (void(*)(void))MatMultTrans_NepOperator);
1679  MatShellSetOperation(EvalJacob, MATOP_DUPLICATE, (void(*)(void))MatDuplicate_NepOperator);
1680  }
1681 
1682  if (var.ExplicitMatrix())
1683  NEPSetFunction(solver, EvalF, EvalF, FormFunctionExpNEP, reinterpret_cast<void*>(&operT));
1684  else
1685  NEPSetFunction(solver, EvalF, EvalF, FormFunctionNEP, NULL);
1686 
1687  if (var.ExplicitMatrix())
1688  NEPSetJacobian(solver, EvalJacob, FormJacobianExpNEP, reinterpret_cast<void*>(&operTp));
1689  else
1690  NEPSetJacobian(solver, EvalJacob, FormJacobianNEP, NULL);
1691  }
1692 
1693  // sorting and selection of spectrum
1694  NEPWhich which(NEP_LARGEST_MAGNITUDE);
1695  switch (var.GetTypeSorting())
1696  {
1697  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_REAL : which = NEP_LARGEST_REAL; break;
1698  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_IMAG : which = NEP_LARGEST_IMAGINARY; break;
1699  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_MODULUS : which = NEP_LARGEST_MAGNITUDE; break;
1700  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_USER : which = NEP_WHICH_USER; break;
1701  }
1702 
1703  if (var.GetTypeSpectrum() == var.SMALL_EIGENVALUES)
1704  {
1705  switch (var.GetTypeSorting())
1706  {
1707  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_REAL : which = NEP_SMALLEST_REAL; break;
1708  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_IMAG : which = NEP_SMALLEST_IMAGINARY; break;
1709  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_MODULUS : which = NEP_SMALLEST_MAGNITUDE; break;
1710  }
1711  }
1712 
1713  if (var.GetTypeSpectrum() == var.CENTERED_EIGENVALUES)
1714  {
1715  switch (var.GetTypeSorting())
1716  {
1717  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_REAL : which = NEP_TARGET_REAL; break;
1718  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_IMAG : which = NEP_TARGET_IMAGINARY; break;
1719  case NonLinearEigenProblem_Base<PetscScalar>::SORTED_MODULUS : which = NEP_TARGET_MAGNITUDE; break;
1720  }
1721 
1722  NEPSetTarget(solver, shift);
1723  }
1724 
1725  NEPSetWhichEigenpairs(solver, which);
1726  if (which == NEP_WHICH_USER)
1727  NEPSetEigenvalueComparison(solver, &NonLinearEigenProblem_Base<PetscScalar>::GetComparisonEigenvalueSlepc, &var);
1728 
1729  // tolerance and number of iterations
1730  double tol = var.GetStoppingCriterion();
1731  int nb_max_iter = var.GetNbMaximumIterations();
1732  NEPSetTolerances(solver, tol, nb_max_iter);
1733  int ncv = PETSC_DEFAULT;
1734  if (var.GetNbArnoldiVectors() > 0)
1735  ncv = var.GetNbArnoldiVectors();
1736 
1737  if (nev > 0)
1738  NEPSetDimensions(solver, nev, ncv, PETSC_DEFAULT);
1739 
1740  // other parameters of NEP
1741  NepSlepcPreconditioning prec;
1742  SetParametersSlepc(var.GetSlepcParameters(), solver, var, prec);
1743 
1744  // the eigenvalue problem is solved
1745  PetscErrorCode ierr = NEPSolve(solver);
1746  if (ierr != 0)
1747  {
1748  cout << "Error during solution of eigensystem = " << ierr << endl;
1749  abort();
1750  }
1751 
1752  int print_level = var.GetPrintLevel();
1753  if (print_level >= 4)
1754  {
1755  PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
1756  NEPConvergedReasonView(solver, PETSC_VIEWER_STDOUT_WORLD);
1757  NEPErrorView(solver, NEP_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD);
1758  PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
1759  }
1760 
1761  // nev is modified if less eigenvalues are converged
1762  NEPGetConverged(solver, &nev);
1763 
1764  // eigenvalues and eigenvectors are extracted
1765  Vec Vr, Vi, Vr_next, Vi_next;
1766  Vector<PetscScalar> Vr_vec, Vi_vec, Vr_vec_next, Vi_vec_next;
1767  AllocatePetscVector(var.GetCommunicator(), Vr, n, nglob, Vr_vec);
1768  AllocatePetscVector(var.GetCommunicator(), Vi, n, nglob, Vi_vec);
1769  AllocatePetscVector(var.GetCommunicator(), Vr_next, n, nglob, Vr_vec_next);
1770  AllocatePetscVector(var.GetCommunicator(), Vi_next, n, nglob, Vi_vec_next);
1771 
1772  eigen_values.Reallocate(nev);
1773  lambda_imag.Reallocate(nev);
1774  eigen_vectors.Reallocate(n, nev);
1775  int num = 0;
1776  PetscScalar Lr, Li, Lr_next, Li_next;
1777  bool eigen_pair = true;
1778  while (num < nev)
1779  {
1780  if (eigen_pair)
1781  NEPGetEigenpair(solver, num, &Lr, &Li, Vr, Vi);
1782 
1783  if (num < nev-1)
1784  NEPGetEigenpair(solver, num+1, &Lr_next, &Li_next, Vr_next, Vi_next);
1785 
1786  eigen_pair = PutEigenpairLapackForm(num, nev, Lr, Li, Vr_vec, Vi_vec,
1787  Lr_next, Li_next, Vr_vec_next, Vi_vec_next,
1788  eigen_values, lambda_imag, eigen_vectors);
1789 
1790  if (eigen_pair)
1791  num += 2;
1792  else
1793  {
1794  Lr = Lr_next; Li = Li_next;
1795  Vr_vec = Vr_vec_next; Vi_vec = Vi_vec_next;
1796  num++;
1797  }
1798  }
1799 
1800  Vr_vec.Nullify();
1801  Vi_vec.Nullify();
1802  Vr_vec_next.Nullify();
1803  Vi_vec_next.Nullify();
1804 
1805  // temporary objects are destroyed
1806  VecDestroy(&Vr);
1807  VecDestroy(&Vi);
1808  NEPDestroy(&solver);
1809  if (var.UseSplitMatrices())
1810  {
1811  for (int i = 0; i < var.GetNbSplitMatrices(); i++)
1812  {
1813  MatDestroy(&EvalF_Split(i));
1814  FNDestroy(&FunctionF_Split(i));
1815  }
1816  }
1817  else
1818  {
1819  MatDestroy(&EvalF);
1820  MatDestroy(&EvalJacob);
1821  }
1822 
1823  var.DistributeEigenvectors(eigen_vectors);
1824  }
1825 
1826 #ifdef SELDON_WITH_VIRTUAL
1827  void FindEigenvaluesSlepc(EigenProblem_Base<PetscScalar>& var,
1828  Vector<PetscScalar>& eigen_values,
1829  Vector<PetscScalar>& lambda_imag,
1830  Matrix<PetscScalar, General, ColMajor>& eigen_vectors)
1831  {
1832  FindEigenvaluesSlepc_(var, eigen_values, lambda_imag, eigen_vectors);
1833  }
1834 
1835  void FindEigenvaluesSlepc(PolynomialEigenProblem_Base<Petsc_Scalar>& var,
1836  Vector<Petsc_Scalar>& eigen_values,
1837  Vector<Petsc_Scalar>& lambda_imag,
1838  Matrix<Petsc_Scalar, General, ColMajor>& eigen_vectors)
1839  {
1840  FindEigenvaluesSlepc_(var, eigen_values, lambda_imag, eigen_vectors);
1841  }
1842 
1843  void FindEigenvaluesSlepc(NonLinearEigenProblem_Base<Petsc_Scalar>& var,
1844  Vector<Petsc_Scalar>& eigen_values,
1845  Vector<Petsc_Scalar>& lambda_imag,
1846  Matrix<Petsc_Scalar, General, ColMajor>& eigen_vectors)
1847  {
1848  FindEigenvaluesSlepc_(var, eigen_values, lambda_imag, eigen_vectors);
1849  }
1850 #else
1851  template<class EigenProblem, class T, class Allocator1,
1852  class Allocator2, class Allocator3>
1853  void FindEigenvaluesSlepc(EigenProblem& var,
1854  Vector<T, VectFull, Allocator1>& eigen_values,
1855  Vector<T, VectFull, Allocator2>& lambda_imag,
1856  Matrix<T, General, ColMajor, Allocator3>& eigen_vectors)
1857  {
1858  cout << "Recompile with SELDON_WITH_VIRTUAL" << endl;
1859  abort();
1860  }
1861 #endif
1862 
1863 }
1864 
1865 #define SELDON_FILE_SLEPC_HXX
1866 #endif
Seldon::GeneralEigenProblem_Base::GetM
int GetM() const
returns number of rows
Definition: VirtualEigenvalueSolver.cxx:359
Seldon::EigenProblem_Base::GetTypeSpectrum
int GetTypeSpectrum() const
returns the spectrum desired (large, small eigenvalues, etc)
Definition: EigenvalueSolver.cxx:253
Seldon::UpdatePrecond_Nep
void UpdatePrecond_Nep(PC pc, bool first_construct)
preconditioning is updated
Definition: Slepc.cxx:1120
Seldon::MatMult_MassMatrix
PetscErrorCode MatMult_MassMatrix(Mat mat, Vec x, Vec y)
matrix-vector product y = mat x (mat : mass operator)
Definition: Slepc.cxx:410
Seldon::MatMult_Matrix
PetscErrorCode MatMult_Matrix(Mat mat, Vec x, Vec y)
matrix-vector product y = mat x (mat : stiffness operator)
Definition: Slepc.cxx:320
Seldon::CopyPointerPetsc
void CopyPointerPetsc(const Vec &x, Vector< PetscScalar > &y)
filling the vector y from pointer contained in x (Petsc vector)
Definition: Slepc.cxx:295
Seldon::PolynomialEigenProblem_Base::GetSlepcParameters
SlepcParamPep & GetSlepcParameters()
returns object storing slepc parameters
Definition: PolynomialEigenvalueSolver.cxx:55
Seldon::NepSlepcPreconditioning
class for preconditioning for NEP solver
Definition: Slepc.cxx:1107
Seldon::PolynomialEigenProblem_Base::SolveMass
virtual void SolveMass(const SeldonTranspose &, const Vector< T > &x, Vector< T > &y)
to overload for non-diagonal mass
Definition: PolynomialEigenvalueSolver.cxx:123
Seldon::Vector< PetscScalar >
Seldon::GeneralEigenProblem_Base::GetTypeSpectrum
int GetTypeSpectrum() const
returns the spectrum desired (large, small eigenvalues, etc)
Definition: VirtualEigenvalueSolver.cxx:294
Seldon::GeneralEigenProblem_Base::GetNbLinearSolves
int GetNbLinearSolves() const
returns the number of linear solves
Definition: VirtualEigenvalueSolver.cxx:408
Seldon::PolynomialEigenProblem_Base::FactorizeMass
virtual void FactorizeMass()
to overload
Definition: PolynomialEigenvalueSolver.cxx:113
Seldon::EigenProblem_Base::GetComputationalMode
int GetComputationalMode() const
returns the spectral transformation used for evaluation of eigenvalues
Definition: EigenvalueSolver.cxx:133
Seldon::MatGetDiagonal_Matrix
PetscErrorCode MatGetDiagonal_Matrix(Mat A, Vec d)
retrieves diagonal of A, d = diag(A)
Definition: Slepc.cxx:439
Seldon::EigenProblem_Base::MltStiffness
void MltStiffness(const Vector< T > &X, Vector< T > &Y)
matrix vector product with stifness matrix Y = K X
Definition: EigenvalueSolver.cxx:709
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::MatSolve_PepOperator
PetscErrorCode MatSolve_PepOperator(Mat mat, Vec x, Vec y)
solving op y = x
Definition: Slepc.cxx:845
Seldon::MatMult_NepOperator
PetscErrorCode MatMult_NepOperator(Mat mat, Vec x, Vec y)
matrix-vector product y = op x
Definition: Slepc.cxx:1478
Seldon::GeneralEigenProblem_Base::GetNbMaximumIterations
int GetNbMaximumIterations() const
returns the maximal number of iterations allowed for the iterative algorithm
Definition: VirtualEigenvalueSolver.cxx:329
Seldon::GeneralEigenProblem::GetShiftValue
T GetShiftValue() const
returns the shift value used
Definition: VirtualEigenvalueSolver.cxx:483
Seldon::IsComplexNumber
bool IsComplexNumber(const T &number)
Returns true for a complex number.
Definition: CommonInline.cxx:293
Seldon::MatMult_PepOperatorTrans
PetscErrorCode MatMult_PepOperatorTrans(Mat mat, Vec x, Vec y)
matrix-vector product y = op^T x
Definition: Slepc.cxx:826
Seldon::MatMultTranspose_Matrix
PetscErrorCode MatMultTranspose_Matrix(Mat A, Vec x, Vec y)
matrix-vector product with transpose, y = mat^T x
Definition: Slepc.cxx:431
Seldon::EigenProblem_Base::MltInvSqrtDiagonalMass
void MltInvSqrtDiagonalMass(Vector< T0 > &X)
multiplication of X by D^-1/2
Definition: EigenvalueSolver.cxx:614
Seldon::EigenProblem_Base::ComputeSolution
void ComputeSolution(const Vector< T0 > &X, Vector< T0 > &Y)
solving the linear system (a M + b K) Y = X
Definition: EigenvalueSolver.cxx:781
Seldon::EigenProblem_Base::IncrementProdMatVect
void IncrementProdMatVect()
increment of the number of matrix vector products
Definition: EigenvalueSolver.cxx:514
Seldon::EigenProblem_Base::SolveCholeskyMass
void SolveCholeskyMass(const TransStatus &transA, Vector< T > &X)
computation of L^-1 X or L^-T x if M = L L^T
Definition: EigenvalueSolver.cxx:819
Seldon::DistributedMatrix::Reallocate
void Reallocate(int m, int n)
changing the size of the local matrix, previous values are lost
Definition: DistributedMatrix.cxx:4832
Seldon::Undefined
Definition: Errors.hxx:62
Seldon::EigenProblem_Base::MltMass
void MltMass(const Vector< T > &X, Vector< T > &Y)
matrix vector product with mass matrix Y = M X
Definition: EigenvalueSolver.cxx:672
Seldon::SlepcParam::type_solver
int type_solver
which solver ?
Definition: VirtualEigenvalueSolver.hxx:68
Seldon::GeneralEigenProblem_Base::IncrementLinearSolves
void IncrementLinearSolves()
increments the number of linear solves
Definition: VirtualEigenvalueSolver.cxx:415
Seldon::EigenProblem_Base::DiagonalMass
bool DiagonalMass() const
returns true if the mass matrix is diagonal
Definition: EigenvalueSolver.cxx:412
Seldon::NepSlepcOperator
class for matshell operator for NEP solver
Definition: Slepc.cxx:1093
Seldon::GeneralEigenProblem_Base::GetPrintLevel
int GetPrintLevel() const
returns level of verbosity
Definition: VirtualEigenvalueSolver.cxx:380
Seldon::GeneralEigenProblem_Base::GetNbAskedEigenvalues
int GetNbAskedEigenvalues() const
returns the number of eigenvalues asked by the user
Definition: VirtualEigenvalueSolver.cxx:280
Seldon::EigenProblem_Base
Base class to solve an eigenvalue problem.
Definition: EigenvalueSolver.hxx:18
Seldon::MatMultTrans_NepOperator
PetscErrorCode MatMultTrans_NepOperator(Mat mat, Vec x, Vec y)
matrix-vector product y = op^T x
Definition: Slepc.cxx:1510
Seldon::PolynomialEigenProblem_Base< PetscScalar >
Seldon::DistributedMatrix
matrix distributed over all the processors
Definition: DistributedMatrix.hxx:506
Seldon::EigenProblem_Base::UseCholeskyFactoForMass
bool UseCholeskyFactoForMass() const
returns true if Cholesky factorisation has to be used for mass matrix
Definition: EigenvalueSolver.cxx:396
Seldon::PolynomialEigenProblem_Base::ComputeOperator
virtual void ComputeOperator(int num, const Vector< T > &coef)
to overload
Definition: PolynomialEigenvalueSolver.cxx:95
Seldon::GeneralEigenProblem_Base::GetStoppingCriterion
double GetStoppingCriterion() const
returns the stopping criterion
Definition: VirtualEigenvalueSolver.cxx:315
Seldon::PolynomialEigenProblem_Base::GetPolynomialDegree
int GetPolynomialDegree() const
returns the polynomial degree
Definition: PolynomialEigenvalueSolver.cxx:71
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::EigenProblem_Base::MltCholeskyMass
void MltCholeskyMass(const TransStatus &transA, Vector< T > &X)
computation of L X or L^T x if M = L L^T
Definition: EigenvalueSolver.cxx:809
Seldon::GeneralEigenProblem_Base::GetGlobalM
int GetGlobalM() const
returns global number of rows
Definition: VirtualEigenvalueSolver.cxx:366
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
Seldon::PolynomialEigenProblem_Base::UseSpectralTransformation
bool UseSpectralTransformation() const
returns true if a spectral transformation has to be used
Definition: PolynomialEigenvalueSolver.cxx:39
NEP_NLEIGS_MATSHELL
Definition: Slepc.cxx:9
Seldon::GeneralEigenProblem::DistributeEigenvectors
virtual void DistributeEigenvectors(Matrix< T, General, ColMajor > &eigen_vec)
changes final eigenvectors if needed
Definition: VirtualEigenvalueSolver.cxx:737
Seldon::MatMult_PepOperator
PetscErrorCode MatMult_PepOperator(Mat mat, Vec x, Vec y)
matrix-vector product y = op x
Definition: Slepc.cxx:807
Seldon::MySlepcOperator
Definition: Slepc.cxx:800
Seldon::PolynomialEigenProblem_Base::MltOperator
virtual void MltOperator(int num, const SeldonTranspose &, const Vector< T > &X, Vector< T > &Y)
to overload
Definition: PolynomialEigenvalueSolver.cxx:104
Seldon::MatSolveTrans_PepOperator
PetscErrorCode MatSolveTrans_PepOperator(Mat mat, Vec x, Vec y)
solving op y = x
Definition: Slepc.cxx:869
Seldon::EigenProblem_Base::MltSqrtDiagonalMass
void MltSqrtDiagonalMass(Vector< T0 > &X)
multiplication of X by D^1/2
Definition: EigenvalueSolver.cxx:624
Seldon::GeneralEigenProblem_Base::GetTypeSorting
int GetTypeSorting() const
returns how eigenvalues are sorted (real, imaginary part or modulus)
Definition: VirtualEigenvalueSolver.cxx:301