"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Kernel/EigenSolve_SLEPC.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

[ To the main getdp source changes report ]

EigenSolve_SLEPC.cpp  (getdp-3.4.0-source.tgz):EigenSolve_SLEPC.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// See the LICENSE.txt file for license information. Please report all // See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
//
// Contributor(s):
// Guillaume Demesy
#include "GetDPConfig.h" #include "GetDPConfig.h"
#if defined(HAVE_SLEPC) #if defined(HAVE_SLEPC)
// SLEPc interface for solving eigenvalue problems // SLEPc interface for solving eigenvalue problems
// //
// SLEPc can solve both linear and quadratic eigenvalue // SLEPc can solve both linear and quadratic eigenvalue
// problems. SLEPc options can be specified in the .petscrc file, or // problems. SLEPc options can be specified in the .petscrc file, or
// directly on the command line. // directly on the command line.
skipping to change at line 46 skipping to change at line 49
#include <sstream> #include <sstream>
#include <string> #include <string>
#include <cstring> #include <cstring>
#include <complex> #include <complex>
#include "ProData.h" #include "ProData.h"
#include "DofData.h" #include "DofData.h"
#include "Cal_Quantity.h" #include "Cal_Quantity.h"
#include "Message.h" #include "Message.h"
#include "MallocUtils.h" #include "MallocUtils.h"
#include <slepceps.h> #include <slepceps.h>
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5)
#include <slepcqep.h> #include <slepcqep.h>
#else #else
#include <slepcpep.h> #include <slepcpep.h>
#include <slepcnep.h> //nleigchange #include <slepcnep.h> //nleigchange
#endif #endif
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 9) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 9)
#define PCFactorSetMatSolverPackage PCFactorSetMatSolverType #define PCFactorSetMatSolverPackage PCFactorSetMatSolverType
#endif #endif
extern struct CurrentData Current ; extern struct CurrentData Current;
extern void _fillseq(Vec &V, Vec &Vseq); // from LinAlg_PETsc.cpp extern void _fillseq(Vec &V, Vec &Vseq); // from LinAlg_PETsc.cpp
static void _try(int ierr) static void _try(int ierr)
{ {
CHKERRCONTINUE(ierr); CHKERRCONTINUE(ierr);
if(PetscUnlikely(ierr)){ if(PetscUnlikely(ierr)) {
const char *text; const char *text;
PetscErrorMessage(ierr, &text, 0); PetscErrorMessage(ierr, &text, 0);
Message::Error("PETSc error: %s", text); Message::Error("PETSc error: %s", text);
Message::SetLastPETScError(ierr); Message::SetLastPETScError(ierr);
} }
} }
static PetscErrorCode _myMonitor(const char *str, int its, int nconv, PetscScala static PetscErrorCode _myMonitor(const char *str, int its, int nconv,
r *eigr, PetscScalar *eigr, PetscScalar *eigi,
PetscScalar *eigi, PetscReal* errest) PetscReal *errest)
{ {
if(!its) return 0; if(!its) return 0;
std::ostringstream sstream; std::ostringstream sstream;
sstream << " " << its << " " << str << " nconv=" << nconv << " first unconver sstream << " " << its << " " << str << " nconv=" << nconv
ged value " << " first unconverged value "
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
<< PetscRealPart(eigr[nconv]) << " + i * (" << PetscImaginaryPart(eigr << PetscRealPart(eigr[nconv]) << " + i * ("
[nconv]) << ")" << PetscImaginaryPart(eigr[nconv]) << ")"
#else #else
<< eigr[nconv] << " + i * (" << eigi[nconv] << ")" << eigr[nconv] << " + i * (" << eigi[nconv] << ")"
#endif #endif
<< " error " << errest[nconv]; << " error " << errest[nconv];
Message::Info("%s", sstream.str().c_str()); Message::Info("%s", sstream.str().c_str());
return 0; return 0;
} }
static PetscErrorCode _myEpsMonitor(EPS eps, int its, int nconv, PetscScalar *ei static PetscErrorCode _myEpsMonitor(EPS eps, int its, int nconv,
gr, PetscScalar *eigr, PetscScalar *eigi,
PetscScalar *eigi, PetscReal* errest, int ne PetscReal *errest, int nest, void *mctx)
st,
void *mctx)
{ {
return _myMonitor("EPS", its, nconv, eigr, eigi, errest); return _myMonitor("EPS", its, nconv, eigr, eigi, errest);
} }
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5)
static PetscErrorCode _myQepMonitor(QEP qep, int its, int nconv, PetscScalar *ei static PetscErrorCode _myQepMonitor(QEP qep, int its, int nconv,
gr, PetscScalar *eigr, PetscScalar *eigi,
PetscScalar *eigi, PetscReal* errest, int ne PetscReal *errest, int nest, void *mctx)
st,
void *mctx)
{ {
return _myMonitor("QEP", its, nconv, eigr, eigi, errest); return _myMonitor("QEP", its, nconv, eigr, eigi, errest);
} }
#else #else
static PetscErrorCode _myPepMonitor(PEP pep, int its, int nconv, PetscScalar *ei static PetscErrorCode _myPepMonitor(PEP pep, int its, int nconv,
gr, PetscScalar *eigr, PetscScalar *eigi,
PetscScalar *eigi, PetscReal* errest, int ne PetscReal *errest, int nest, void *mctx)
st,
void *mctx)
{ {
return _myMonitor("PEP", its, nconv, eigr, eigi, errest); return _myMonitor("PEP", its, nconv, eigr, eigi, errest);
} }
#endif #endif
static void _storeEigenVectors(struct DofData *DofData_P, int nconv, EPS eps, static void _storeEigenVectors(struct DofData *DofData_P, int nconv, EPS eps,
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5)
QEP qep, QEP qep,
#else #else
PEP pep, PEP pep, NEP nep,
NEP nep,
#endif #endif
int filterExpressionIndex) int filterExpressionIndex)
{ {
if (nconv <= 0) return; if(nconv <= 0) return;
// temporary (parallel) vectors to store real and imaginary part of eigenvecto // temporary (parallel) vectors to store real and imaginary part of
rs // eigenvectors
Vec xr, xi; Vec xr, xi;
if(!nep){ if(!nep) {
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6)
_try(MatGetVecs(DofData_P->M1.M, PETSC_NULL, &xr)); _try(MatGetVecs(DofData_P->M1.M, PETSC_NULL, &xr));
_try(MatGetVecs(DofData_P->M1.M, PETSC_NULL, &xi)); _try(MatGetVecs(DofData_P->M1.M, PETSC_NULL, &xi));
#else #else
_try(MatCreateVecs(DofData_P->M1.M, PETSC_NULL, &xr)); _try(MatCreateVecs(DofData_P->M1.M, PETSC_NULL, &xr));
_try(MatCreateVecs(DofData_P->M1.M, PETSC_NULL, &xi)); _try(MatCreateVecs(DofData_P->M1.M, PETSC_NULL, &xi));
#endif #endif
} }
else{ else {
_try(MatCreateVecs(DofData_P->M7.M, PETSC_NULL, &xr)); _try(MatCreateVecs(DofData_P->M7.M, PETSC_NULL, &xr));
_try(MatCreateVecs(DofData_P->M7.M, PETSC_NULL, &xi)); _try(MatCreateVecs(DofData_P->M7.M, PETSC_NULL, &xi));
} }
// temporary sequential vectors to transfer eigenvectors to getdp // temporary sequential vectors to transfer eigenvectors to getdp
Vec xr_seq, xi_seq; Vec xr_seq, xi_seq;
if(Message::GetCommSize() > 1){ if(Message::GetCommSize() > 1) {
PetscInt n; PetscInt n;
_try(VecGetSize(xr, &n)); _try(VecGetSize(xr, &n));
_try(VecCreateSeq(PETSC_COMM_SELF, n, &xr_seq)); _try(VecCreateSeq(PETSC_COMM_SELF, n, &xr_seq));
_try(VecCreateSeq(PETSC_COMM_SELF, n, &xi_seq)); _try(VecCreateSeq(PETSC_COMM_SELF, n, &xi_seq));
} }
Message::Info(" %-24s%-24s%-12s", "Re", "Im", "Relative error"); Message::Info(" %-24s%-24s%-12s", "Re", "Im", "Relative error");
bool newsol = false; bool newsol = false;
for (int i = 0; i < nconv; i++){ for(int i = 0; i < nconv; i++) {
PetscScalar kr, ki; PetscScalar kr, ki;
PetscReal error; PetscReal error;
if(eps){ if(eps) {
_try(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi)); _try(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi));
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6)
_try(EPSComputeRelativeError(eps, i, &error)); _try(EPSComputeRelativeError(eps, i, &error));
#else #else
_try(EPSComputeError(eps, i, EPS_ERROR_RELATIVE, &error)); _try(EPSComputeError(eps, i, EPS_ERROR_RELATIVE, &error));
#endif #endif
} }
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5)
else if (qep){ else if(qep) {
_try(QEPGetEigenpair(qep, i, &kr, &ki, xr, xi)); _try(QEPGetEigenpair(qep, i, &kr, &ki, xr, xi));
_try(QEPComputeRelativeError(qep, i, &error)); _try(QEPComputeRelativeError(qep, i, &error));
} }
#else #else
else if (pep){ else if(pep) {
_try(PEPGetEigenpair(pep, i, &kr, &ki, xr, xi)); _try(PEPGetEigenpair(pep, i, &kr, &ki, xr, xi));
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6)
_try(PEPComputeRelativeError(pep, i, &error)); _try(PEPComputeRelativeError(pep, i, &error));
#else #else
_try(PEPComputeError(pep, i, PEP_ERROR_BACKWARD, &error)); _try(PEPComputeError(pep, i, PEP_ERROR_BACKWARD, &error));
#endif #endif
} }
#endif #endif
else if(nep){ else if(nep) {
_try(NEPGetEigenpair(nep, i, &kr, &ki, xr, xi)); _try(NEPGetEigenpair(nep, i, &kr, &ki, xr, xi));
_try(NEPComputeError(nep, i, NEP_ERROR_RELATIVE, &error)); _try(NEPComputeError(nep, i, NEP_ERROR_RELATIVE, &error));
} }
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
PetscReal re = PetscRealPart(kr), im = PetscImaginaryPart(kr); PetscReal re = PetscRealPart(kr), im = PetscImaginaryPart(kr);
#else #else
PetscReal re = kr, im = ki; PetscReal re = kr, im = ki;
#endif #endif
double ore=0, oim=0; double ore = 0, oim = 0;
if(eps){ if(eps) {
Message::Info("EIG %03d w^2 = %s%.16e %s%.16e %3.6e", Message::Info("EIG %03d w^2 = %s%.16e %s%.16e %3.6e", i,
i, (re < 0) ? "" : " ", re, (im < 0) ? "" : " ", im, error); (re < 0) ? "" : " ", re, (im < 0) ? "" : " ", im, error);
double abs = sqrt(re * re + im * im), arg = atan2(im, re); double abs = sqrt(re * re + im * im), arg = atan2(im, re);
ore = sqrt(abs) * cos(0.5*arg); ore = sqrt(abs) * cos(0.5 * arg);
oim = sqrt(abs) * sin(0.5*arg); oim = sqrt(abs) * sin(0.5 * arg);
double fre = ore / 2. / M_PI, fim = oim / 2. / M_PI; double fre = ore / 2. / M_PI, fim = oim / 2. / M_PI;
Message::Info(" w = %s%.16e %s%.16e", Message::Info(" w = %s%.16e %s%.16e", (ore < 0) ? "" : " ", ore,
(ore < 0) ? "" : " ", ore, (oim < 0) ? "" : " ", oim); (oim < 0) ? "" : " ", oim);
Message::Info(" f = %s%.16e %s%.16e", Message::Info(" f = %s%.16e %s%.16e", (fre < 0) ? "" : " ", fre,
(fre < 0) ? "" : " ", fre, (fim < 0) ? "" : " ", fim); (fim < 0) ? "" : " ", fim);
} }
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5)
else if(qep){ else if(qep) {
// lambda = iw // lambda = iw
ore = im; ore = im;
oim = -re; oim = -re;
Message::Info("EIG %03d w = %s%.16e %s%.16e %3.6e", Message::Info("EIG %03d w = %s%.16e %s%.16e %3.6e", i,
i, (ore < 0) ? "" : " ", ore, (oim < 0) ? "" : " ", oim, err (ore < 0) ? "" : " ", ore, (oim < 0) ? "" : " ", oim,
or); error);
double fre = ore / 2. / M_PI, fim = oim / 2. / M_PI; double fre = ore / 2. / M_PI, fim = oim / 2. / M_PI;
Message::Info(" f = %s%.16e %s%.16e", Message::Info(" f = %s%.16e %s%.16e", (fre < 0) ? "" : " ", fre,
(fre < 0) ? "" : " ", fre, (fim < 0) ? "" : " ", fim); (fim < 0) ? "" : " ", fim);
} }
#else #else
else if (pep){ else if(pep) {
// lambda = iw // lambda = iw
ore = im; ore = im;
oim = -re; oim = -re;
Message::Info("EIG %03d w = %s%.16e %s%.16e %3.6e", Message::Info("EIG %03d w = %s%.16e %s%.16e %3.6e", i,
i, (ore < 0) ? "" : " ", ore, (oim < 0) ? "" : " ", oim, err (ore < 0) ? "" : " ", ore, (oim < 0) ? "" : " ", oim,
or); error);
double fre = ore / 2. / M_PI, fim = oim / 2. / M_PI; double fre = ore / 2. / M_PI, fim = oim / 2. / M_PI;
Message::Info(" f = %s%.16e %s%.16e", Message::Info(" f = %s%.16e %s%.16e", (fre < 0) ? "" : " ", fre,
(fre < 0) ? "" : " ", fre, (fim < 0) ? "" : " ", fim); (fim < 0) ? "" : " ", fim);
} }
#endif #endif
else if (nep){ else if(nep) {
// lambda = iw // lambda = iw
ore = im; ore = im;
oim = -re; oim = -re;
Message::Info("EIG %03d w = %s%.16e %s%.16e %3.6e", Message::Info("EIG %03d w = %s%.16e %s%.16e %3.6e", i,
i, (ore < 0) ? "" : " ", ore, (oim < 0) ? "" : " ", oim, err (ore < 0) ? "" : " ", ore, (oim < 0) ? "" : " ", oim,
or); error);
double fre = ore / 2. / M_PI, fim = oim / 2. / M_PI; double fre = ore / 2. / M_PI, fim = oim / 2. / M_PI;
Message::Info(" f = %s%.16e %s%.16e", Message::Info(" f = %s%.16e %s%.16e", (fre < 0) ? "" : " ", fre,
(fre < 0) ? "" : " ", fre, (fim < 0) ? "" : " ", fim); (fim < 0) ? "" : " ", fim);
} }
// update the current value of Time and TimeImag so that // update the current value of Time and TimeImag so that
// $EigenvalueReal and $EigenvalueImag are up-to-date // $EigenvalueReal and $EigenvalueImag are up-to-date
Current.Time = ore; Current.Time = ore;
Current.TimeImag = oim; Current.TimeImag = oim;
// test filter expression and continue without storing if false // test filter expression and continue without storing if false
if(filterExpressionIndex >= 0){ if(filterExpressionIndex >= 0) {
struct Value val; struct Value val;
Get_ValueOfExpressionByIndex(filterExpressionIndex, NULL, 0., 0., 0., &val Get_ValueOfExpressionByIndex(filterExpressionIndex, NULL, 0., 0., 0.,
); &val);
if(!val.Val[0]){ if(!val.Val[0]) {
Message::Debug("Skipping eigenvalue %g + i * %g", ore, oim); Message::Debug("Skipping eigenvalue %g + i * %g", ore, oim);
continue; continue;
} }
} }
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Re(Omega)" Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
, "/Re(Omega)",
std::vector<double>(1, ore)); std::vector<double>(1, ore));
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Im(Omega)" Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
, "/Im(Omega)",
std::vector<double>(1, oim)); std::vector<double>(1, oim));
// create new solution vector if necessary // create new solution vector if necessary
if(newsol) { if(newsol) {
struct Solution Solution_S; struct Solution Solution_S;
Solution_S.TimeFunctionValues = NULL; Solution_S.TimeFunctionValues = NULL;
LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof);
List_Add(DofData_P->Solutions, &Solution_S); List_Add(DofData_P->Solutions, &Solution_S);
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1); DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
} }
newsol = true; newsol = true;
DofData_P->CurrentSolution->Time = ore; DofData_P->CurrentSolution->Time = ore;
DofData_P->CurrentSolution->TimeImag = oim; DofData_P->CurrentSolution->TimeImag = oim;
DofData_P->CurrentSolution->TimeStep = (int)Current.TimeStep; DofData_P->CurrentSolution->TimeStep = (int)Current.TimeStep;
Free(DofData_P->CurrentSolution->TimeFunctionValues); Free(DofData_P->CurrentSolution->TimeFunctionValues);
DofData_P->CurrentSolution->TimeFunctionValues = NULL; DofData_P->CurrentSolution->TimeFunctionValues = NULL;
DofData_P->CurrentSolution->SolutionExist = 1; DofData_P->CurrentSolution->SolutionExist = 1;
// store eigenvector // store eigenvector
PetscScalar *tmpr, *tmpi; PetscScalar *tmpr, *tmpi;
if(Message::GetCommSize() == 1){ if(Message::GetCommSize() == 1) {
_try(VecGetArray(xr, &tmpr)); _try(VecGetArray(xr, &tmpr));
_try(VecGetArray(xi, &tmpi)); _try(VecGetArray(xi, &tmpi));
} }
else{ else {
_fillseq(xr, xr_seq); _fillseq(xr, xr_seq);
_fillseq(xi, xi_seq); _fillseq(xi, xi_seq);
_try(VecGetArray(xr_seq, &tmpr)); _try(VecGetArray(xr_seq, &tmpr));
_try(VecGetArray(xi_seq, &tmpi)); _try(VecGetArray(xi_seq, &tmpi));
} }
int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1; int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1;
for(int l = 0; l < DofData_P->NbrDof; l += incr){ for(int l = 0; l < DofData_P->NbrDof; l += incr) {
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
double var_r = (double)PetscRealPart(tmpr[l]); double var_r = (double)PetscRealPart(tmpr[l]);
double var_i = (double)PetscImaginaryPart(tmpr[l]); double var_i = (double)PetscImaginaryPart(tmpr[l]);
#else #else
double var_r = (double)tmpr[l]; double var_r = (double)tmpr[l];
double var_i = (double)tmpi[l]; double var_i = (double)tmpi[l];
#endif #endif
if(Current.NbrHar == 2) if(Current.NbrHar == 2)
LinAlg_SetComplexInVector(var_r, var_i, &DofData_P->CurrentSolution->x, LinAlg_SetComplexInVector(var_r, var_i, &DofData_P->CurrentSolution->x,
l, l+1); l, l + 1);
else else
LinAlg_SetDoubleInVector(var_r, &DofData_P->CurrentSolution->x, l); LinAlg_SetDoubleInVector(var_r, &DofData_P->CurrentSolution->x, l);
} }
if(Message::GetCommSize() == 1){ if(Message::GetCommSize() == 1) {
_try(VecRestoreArray(xr, &tmpr)); _try(VecRestoreArray(xr, &tmpr));
_try(VecRestoreArray(xi, &tmpi)); _try(VecRestoreArray(xi, &tmpi));
} }
else{ else {
_try(VecRestoreArray(xr_seq, &tmpr)); _try(VecRestoreArray(xr_seq, &tmpr));
_try(VecRestoreArray(xi_seq, &tmpi)); _try(VecRestoreArray(xi_seq, &tmpi));
} }
LinAlg_AssembleVector(&DofData_P->CurrentSolution->x); LinAlg_AssembleVector(&DofData_P->CurrentSolution->x);
// SLEPc returns eigenvectors normalized in L-2 norm. Renormalize them in // SLEPc returns eigenvectors normalized in L-2 norm. Renormalize them in
// L-infty norm so that the absolute value of the largest element is 1 // L-infty norm so that the absolute value of the largest element is 1
double norm = 0; double norm = 0;
LinAlg_VectorNormInf(&DofData_P->CurrentSolution->x, &norm); LinAlg_VectorNormInf(&DofData_P->CurrentSolution->x, &norm);
if(norm > 1.e-16) if(norm > 1.e-16)
LinAlg_ProdVectorDouble(&DofData_P->CurrentSolution->x, 1. / norm, LinAlg_ProdVectorDouble(&DofData_P->CurrentSolution->x, 1. / norm,
&DofData_P->CurrentSolution->x); &DofData_P->CurrentSolution->x);
// increment the global timestep counter so that a future // increment the global timestep counter so that a future
// GenerateSystem knows which solutions exist // GenerateSystem knows which solutions exist
Current.TimeStep += 1.; Current.TimeStep += 1.;
} }
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)
_try(VecDestroy(&xr)); _try(VecDestroy(&xr));
_try(VecDestroy(&xi)); _try(VecDestroy(&xi));
if(Message::GetCommSize() > 1){ if(Message::GetCommSize() > 1) {
_try(VecDestroy(&xr_seq)); _try(VecDestroy(&xr_seq));
_try(VecDestroy(&xi_seq)); _try(VecDestroy(&xi_seq));
} }
#else #else
_try(VecDestroy(xr)); _try(VecDestroy(xr));
_try(VecDestroy(xi)); _try(VecDestroy(xi));
if(Message::GetCommSize() > 1){ if(Message::GetCommSize() > 1) {
_try(VecDestroy(xr_seq)); _try(VecDestroy(xr_seq));
_try(VecDestroy(xi_seq)); _try(VecDestroy(xi_seq));
} }
#endif #endif
} }
static void _storeExpansion(struct DofData *DofData_P, int nbApplyResolventRealF static void _storeExpansion(struct DofData *DofData_P,
reqs, int nbApplyResolventRealFreqs,
std::vector<PetscReal> tabApplyResolventRealFreqs, std::vector<PetscReal> tabApplyResolventRealFreqs,
std::vector<Vec> &ListOfExpansionResults) std::vector<Vec> &ListOfExpansionResults)
{ {
Vec x; Vec x;
_try(MatCreateVecs(DofData_P->M7.M, PETSC_NULL, &x)); _try(MatCreateVecs(DofData_P->M7.M, PETSC_NULL, &x));
// temporary sequential vectors to transfer eigenvectors to getdp // temporary sequential vectors to transfer eigenvectors to getdp
Vec x_seq; Vec x_seq;
if(Message::GetCommSize() > 1){ if(Message::GetCommSize() > 1) {
PetscInt n; PetscInt n;
_try(VecGetSize(x, &n)); _try(VecGetSize(x, &n));
_try(VecCreateSeq(PETSC_COMM_SELF, n, &x_seq)); _try(VecCreateSeq(PETSC_COMM_SELF, n, &x_seq));
} }
bool newsol = true; bool newsol = true;
// char fname[100]; // char fname[100];
// static PetscViewer myviewer; // static PetscViewer myviewer;
for (int i = 0; i < nbApplyResolventRealFreqs; i++){ for(int i = 0; i < nbApplyResolventRealFreqs; i++) {
x = ListOfExpansionResults[i]; x = ListOfExpansionResults[i];
// // store vectors in petsc/matlab format // // store vectors in petsc/matlab format
// sprintf(fname,"applyresolvent_vec_result_%03d.m",i); // sprintf(fname,"applyresolvent_vec_result_%03d.m",i);
// PetscViewerASCIIOpen(PETSC_COMM_WORLD, fname , &myviewer); // PetscViewerASCIIOpen(PETSC_COMM_WORLD, fname , &myviewer);
// PetscViewerPushFormat(myviewer, PETSC_VIEWER_ASCII_MATLAB); // PetscViewerPushFormat(myviewer, PETSC_VIEWER_ASCII_MATLAB);
// VecView(x,myviewer); // VecView(x,myviewer);
// PetscViewerPopFormat(myviewer); // PetscViewerPopFormat(myviewer);
Current.Time = tabApplyResolventRealFreqs[i]; Current.Time = tabApplyResolventRealFreqs[i];
Current.TimeImag = 0; Current.TimeImag = 0;
// create new solution vector if necessary // create new solution vector if necessary
if(newsol) { if(newsol) {
struct Solution Solution_S; struct Solution Solution_S;
Solution_S.TimeFunctionValues = NULL; Solution_S.TimeFunctionValues = NULL;
LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof);
List_Add(DofData_P->Solutions, &Solution_S); List_Add(DofData_P->Solutions, &Solution_S);
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1); DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
} }
DofData_P->CurrentSolution->Time = tabApplyResolventRealFreqs[i]; DofData_P->CurrentSolution->Time = tabApplyResolventRealFreqs[i];
DofData_P->CurrentSolution->TimeImag = 0; DofData_P->CurrentSolution->TimeImag = 0;
DofData_P->CurrentSolution->TimeStep = (int)Current.TimeStep; DofData_P->CurrentSolution->TimeStep = (int)Current.TimeStep;
Free(DofData_P->CurrentSolution->TimeFunctionValues); Free(DofData_P->CurrentSolution->TimeFunctionValues);
DofData_P->CurrentSolution->TimeFunctionValues = NULL; DofData_P->CurrentSolution->TimeFunctionValues = NULL;
DofData_P->CurrentSolution->SolutionExist = 1; DofData_P->CurrentSolution->SolutionExist = 1;
// store eigenvector // store eigenvector
PetscScalar *tmp; PetscScalar *tmp;
if(Message::GetCommSize() == 1){ if(Message::GetCommSize() == 1) { _try(VecGetArray(x, &tmp)); }
_try(VecGetArray(x, &tmp)); else {
}
else{
_fillseq(x, x_seq); _fillseq(x, x_seq);
_try(VecGetArray(x_seq, &tmp)); _try(VecGetArray(x_seq, &tmp));
} }
int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1; int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1;
for(int l = 0; l < DofData_P->NbrDof; l += incr){ for(int l = 0; l < DofData_P->NbrDof; l += incr) {
double var_r = (double)PetscRealPart(tmp[l]); double var_r = (double)PetscRealPart(tmp[l]);
double var_i = (double)PetscImaginaryPart(tmp[l]); double var_i = (double)PetscImaginaryPart(tmp[l]);
if(Current.NbrHar == 2) if(Current.NbrHar == 2)
LinAlg_SetComplexInVector(var_r, var_i, &DofData_P->CurrentSolution->x, LinAlg_SetComplexInVector(var_r, var_i, &DofData_P->CurrentSolution->x,
l, l+1); l, l + 1);
else else
LinAlg_SetDoubleInVector(var_r, &DofData_P->CurrentSolution->x, l); LinAlg_SetDoubleInVector(var_r, &DofData_P->CurrentSolution->x, l);
} }
if(Message::GetCommSize() == 1){ if(Message::GetCommSize() == 1) { _try(VecRestoreArray(x, &tmp)); }
_try(VecRestoreArray(x, &tmp)); else {
}
else{
_try(VecRestoreArray(x_seq, &tmp)); _try(VecRestoreArray(x_seq, &tmp));
} }
LinAlg_AssembleVector(&DofData_P->CurrentSolution->x); LinAlg_AssembleVector(&DofData_P->CurrentSolution->x);
// increment the global timestep counter so that a future // increment the global timestep counter so that a future
// GenerateSystem knows which solutions exist // GenerateSystem knows which solutions exist
Current.TimeStep += 1.; Current.TimeStep += 1.;
} }
// _try(PetscViewerDestroy(&myviewer)); // _try(PetscViewerDestroy(&myviewer));
_try(VecDestroy(&x)); _try(VecDestroy(&x));
if(Message::GetCommSize() > 1){ if(Message::GetCommSize() > 1) { _try(VecDestroy(&x_seq)); }
_try(VecDestroy(&x_seq));
}
} }
static void _linearEVP(struct DofData * DofData_P, int numEigenValues, static void _linearEVP(struct DofData *DofData_P, int numEigenValues,
double shift_r, double shift_i, int filterExpressionIndex double shift_r, double shift_i,
) int filterExpressionIndex)
{ {
Message::Info("Solving linear eigenvalue problem"); Message::Info("Solving linear eigenvalue problem");
// GetDP notation: -w^2 M3 x (+ iw M2 x) + M1 x = 0 // GetDP notation: -w^2 M3 x (+ iw M2 x) + M1 x = 0
// SLEPC notation for generalized linear EVP: A x - \lambda B x = 0 // SLEPC notation for generalized linear EVP: A x - \lambda B x = 0
Mat A = DofData_P->M1.M; Mat A = DofData_P->M1.M;
Mat B = DofData_P->M3.M; Mat B = DofData_P->M3.M;
EPS eps; EPS eps;
_try(EPSCreate(PETSC_COMM_WORLD, &eps)); _try(EPSCreate(PETSC_COMM_WORLD, &eps));
_try(EPSSetOperators(eps, A, B)); _try(EPSSetOperators(eps, A, B));
skipping to change at line 460 skipping to change at line 470
// force options specified directly as arguments // force options specified directly as arguments
if(numEigenValues) if(numEigenValues)
_try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); _try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
// apply shift-and-invert transformation // apply shift-and-invert transformation
ST st; ST st;
_try(EPSGetST(eps, &st)); _try(EPSGetST(eps, &st));
_try(STSetType(st, STSINVERT)); _try(STSetType(st, STSINVERT));
_try(EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE)); _try(EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE));
if(shift_r || shift_i){ if(shift_r || shift_i) {
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
PetscScalar shift = shift_r + PETSC_i * shift_i; PetscScalar shift = shift_r + PETSC_i * shift_i;
#else #else
PetscScalar shift = shift_r; PetscScalar shift = shift_r;
if(shift_i) if(shift_i)
Message::Warning("Imaginary part of shift discarded: use PETSc with comple Message::Warning(
x numbers"); "Imaginary part of shift discarded: use PETSc with complex numbers");
#endif #endif
//_try(STSetShift(st, shift)); //_try(STSetShift(st, shift));
_try(EPSSetTarget(eps, shift)); _try(EPSSetTarget(eps, shift));
} }
KSP ksp; KSP ksp;
_try(STGetKSP(st, &ksp)); _try(STGetKSP(st, &ksp));
_try(KSPSetType(ksp, "preonly")); _try(KSPSetType(ksp, "preonly"));
PC pc; PC pc;
_try(KSPGetPC(ksp, &pc)); _try(KSPGetPC(ksp, &pc));
_try(PCSetType(pc, PCLU)); _try(PCSetType(pc, PCLU));
#if (PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS) #if(PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS)
_try(PCFactorSetMatSolverPackage(pc, "mumps")); _try(PCFactorSetMatSolverPackage(pc, "mumps"));
#endif #endif
// print info // print info
#if (PETSC_VERSION_RELASE == 0) || ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION #if(PETSC_VERSION_RELASE == 0) || \
_MINOR == 4)) ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR == 4))
const char *type = ""; const char *type = "";
#else #else
const EPSType type; const EPSType type;
#endif #endif
_try(EPSGetType(eps, &type)); _try(EPSGetType(eps, &type));
Message::Info("SLEPc solution method: %s", type); Message::Info("SLEPc solution method: %s", type);
PetscInt nev; PetscInt nev;
_try(EPSGetDimensions(eps, &nev, PETSC_NULL, PETSC_NULL)); _try(EPSGetDimensions(eps, &nev, PETSC_NULL, PETSC_NULL));
Message::Info("SLEPc number of requested eigenvalues: %d", nev); Message::Info("SLEPc number of requested eigenvalues: %d", nev);
PetscReal tol; PetscReal tol;
skipping to change at line 527 skipping to change at line 539
// get number of converged approximate eigenpairs // get number of converged approximate eigenpairs
PetscInt nconv; PetscInt nconv;
_try(EPSGetConverged(eps, &nconv)); _try(EPSGetConverged(eps, &nconv));
Message::Info("SLEPc number of converged eigenpairs: %d", nconv); Message::Info("SLEPc number of converged eigenpairs: %d", nconv);
// ignore additional eigenvalues if we get more than what we asked // ignore additional eigenvalues if we get more than what we asked
if(nconv > nev) nconv = nev; if(nconv > nev) nconv = nev;
// print eigenvalues and store eigenvectors in DofData // print eigenvalues and store eigenvectors in DofData
_storeEigenVectors(DofData_P, nconv, eps, PETSC_NULL,PETSC_NULL , filterExpres _storeEigenVectors(DofData_P, nconv, eps, PETSC_NULL, PETSC_NULL,
sionIndex); filterExpressionIndex);
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)
_try(EPSDestroy(&eps)); _try(EPSDestroy(&eps));
#else #else
_try(EPSDestroy(eps)); _try(EPSDestroy(eps));
#endif #endif
} }
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5)
// SLEPc < 3.5 interface using QEP // SLEPc < 3.5 interface using QEP
static void _quadraticEVP(struct DofData * DofData_P, int numEigenValues, static void _quadraticEVP(struct DofData *DofData_P, int numEigenValues,
double shift_r, double shift_i, int filterExpressionIn double shift_r, double shift_i,
dex) int filterExpressionIndex)
{ {
Message::Info("Solving quadratic eigenvalue problem using QEP"); Message::Info("Solving quadratic eigenvalue problem using QEP");
// GetDP notation: -w^2 M3 x + iw M2 x + M1 x = 0 // GetDP notation: -w^2 M3 x + iw M2 x + M1 x = 0
// SLEPC notations for quadratic EVP: (\lambda^2 M + \lambda C + K) x = 0 // SLEPC notations for quadratic EVP: (\lambda^2 M + \lambda C + K) x = 0
Mat M = DofData_P->M3.M; Mat M = DofData_P->M3.M;
Mat C = DofData_P->M2.M; Mat C = DofData_P->M2.M;
Mat K = DofData_P->M1.M; Mat K = DofData_P->M1.M;
QEP qep; QEP qep;
skipping to change at line 563 skipping to change at line 577
_try(QEPSetProblemType(qep, QEP_GENERAL)); _try(QEPSetProblemType(qep, QEP_GENERAL));
// set some default options // set some default options
_try(QEPSetDimensions(qep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); _try(QEPSetDimensions(qep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
_try(QEPSetTolerances(qep, 1.e-6, 100)); _try(QEPSetTolerances(qep, 1.e-6, 100));
_try(QEPSetType(qep, QEPLINEAR)); _try(QEPSetType(qep, QEPLINEAR));
_try(QEPSetWhichEigenpairs(qep, QEP_SMALLEST_MAGNITUDE)); _try(QEPSetWhichEigenpairs(qep, QEP_SMALLEST_MAGNITUDE));
_try(QEPMonitorSet(qep, _myQepMonitor, PETSC_NULL, PETSC_NULL)); _try(QEPMonitorSet(qep, _myQepMonitor, PETSC_NULL, PETSC_NULL));
// if we linearize we can set additional options // if we linearize we can set additional options
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR == 4) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR == 4)
const char *type = ""; const char *type = "";
#else #else
const QEPType type; const QEPType type;
#endif #endif
_try(QEPGetType(qep, &type)); _try(QEPGetType(qep, &type));
if(!strcmp(type, QEPLINEAR)){ if(!strcmp(type, QEPLINEAR)) {
EPS eps; EPS eps;
_try(QEPLinearGetEPS(qep, &eps)); _try(QEPLinearGetEPS(qep, &eps));
_try(EPSSetType(eps, EPSKRYLOVSCHUR)); _try(EPSSetType(eps, EPSKRYLOVSCHUR));
// apply shift-and-invert transformation // apply shift-and-invert transformation
ST st; ST st;
_try(EPSGetST(eps, &st)); _try(EPSGetST(eps, &st));
_try(STSetType(st, STSINVERT)); _try(STSetType(st, STSINVERT));
_try(EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE)); _try(EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE));
if(shift_r || shift_i){ if(shift_r || shift_i) {
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
PetscScalar shift = shift_r + PETSC_i * shift_i; PetscScalar shift = shift_r + PETSC_i * shift_i;
#else #else
PetscScalar shift = shift_r; PetscScalar shift = shift_r;
if(shift_i) if(shift_i)
Message::Warning("Imaginary part of shift discarded: use PETSc with comp Message::Warning(
lex numbers"); "Imaginary part of shift discarded: use PETSc with complex numbers");
#endif #endif
_try(EPSSetTarget(eps, shift)); _try(EPSSetTarget(eps, shift));
} }
_try(QEPLinearSetExplicitMatrix(qep, PETSC_TRUE)); _try(QEPLinearSetExplicitMatrix(qep, PETSC_TRUE));
Message::Info("SLEPc forcing explicit construction of matrix"); Message::Info("SLEPc forcing explicit construction of matrix");
KSP ksp; KSP ksp;
_try(STGetKSP(st, &ksp)); _try(STGetKSP(st, &ksp));
_try(KSPSetType(ksp, "preonly")); _try(KSPSetType(ksp, "preonly"));
PC pc; PC pc;
_try(KSPGetPC(ksp, &pc)); _try(KSPGetPC(ksp, &pc));
_try(PCSetType(pc, PCLU)); _try(PCSetType(pc, PCLU));
#if (PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS) #if(PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS)
_try(PCFactorSetMatSolverPackage(pc, "mumps")); _try(PCFactorSetMatSolverPackage(pc, "mumps"));
#endif #endif
} }
// override these options at runtime, if necessary // override these options at runtime, if necessary
_try(QEPSetFromOptions(qep)); _try(QEPSetFromOptions(qep));
// force options specified directly as arguments // force options specified directly as arguments
if(numEigenValues) if(numEigenValues)
_try(QEPSetDimensions(qep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); _try(QEPSetDimensions(qep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
skipping to change at line 645 skipping to change at line 660
// get number of converged approximate eigenpairs // get number of converged approximate eigenpairs
PetscInt nconv; PetscInt nconv;
_try(QEPGetConverged(qep, &nconv)); _try(QEPGetConverged(qep, &nconv));
Message::Info("SLEPc number of converged eigenpairs: %d", nconv); Message::Info("SLEPc number of converged eigenpairs: %d", nconv);
// ignore additional eigenvalues if we get more than what we asked // ignore additional eigenvalues if we get more than what we asked
if(nconv > nev) nconv = nev; if(nconv > nev) nconv = nev;
// print eigenvalues and store eigenvectors in DofData // print eigenvalues and store eigenvectors in DofData
_storeEigenVectors(DofData_P, nconv, PETSC_NULL, qep, PETSC_NULL, filterExpres _storeEigenVectors(DofData_P, nconv, PETSC_NULL, qep, PETSC_NULL,
sionIndex); filterExpressionIndex);
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)
_try(QEPDestroy(&qep)); _try(QEPDestroy(&qep));
#else #else
_try(QEPDestroy(qep)); _try(QEPDestroy(qep));
#endif #endif
} }
#else // SLEPc >= 3.5 interface using PEP #else // SLEPc >= 3.5 interface using PEP
static void _quadraticEVP(struct DofData * DofData_P, int numEigenValues, static void _quadraticEVP(struct DofData *DofData_P, int numEigenValues,
double shift_r, double shift_i, int filterExpressionIn double shift_r, double shift_i,
dex) int filterExpressionIndex)
{ {
Message::Info("Solving quadratic eigenvalue problem using PEP"); Message::Info("Solving quadratic eigenvalue problem using PEP");
// GetDP notation: -w^2 M3 x + iw M2 x + M1 x = 0 // GetDP notation: -w^2 M3 x + iw M2 x + M1 x = 0
// SLEPC notations for quadratic EVP: (\lambda^2 A[2] + \lambda A[1] + A[0]) x // SLEPC notations for quadratic EVP: (\lambda^2 A[2] + \lambda A[1] + A[0]) x
= 0 // = 0
PEP pep; PEP pep;
ST st; ST st;
KSP ksp; KSP ksp;
PC pc; PC pc;
PEPType type; PEPType type;
Mat A[3] = {DofData_P->M1.M, DofData_P->M2.M, DofData_P->M3.M}; Mat A[3] = {DofData_P->M1.M, DofData_P->M2.M, DofData_P->M3.M};
_try(PEPCreate(PETSC_COMM_WORLD, &pep)); _try(PEPCreate(PETSC_COMM_WORLD, &pep));
_try(PEPSetOperators(pep, 3, A)); _try(PEPSetOperators(pep, 3, A));
_try(PEPSetProblemType(pep, PEP_GENERAL)); _try(PEPSetProblemType(pep, PEP_GENERAL));
skipping to change at line 688 skipping to change at line 706
_try(PEPSetTolerances(pep, 1.e-6, 100)); _try(PEPSetTolerances(pep, 1.e-6, 100));
// _try(PEPSetType(pep, PEPLINEAR)); // _try(PEPSetType(pep, PEPLINEAR));
_try(PEPSetType(pep, PEPTOAR)); _try(PEPSetType(pep, PEPTOAR));
_try(PEPGetST(pep, &st)); _try(PEPGetST(pep, &st));
_try(STSetType(st, STSINVERT)); _try(STSetType(st, STSINVERT));
_try(PEPSetWhichEigenpairs(pep, PEP_TARGET_MAGNITUDE)); _try(PEPSetWhichEigenpairs(pep, PEP_TARGET_MAGNITUDE));
_try(STGetKSP(st, &ksp)); _try(STGetKSP(st, &ksp));
_try(KSPSetType(ksp, KSPPREONLY)); _try(KSPSetType(ksp, KSPPREONLY));
_try(KSPGetPC(ksp, &pc)); _try(KSPGetPC(ksp, &pc));
_try(PCSetType(pc, PCLU)); _try(PCSetType(pc, PCLU));
#if (PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS) #if(PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS)
_try(PCFactorSetMatSolverPackage(pc, "mumps")); _try(PCFactorSetMatSolverPackage(pc, "mumps"));
#endif #endif
_try(PEPMonitorSet(pep, _myPepMonitor, PETSC_NULL, PETSC_NULL)); _try(PEPMonitorSet(pep, _myPepMonitor, PETSC_NULL, PETSC_NULL));
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
PetscScalar shift = shift_r + PETSC_i * shift_i; PetscScalar shift = shift_r + PETSC_i * shift_i;
#else #else
PetscScalar shift = shift_r; PetscScalar shift = shift_r;
Message::Warning("Imaginary part of shift discarded: use PETSc with complex nu Message::Warning(
mbers"); "Imaginary part of shift discarded: use PETSc with complex numbers");
#endif #endif
_try(PEPSetTarget(pep, shift)); _try(PEPSetTarget(pep, shift));
// _try(PEPMonitorSet(pep, _myPepMonitor, PETSC_NULL, PETSC_NULL)); // _try(PEPMonitorSet(pep, _myPepMonitor, PETSC_NULL, PETSC_NULL));
// //
// // shift // // shift
// PetscScalar shift = 0; // PetscScalar shift = 0;
// if(shift_r || shift_i){ // if(shift_r || shift_i){
// #if defined(PETSC_USE_COMPLEX) // #if defined(PETSC_USE_COMPLEX)
// shift = shift_r + PETSC_i * shift_i; // shift = shift_r + PETSC_i * shift_i;
// #else // #else
// shift = shift_r; // shift = shift_r;
// if(shift_i) // if(shift_i)
// Message::Warning("Imaginary part of shift discarded: use PETSc with com // Message::Warning("Imaginary part of shift discarded: use PETSc with
plex numbers"); // complex numbers");
// #endif // #endif
// } // }
// //
// if(shift_r || shift_i){ // if(shift_r || shift_i){
// _try(PEPSetTarget(pep, shift)); // _try(PEPSetTarget(pep, shift));
// _try(PEPSetWhichEigenpairs(pep, PEP_TARGET_MAGNITUDE)); // _try(PEPSetWhichEigenpairs(pep, PEP_TARGET_MAGNITUDE));
// } // }
// //
// // if we linearize we can set additional options // // if we linearize we can set additional options
// const char *type = ""; // const char *type = "";
// _try(PEPGetType(pep, &type)); // _try(PEPGetType(pep, &type));
// if(!strcmp(type, PEPLINEAR)){ // if(!strcmp(type, PEPLINEAR)){
// EPS eps; // EPS eps;
// _try(PEPLinearGetEPS(pep, &eps)); // _try(PEPLinearGetEPS(pep, &eps));
// _try(EPSSetType(eps, EPSKRYLOVSCHUR)); // _try(EPSSetType(eps, EPSKRYLOVSCHUR));
// ST st; // ST st;
// _try(EPSGetST(eps, &st)); // _try(EPSGetST(eps, &st));
// _try(STSetType(st, STSINVERT)); // _try(STSetType(st, STSINVERT));
// _try(EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE)); // _try(EPSSetWhichEigenpairs(eps, EPS_TARGET_MAGNITUDE));
// if(shift_r || shift_i){ // if(shift_r || shift_i){
// _try(EPSSetTarget(eps, shift)); // _try(EPSSetTarget(eps, shift));
// } // }
// _try(PEPLinearSetExplicitMatrix(pep, PETSC_TRUE)); // _try(PEPLinearSetExplicitMatrix(pep, PETSC_TRUE));
// Message::Info("SLEPc forcing explicit construction of matrix"); // Message::Info("SLEPc forcing explicit construction of matrix");
// KSP ksp; // KSP ksp;
// _try(STGetKSP(st, &ksp)); // _try(STGetKSP(st, &ksp));
// _try(KSPSetType(ksp, "preonly")); // _try(KSPSetType(ksp, "preonly"));
// PC pc; // PC pc;
// _try(KSPGetPC(ksp, &pc)); // _try(KSPGetPC(ksp, &pc));
// _try(PCSetType(pc, PCLU)); // _try(PCSetType(pc, PCLU));
// #if (PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS) // #if (PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS)
// _try(PCFactorSetMatSolverPackage(pc, "mumps")); // _try(PCFactorSetMatSolverPackage(pc, "mumps"));
// #endif // #endif
// _try(EPSSetFromOptions(eps)); // _try(EPSSetFromOptions(eps));
// } // }
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6)
_try(PEPSetScale(pep, PEP_SCALE_SCALAR, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECI _try(PEPSetScale(pep, PEP_SCALE_SCALAR, PETSC_DECIDE, PETSC_DECIDE,
DE)); PETSC_DECIDE));
#else #else
_try(PEPSetScale(pep, PEP_SCALE_SCALAR, PETSC_DECIDE, PETSC_NULL, PETSC_NULL, _try(PEPSetScale(pep, PEP_SCALE_SCALAR, PETSC_DECIDE, PETSC_NULL, PETSC_NULL,
PETSC_DECIDE, PETSC_DECIDE)); PETSC_DECIDE, PETSC_DECIDE));
#endif #endif
// override these options at runtime, petsc-style // override these options at runtime, petsc-style
_try(PEPSetFromOptions(pep)); _try(PEPSetFromOptions(pep));
// force options specified directly as arguments // force options specified directly as arguments
if(numEigenValues){ if(numEigenValues) {
_try(PEPSetDimensions(pep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); _try(PEPSetDimensions(pep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
} }
// print info // print info
_try(PEPGetType(pep, &type)); _try(PEPGetType(pep, &type));
Message::Info("SLEPc solution method: %s", type); Message::Info("SLEPc solution method: %s", type);
PetscInt nev; PetscInt nev;
_try(PEPGetDimensions(pep, &nev, PETSC_NULL, PETSC_NULL)); _try(PEPGetDimensions(pep, &nev, PETSC_NULL, PETSC_NULL));
Message::Info("SLEPc number of requested eigenvalues: %d", nev); Message::Info("SLEPc number of requested eigenvalues: %d", nev);
PetscReal tol; PetscReal tol;
skipping to change at line 797 skipping to change at line 818
// get number of converged approximate eigenpairs // get number of converged approximate eigenpairs
PetscInt nconv; PetscInt nconv;
_try(PEPGetConverged(pep, &nconv)); _try(PEPGetConverged(pep, &nconv));
Message::Info("SLEPc number of converged eigenpairs: %d", nconv); Message::Info("SLEPc number of converged eigenpairs: %d", nconv);
// ignore additional eigenvalues if we get more than what we asked // ignore additional eigenvalues if we get more than what we asked
if(nconv > nev) nconv = nev; if(nconv > nev) nconv = nev;
// print eigenvalues and store eigenvectors in DofData // print eigenvalues and store eigenvectors in DofData
_storeEigenVectors(DofData_P, nconv, PETSC_NULL, pep, PETSC_NULL, filterExpres _storeEigenVectors(DofData_P, nconv, PETSC_NULL, pep, PETSC_NULL,
sionIndex); filterExpressionIndex);
_try(PEPDestroy(&pep)); _try(PEPDestroy(&pep));
} }
static void _polynomialEVP(struct DofData * DofData_P, int numEigenValues, static void _polynomialEVP(struct DofData *DofData_P, int numEigenValues,
double shift_r, double shift_i, int filterExpressionI double shift_r, double shift_i,
ndex) int filterExpressionIndex)
{ {
PEP pep; PEP pep;
ST st; ST st;
KSP ksp; KSP ksp;
PC pc; PC pc;
Message::Info("Solving polynomial eigenvalue problem using PEP"); Message::Info("Solving polynomial eigenvalue problem using PEP");
_try(PEPCreate(PETSC_COMM_WORLD, &pep)); _try(PEPCreate(PETSC_COMM_WORLD, &pep));
if(DofData_P->Flag_Init[6]){ if(DofData_P->Flag_Init[6]) {
Message::Info("Solving polynomial i*w^5 M6 x + w^4 M5 x + -iw^3 M4 x +" Message::Info(
" -w^2 M3 x + iw M2 x + M1 x = 0 eigenvalue problem using PEP" "Solving polynomial i*w^5 M6 x + w^4 M5 x + -iw^3 M4 x +"
); " -w^2 M3 x + iw M2 x + M1 x = 0 eigenvalue problem using PEP");
Mat A[6] = {DofData_P->M1.M, DofData_P->M2.M, DofData_P->M3.M, Mat A[6] = {DofData_P->M1.M, DofData_P->M2.M, DofData_P->M3.M,
DofData_P->M4.M, DofData_P->M5.M, DofData_P->M6.M}; DofData_P->M4.M, DofData_P->M5.M, DofData_P->M6.M};
_try(PEPSetOperators(pep, 6, A)); _try(PEPSetOperators(pep, 6, A));
} }
if(DofData_P->Flag_Init[5] && !DofData_P->Flag_Init[6]){ if(DofData_P->Flag_Init[5] && !DofData_P->Flag_Init[6]) {
Message::Info("Solving polynomial w^4 M5 x + -iw^3 M4 x + -w^2 M3 x + " Message::Info("Solving polynomial w^4 M5 x + -iw^3 M4 x + -w^2 M3 x + "
"iw M2 x + M1 x = 0 eigenvalue problem using PEP"); "iw M2 x + M1 x = 0 eigenvalue problem using PEP");
Mat A[5] = {DofData_P->M1.M, DofData_P->M2.M, DofData_P->M3.M, Mat A[5] = {DofData_P->M1.M, DofData_P->M2.M, DofData_P->M3.M,
DofData_P->M4.M, DofData_P->M5.M}; DofData_P->M4.M, DofData_P->M5.M};
_try(PEPSetOperators(pep, 5, A)); _try(PEPSetOperators(pep, 5, A));
} }
if(!DofData_P->Flag_Init[5] && !DofData_P->Flag_Init[6]){ if(!DofData_P->Flag_Init[5] && !DofData_P->Flag_Init[6]) {
Message::Info("Solving polynomial -iw^3 M4 x + -w^2 M3 x + iw M2 x + " Message::Info("Solving polynomial -iw^3 M4 x + -w^2 M3 x + iw M2 x + "
"M1 x = 0 eigenvalue problem using PEP"); "M1 x = 0 eigenvalue problem using PEP");
Mat A[4] = {DofData_P->M1.M, DofData_P->M2.M, DofData_P->M3.M, Mat A[4] = {DofData_P->M1.M, DofData_P->M2.M, DofData_P->M3.M,
DofData_P->M4.M}; DofData_P->M4.M};
_try(PEPSetOperators(pep, 4, A)); _try(PEPSetOperators(pep, 4, A));
} }
_try(PEPSetProblemType(pep, PEP_GENERAL)); _try(PEPSetProblemType(pep, PEP_GENERAL));
// set some default options // set some default options
_try(PEPSetDimensions(pep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); _try(PEPSetDimensions(pep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
_try(PEPSetTolerances(pep, 1.e-6, 100)); _try(PEPSetTolerances(pep, 1.e-6, 100));
_try(PEPSetType(pep, PEPTOAR)); _try(PEPSetType(pep, PEPTOAR));
_try(PEPGetST(pep, &st)); _try(PEPGetST(pep, &st));
_try(STSetType(st, STSINVERT)); _try(STSetType(st, STSINVERT));
_try(PEPSetWhichEigenpairs(pep, PEP_TARGET_MAGNITUDE)); _try(PEPSetWhichEigenpairs(pep, PEP_TARGET_MAGNITUDE));
_try(STGetKSP(st, &ksp)); _try(STGetKSP(st, &ksp));
_try(KSPSetType(ksp, KSPPREONLY)); _try(KSPSetType(ksp, KSPPREONLY));
_try(KSPGetPC(ksp, &pc)); _try(KSPGetPC(ksp, &pc));
_try(PCSetType(pc, PCLU)); _try(PCSetType(pc, PCLU));
#if (PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS) #if(PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS)
_try(PCFactorSetMatSolverPackage(pc, "mumps")); _try(PCFactorSetMatSolverPackage(pc, "mumps"));
#endif #endif
_try(PEPMonitorSet(pep, _myPepMonitor, PETSC_NULL, PETSC_NULL)); _try(PEPMonitorSet(pep, _myPepMonitor, PETSC_NULL, PETSC_NULL));
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
PetscScalar shift = shift_r + PETSC_i * shift_i; PetscScalar shift = shift_r + PETSC_i * shift_i;
#else #else
PetscScalar shift = shift_r; PetscScalar shift = shift_r;
#endif #endif
_try(PEPSetTarget(pep, shift)); _try(PEPSetTarget(pep, shift));
// _try(PEPSetScale(pep,PEP_SCALE_BOTH,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE) ); // _try(PEPSetScale(pep,PEP_SCALE_BOTH,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE) );
_try(PEPSetFromOptions(pep)); _try(PEPSetFromOptions(pep));
// force options specified directly as arguments // force options specified directly as arguments
if(numEigenValues){ if(numEigenValues) {
_try(PEPSetDimensions(pep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); _try(PEPSetDimensions(pep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
} }
// print info // print info
// Message::Info("SLEPc solution method: %s", type); // Message::Info("SLEPc solution method: %s", type);
PetscInt nev; PetscInt nev;
_try(PEPGetDimensions(pep, &nev, PETSC_NULL, PETSC_NULL)); _try(PEPGetDimensions(pep, &nev, PETSC_NULL, PETSC_NULL));
Message::Info("SLEPc number of requested eigenvalues: %d", nev); Message::Info("SLEPc number of requested eigenvalues: %d", nev);
PetscReal tol; PetscReal tol;
PetscInt maxit; PetscInt maxit;
_try(PEPGetTolerances(pep, &tol, &maxit)); _try(PEPGetTolerances(pep, &tol, &maxit));
Message::Info("SLEPc stopping condition: tol=%g, maxit=%d", tol, maxit); Message::Info("SLEPc stopping condition: tol=%g, maxit=%d", tol, maxit);
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6)
_try(PEPSetScale(pep, PEP_SCALE_SCALAR, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECI _try(PEPSetScale(pep, PEP_SCALE_SCALAR, PETSC_DECIDE, PETSC_DECIDE,
DE)); PETSC_DECIDE));
#else #else
_try(PEPSetScale(pep, PEP_SCALE_SCALAR, PETSC_DECIDE, PETSC_NULL, PETSC_NULL, _try(PEPSetScale(pep, PEP_SCALE_SCALAR, PETSC_DECIDE, PETSC_NULL, PETSC_NULL,
PETSC_DECIDE, PETSC_DECIDE)); PETSC_DECIDE, PETSC_DECIDE));
#endif #endif
// solve // solve
_try(PEPSolve(pep)); _try(PEPSolve(pep));
// check convergence // check convergence
int its; int its;
skipping to change at line 906 skipping to change at line 931
// get number of converged approximate eigenpairs // get number of converged approximate eigenpairs
PetscInt nconv; PetscInt nconv;
_try(PEPGetConverged(pep, &nconv)); _try(PEPGetConverged(pep, &nconv));
Message::Info("SLEPc number of converged eigenpairs: %d", nconv); Message::Info("SLEPc number of converged eigenpairs: %d", nconv);
// ignore additional eigenvalues if we get more than what we asked // ignore additional eigenvalues if we get more than what we asked
if(nconv > nev) nconv = nev; if(nconv > nev) nconv = nev;
// print eigenvalues and store eigenvectors in DofData // print eigenvalues and store eigenvectors in DofData
_storeEigenVectors(DofData_P, nconv, PETSC_NULL, pep, PETSC_NULL, filterExpres _storeEigenVectors(DofData_P, nconv, PETSC_NULL, pep, PETSC_NULL,
sionIndex); filterExpressionIndex);
_try(PEPDestroy(&pep)); _try(PEPDestroy(&pep));
} }
#endif #endif
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
static PetscErrorCode _myNepMonitor(NEP nep, int its, int nconv, PetscScalar *ei static PetscErrorCode _myNepMonitor(NEP nep, int its, int nconv,
gr, PetscScalar *eigr, PetscScalar *eigi,
PetscScalar *eigi, PetscReal* errest, int ne PetscReal *errest, int nest, void *mctx)
st,
void *mctx)
{ {
return _myMonitor("NEP", its, nconv, eigr, eigi, errest); return _myMonitor("NEP", its, nconv, eigr, eigi, errest);
} }
static void _rationalEVP(struct DofData * DofData_P, int numEigenValues, static void _rationalEVP(struct DofData *DofData_P, int numEigenValues,
double shift_r, double shift_i, int filterExpressionInd double shift_r, double shift_i,
ex, int filterExpressionIndex, List_T *RationalCoefsNum,
List_T *RationalCoefsNum, List_T *RationalCoefsDen, List_T *RationalCoefsDen,
List_T *ApplyResolventRealFreqs, struct DofData * DofDa List_T *ApplyResolventRealFreqs,
ta_P2) struct DofData *DofData_P2)
{ {
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 8) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 8)
NEP nep; NEP nep;
NEPType type; NEPType type;
int max_Nchar = 1000; int max_Nchar = 1000;
char str_coefsNum[6][max_Nchar]; char str_coefsNum[6][max_Nchar];
char str_coefsDen[6][max_Nchar]; char str_coefsDen[6][max_Nchar];
char str_buff[50]; char str_buff[50];
int NumOperators = 2; int NumOperators = 2;
int Flag_ApplyResolvent=0; int Flag_ApplyResolvent = 0;
int nbApplyResolventRealFreqs = (int)List_Nbr(ApplyResolventRealFreqs); int nbApplyResolventRealFreqs = (int)List_Nbr(ApplyResolventRealFreqs);
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
PetscScalar shift = shift_r + PETSC_i * shift_i; PetscScalar shift = shift_r + PETSC_i * shift_i;
#else #else
Message::Warning("Discarding imaginary part in shift"); Message::Warning("Discarding imaginary part in shift");
PetscScalar shift = shift_r; PetscScalar shift = shift_r;
#endif #endif
if(List_Nbr(RationalCoefsNum) > 6 || List_Nbr(RationalCoefsDen) > 6){ if(List_Nbr(RationalCoefsNum) > 6 || List_Nbr(RationalCoefsDen) > 6) {
Message::Error("Too many rational terms in nonlinear eigenvalue problem (max Message::Error(
= 6)"); "Too many rational terms in nonlinear eigenvalue problem (max = 6)");
return; return;
} }
if (nbApplyResolventRealFreqs>0){ if(nbApplyResolventRealFreqs > 0) {
Flag_ApplyResolvent = 1; Flag_ApplyResolvent = 1;
if (nbApplyResolventRealFreqs!=DofData_P2->CounterOfRHS) if(nbApplyResolventRealFreqs != DofData_P2->CounterOfRHS)
Message::Error("Please provide a number of RHS terms equal to the number o Message::Error("Please provide a number of RHS terms equal to the number "
f real pulsations"); "of real pulsations");
} }
std::vector<PetscScalar> tabCoefsNum[6], tabCoefsDen[6]; std::vector<PetscScalar> tabCoefsNum[6], tabCoefsDen[6];
for(int i = 0; i < List_Nbr(RationalCoefsNum); i++){ for(int i = 0; i < List_Nbr(RationalCoefsNum); i++) {
List_T *coefs; List_T *coefs;
List_Read(RationalCoefsNum, i, &coefs); List_Read(RationalCoefsNum, i, &coefs);
for(int j = 0; j < List_Nbr(coefs); j++){ for(int j = 0; j < List_Nbr(coefs); j++) {
double c; List_Read(coefs, j, &c); double c;
List_Read(coefs, j, &c);
tabCoefsNum[i].push_back(c); tabCoefsNum[i].push_back(c);
} }
} }
for(int i = 0; i < List_Nbr(RationalCoefsDen); i++){ for(int i = 0; i < List_Nbr(RationalCoefsDen); i++) {
List_T *coefs; List_T *coefs;
List_Read(RationalCoefsDen, i, &coefs); List_Read(RationalCoefsDen, i, &coefs);
for(int j = 0; j < List_Nbr(coefs); j++){ for(int j = 0; j < List_Nbr(coefs); j++) {
double c; List_Read(coefs, j, &c); double c;
List_Read(coefs, j, &c);
tabCoefsDen[i].push_back(c); tabCoefsDen[i].push_back(c);
} }
} }
std::vector<PetscReal> tabApplyResolventRealFreqs; std::vector<PetscReal> tabApplyResolventRealFreqs;
if (Flag_ApplyResolvent){ if(Flag_ApplyResolvent) {
for(int i = 0; i < List_Nbr(ApplyResolventRealFreqs); i++){ for(int i = 0; i < List_Nbr(ApplyResolventRealFreqs); i++) {
double c; List_Read(ApplyResolventRealFreqs, i, &c); double c;
tabApplyResolventRealFreqs.push_back(c); List_Read(ApplyResolventRealFreqs, i, &c);
tabApplyResolventRealFreqs.push_back(c);
} }
} }
_try(NEPCreate(PETSC_COMM_WORLD, &nep)); _try(NEPCreate(PETSC_COMM_WORLD, &nep));
// NLEig1Dof and NLEig2Dof // NLEig1Dof and NLEig2Dof
if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] && if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] &&
!DofData_P->Flag_Init[5] && !DofData_P->Flag_Init[4] && !DofData_P->Flag_Init[5] && !DofData_P->Flag_Init[4] &&
!DofData_P->Flag_Init[3] && !DofData_P->Flag_Init[2]){ !DofData_P->Flag_Init[3] && !DofData_P->Flag_Init[2]) {
NumOperators = 2; NumOperators = 2;
Mat A[2] = {DofData_P->M7.M, DofData_P->M6.M}; Mat A[2] = {DofData_P->M7.M, DofData_P->M6.M};
FN funs[2]; FN funs[2];
for(int i=0;i<2;i++){ for(int i = 0; i < 2; i++) {
_try(FNCreate(PETSC_COMM_WORLD, &funs[i])); _try(FNCreate(PETSC_COMM_WORLD, &funs[i]));
_try(FNSetType(funs[i], FNRATIONAL)); _try(FNSetType(funs[i], FNRATIONAL));
_try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(), &tabCoefsNum[i _try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(),
][0])); &tabCoefsNum[i][0]));
_try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(), &tabCoefsDen _try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(),
[i][0])); &tabCoefsDen[i][0]));
} }
_try(NEPSetSplitOperator(nep,2,A,funs,SUBSET_NONZERO_PATTERN)); _try(NEPSetSplitOperator(nep, 2, A, funs, SUBSET_NONZERO_PATTERN));
} }
// NLEig1Dof, NLEig2Dof and NLEig3Dof // NLEig1Dof, NLEig2Dof and NLEig3Dof
else if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] && else if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] &&
DofData_P->Flag_Init[5] && !DofData_P->Flag_Init[4] && DofData_P->Flag_Init[5] && !DofData_P->Flag_Init[4] &&
!DofData_P->Flag_Init[3] && !DofData_P->Flag_Init[2]){ !DofData_P->Flag_Init[3] && !DofData_P->Flag_Init[2]) {
NumOperators = 3; NumOperators = 3;
Mat A[3] = {DofData_P->M7.M, DofData_P->M6.M, DofData_P->M5.M}; Mat A[3] = {DofData_P->M7.M, DofData_P->M6.M, DofData_P->M5.M};
FN funs[3]; FN funs[3];
for(int i=0;i<3;i++){ for(int i = 0; i < 3; i++) {
_try(FNCreate(PETSC_COMM_WORLD,&funs[i])); _try(FNCreate(PETSC_COMM_WORLD, &funs[i]));
_try(FNSetType(funs[i],FNRATIONAL)); _try(FNSetType(funs[i], FNRATIONAL));
_try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(), &tabCoefsNum[i _try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(),
][0])); &tabCoefsNum[i][0]));
_try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(), &tabCoefsDen _try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(),
[i][0])); &tabCoefsDen[i][0]));
} }
_try(NEPSetSplitOperator(nep,3,A,funs,SUBSET_NONZERO_PATTERN)); _try(NEPSetSplitOperator(nep, 3, A, funs, SUBSET_NONZERO_PATTERN));
} }
// NLEig1Dof, NLEig2Dof, NLEig3Dof and NLEig4Dof // NLEig1Dof, NLEig2Dof, NLEig3Dof and NLEig4Dof
else if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] && else if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] &&
DofData_P->Flag_Init[5] && DofData_P->Flag_Init[4] && DofData_P->Flag_Init[5] && DofData_P->Flag_Init[4] &&
!DofData_P->Flag_Init[3] && !DofData_P->Flag_Init[2]){ !DofData_P->Flag_Init[3] && !DofData_P->Flag_Init[2]) {
NumOperators = 4; NumOperators = 4;
Mat A[4] = {DofData_P->M7.M, DofData_P->M6.M, DofData_P->M5.M, Mat A[4] = {DofData_P->M7.M, DofData_P->M6.M, DofData_P->M5.M,
DofData_P->M4.M}; DofData_P->M4.M};
FN funs[4]; FN funs[4];
for(int i=0;i<4;i++){ for(int i = 0; i < 4; i++) {
_try(FNCreate(PETSC_COMM_WORLD,&funs[i])); _try(FNCreate(PETSC_COMM_WORLD, &funs[i]));
_try(FNSetType(funs[i],FNRATIONAL)); _try(FNSetType(funs[i], FNRATIONAL));
_try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(), &tabCoefsNum[i _try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(),
][0])); &tabCoefsNum[i][0]));
_try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(), &tabCoefsDen _try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(),
[i][0])); &tabCoefsDen[i][0]));
} }
_try(NEPSetSplitOperator(nep,4,A,funs,SUBSET_NONZERO_PATTERN)); _try(NEPSetSplitOperator(nep, 4, A, funs, SUBSET_NONZERO_PATTERN));
} }
// NLEig1Dof, NLEig2Dof, NLEig3Dof, NLEig4Dof and NLEig5Dof // NLEig1Dof, NLEig2Dof, NLEig3Dof, NLEig4Dof and NLEig5Dof
else if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] && else if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] &&
DofData_P->Flag_Init[5] && DofData_P->Flag_Init[4] && DofData_P->Flag_Init[5] && DofData_P->Flag_Init[4] &&
DofData_P->Flag_Init[3] && !DofData_P->Flag_Init[2]){ DofData_P->Flag_Init[3] && !DofData_P->Flag_Init[2]) {
NumOperators = 5; NumOperators = 5;
Mat A[5] = {DofData_P->M7.M, DofData_P->M6.M, DofData_P->M5.M, Mat A[5] = {DofData_P->M7.M, DofData_P->M6.M, DofData_P->M5.M,
DofData_P->M4.M, DofData_P->M3.M}; DofData_P->M4.M, DofData_P->M3.M};
FN funs[5]; FN funs[5];
for(int i=0;i<5;i++){ for(int i = 0; i < 5; i++) {
_try(FNCreate(PETSC_COMM_WORLD,&funs[i])); _try(FNCreate(PETSC_COMM_WORLD, &funs[i]));
_try(FNSetType(funs[i],FNRATIONAL)); _try(FNSetType(funs[i], FNRATIONAL));
_try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(), &tabCoefsNum[i _try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(),
][0])); &tabCoefsNum[i][0]));
_try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(), &tabCoefsDen _try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(),
[i][0])); &tabCoefsDen[i][0]));
} }
_try(NEPSetSplitOperator(nep,5,A,funs,SUBSET_NONZERO_PATTERN)); _try(NEPSetSplitOperator(nep, 5, A, funs, SUBSET_NONZERO_PATTERN));
} }
// NLEig1Dof, NLEig2Dof, NLEig3Dof, NLEig4Dof, NLEig5Dof and NLEig6Dof // NLEig1Dof, NLEig2Dof, NLEig3Dof, NLEig4Dof, NLEig5Dof and NLEig6Dof
else if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] && else if(DofData_P->Flag_Init[7] && DofData_P->Flag_Init[6] &&
DofData_P->Flag_Init[5] && DofData_P->Flag_Init[4] && DofData_P->Flag_Init[5] && DofData_P->Flag_Init[4] &&
DofData_P->Flag_Init[3] && DofData_P->Flag_Init[2]){ DofData_P->Flag_Init[3] && DofData_P->Flag_Init[2]) {
NumOperators = 6; NumOperators = 6;
Mat A[6] = {DofData_P->M7.M, DofData_P->M6.M, DofData_P->M5.M, Mat A[6] = {DofData_P->M7.M, DofData_P->M6.M, DofData_P->M5.M,
DofData_P->M4.M, DofData_P->M3.M, DofData_P->M2.M}; DofData_P->M4.M, DofData_P->M3.M, DofData_P->M2.M};
FN funs[6]; FN funs[6];
for(int i=0;i<6;i++){ for(int i = 0; i < 6; i++) {
_try(FNCreate(PETSC_COMM_WORLD,&funs[i])); _try(FNCreate(PETSC_COMM_WORLD, &funs[i]));
_try(FNSetType(funs[i],FNRATIONAL)); _try(FNSetType(funs[i], FNRATIONAL));
_try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(), &tabCoefsNum[i _try(FNRationalSetNumerator(funs[i], tabCoefsNum[i].size(),
][0])); &tabCoefsNum[i][0]));
_try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(), &tabCoefsDen _try(FNRationalSetDenominator(funs[i], tabCoefsDen[i].size(),
[i][0])); &tabCoefsDen[i][0]));
} }
_try(NEPSetSplitOperator(nep,6,A,funs,SUBSET_NONZERO_PATTERN)); _try(NEPSetSplitOperator(nep, 6, A, funs, SUBSET_NONZERO_PATTERN));
} }
else{ else {
Message::Info("Illegal use of Rational..."); Message::Info("Illegal use of Rational...");
} }
Message::Info("Solving non-linear eigenvalue problem using slepc NEP"); Message::Info("Solving non-linear eigenvalue problem using slepc NEP");
Message::Info("Number of Operators %d - Setting %d Rational functions :", Message::Info("Number of Operators %d - Setting %d Rational functions :",
NumOperators, NumOperators); NumOperators, NumOperators);
for(int k = 0; k < NumOperators; k++){ for(int k = 0; k < NumOperators; k++) {
sprintf(str_coefsNum[k],"num%d(iw)=", k + 1); sprintf(str_coefsNum[k], "num%d(iw)=", k + 1);
sprintf(str_coefsDen[k],"den%d(iw)=", k + 1); sprintf(str_coefsDen[k], "den%d(iw)=", k + 1);
for(unsigned int i = 0; i < tabCoefsNum[k].size() - 1; i++){ for(unsigned int i = 0; i < tabCoefsNum[k].size() - 1; i++) {
sprintf(str_buff," (%+.2e)*(iw)^%d +", sprintf(str_buff, " (%+.2e)*(iw)^%d +", PetscRealPart(tabCoefsNum[k][i]),
PetscRealPart(tabCoefsNum[k][i]), int(tabCoefsNum[k].size()-1-i)); int(tabCoefsNum[k].size() - 1 - i));
strcat(str_coefsNum[k],str_buff); strcat(str_coefsNum[k], str_buff);
} }
sprintf(str_buff," (%+.2e)", sprintf(str_buff, " (%+.2e)",
PetscRealPart(tabCoefsNum[k][tabCoefsNum[k].size()-1])); PetscRealPart(tabCoefsNum[k][tabCoefsNum[k].size() - 1]));
strcat(str_coefsNum[k],str_buff); strcat(str_coefsNum[k], str_buff);
for(unsigned int i = 0; i < tabCoefsDen[k].size() - 1; i++){ for(unsigned int i = 0; i < tabCoefsDen[k].size() - 1; i++) {
sprintf(str_buff," (%+.2e)*(iw)^%d +", sprintf(str_buff, " (%+.2e)*(iw)^%d +", PetscRealPart(tabCoefsDen[k][i]),
PetscRealPart(tabCoefsDen[k][i]), int(tabCoefsDen[k].size()-1-i)); int(tabCoefsDen[k].size() - 1 - i));
strcat(str_coefsDen[k],str_buff); strcat(str_coefsDen[k], str_buff);
} }
sprintf(str_buff," (%+.2e)", sprintf(str_buff, " (%+.2e)",
PetscRealPart(tabCoefsDen[k][tabCoefsDen[k].size() - 1])); PetscRealPart(tabCoefsDen[k][tabCoefsDen[k].size() - 1]));
strcat(str_coefsDen[k], str_buff); strcat(str_coefsDen[k], str_buff);
Message::Info(str_coefsNum[k]); Message::Info(str_coefsNum[k]);
Message::Info(str_coefsDen[k]); Message::Info(str_coefsDen[k]);
} }
_try(NEPSetDimensions(nep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); _try(NEPSetDimensions(nep, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
_try(NEPSetTolerances(nep, 1.e-6, PETSC_DEFAULT)); _try(NEPSetTolerances(nep, 1.e-6, PETSC_DEFAULT));
_try(NEPSetType(nep, NEPNLEIGS)); _try(NEPSetType(nep, NEPNLEIGS));
_try(NEPSetProblemType(nep, NEP_RATIONAL)); _try(NEPSetProblemType(nep, NEP_RATIONAL));
_try(NEPSetWhichEigenpairs(nep, NEP_TARGET_MAGNITUDE)); _try(NEPSetWhichEigenpairs(nep, NEP_TARGET_MAGNITUDE));
_try(NEPMonitorSet(nep, _myNepMonitor, PETSC_NULL, PETSC_NULL)); _try(NEPMonitorSet(nep, _myNepMonitor, PETSC_NULL, PETSC_NULL));
_try(NEPSetTarget(nep, shift)); _try(NEPSetTarget(nep, shift));
_try(NEPSetConvergenceTest(nep, NEP_CONV_REL)); _try(NEPSetConvergenceTest(nep, NEP_CONV_REL));
if(Flag_ApplyResolvent){ if(Flag_ApplyResolvent) {
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 12) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 12)
Message::Info("Using full basis variant (required for ApplyResolvent)"); Message::Info("Using full basis variant (required for ApplyResolvent)");
Message::Info("Using two sided variant (required for ApplyResolvent)"); Message::Info("Using two sided variant (required for ApplyResolvent)");
_try(NEPNLEIGSSetFullBasis(nep,PETSC_TRUE)); _try(NEPNLEIGSSetFullBasis(nep, PETSC_TRUE));
_try(NEPNLEIGSSetRestart(nep,0.5)); _try(NEPNLEIGSSetRestart(nep, 0.5));
_try(NEPSetTwoSided(nep,PETSC_TRUE)); _try(NEPSetTwoSided(nep, PETSC_TRUE));
#else #else
Message::Error("SLEPC >= 3.12 required for ApplyResolvant"); Message::Error("SLEPC >= 3.12 required for ApplyResolvant");
#endif #endif
} }
_try(NEPSetFromOptions(nep)); _try(NEPSetFromOptions(nep));
// print info // print info
_try(NEPGetType(nep, &type)); _try(NEPGetType(nep, &type));
Message::Info("SLEPc solution method: %s", type); Message::Info("SLEPc solution method: %s", type);
skipping to change at line 1154 skipping to change at line 1197
// TODO : we could also check left residual here // TODO : we could also check left residual here
// get number of converged approximate eigenpairs // get number of converged approximate eigenpairs
PetscInt nconv; PetscInt nconv;
_try(NEPGetConverged(nep, &nconv)); _try(NEPGetConverged(nep, &nconv));
Message::Info("SLEPc number of converged eigenpairs: %d", nconv); Message::Info("SLEPc number of converged eigenpairs: %d", nconv);
// ignore additional eigenvalues if we get more than what we asked // ignore additional eigenvalues if we get more than what we asked
if(nconv > nev) nconv = nev; if(nconv > nev) nconv = nev;
_storeEigenVectors(DofData_P, nconv, PETSC_NULL, PETSC_NULL, nep, filterExpres _storeEigenVectors(DofData_P, nconv, PETSC_NULL, PETSC_NULL, nep,
sionIndex); filterExpressionIndex);
// TODO ? : we could print the left eigenvectors as well // TODO ? : we could print the left eigenvectors as well
PetscComplex Lambda; PetscComplex Lambda;
// Vec VRHS; // Vec VRHS;
std::vector<Vec> ListOfExpansionResults; std::vector<Vec> ListOfExpansionResults;
// char fname[100]; // char fname[100];
// static PetscViewer myviewer; // static PetscViewer myviewer;
if(Flag_ApplyResolvent){ if(Flag_ApplyResolvent) {
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 12) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 12)
Message::Info("A RHS term is available for ApplyResolvent!"); Message::Info("A RHS term is available for ApplyResolvent!");
ListOfExpansionResults.reserve(nbApplyResolventRealFreqs); ListOfExpansionResults.reserve(nbApplyResolventRealFreqs);
for(int i = 0; i < nbApplyResolventRealFreqs; i++){ for(int i = 0; i < nbApplyResolventRealFreqs; i++) {
// // debug : print RHS in matlab format // // debug : print RHS in matlab format
// sprintf(fname,"applyresolvent_vec_input_%03d.m",i); // sprintf(fname,"applyresolvent_vec_input_%03d.m",i);
// PetscViewerASCIIOpen(PETSC_COMM_WORLD, fname , &myviewer); // PetscViewerASCIIOpen(PETSC_COMM_WORLD, fname , &myviewer);
// PetscViewerPushFormat(myviewer, PETSC_VIEWER_ASCII_MATLAB); // PetscViewerPushFormat(myviewer, PETSC_VIEWER_ASCII_MATLAB);
// VecView(DofData_P2->ListOfRHS[i].V,myviewer); // VecView(DofData_P2->ListOfRHS[i].V,myviewer);
// PetscViewerPopFormat(myviewer); // PetscViewerPopFormat(myviewer);
_try(MatCreateVecs(DofData_P->M7.M,PETSC_NULL,&ListOfExpansionResults[i])) _try(
; MatCreateVecs(DofData_P->M7.M, PETSC_NULL, &ListOfExpansionResults[i]));
Lambda = PETSC_i*tabApplyResolventRealFreqs[i]; Lambda = PETSC_i * tabApplyResolventRealFreqs[i];
Message::Info("Applying Resolvent with real angular frequency : %f",tabApp Message::Info("Applying Resolvent with real angular frequency : %f",
lyResolventRealFreqs[i]); tabApplyResolventRealFreqs[i]);
_try(NEPApplyResolvent(nep,NULL,Lambda, _try(NEPApplyResolvent(nep, NULL, Lambda, DofData_P2->ListOfRHS[i].V,
DofData_P2->ListOfRHS[i].V,ListOfExpansionResults[i])); ListOfExpansionResults[i]));
} }
_storeExpansion(DofData_P, nbApplyResolventRealFreqs, tabApplyResolventRealF _storeExpansion(DofData_P, nbApplyResolventRealFreqs,
reqs, ListOfExpansionResults); tabApplyResolventRealFreqs, ListOfExpansionResults);
#else #else
Message::Error("SLEPC >= 3.12 required for ApplyResolvant"); Message::Error("SLEPC >= 3.12 required for ApplyResolvant");
#endif #endif
} }
// _try(PetscViewerDestroy(&myviewer)); // _try(PetscViewerDestroy(&myviewer));
_try(NEPDestroy(&nep)); _try(NEPDestroy(&nep));
#else #else
Message::Error("Nonlinear eigenvalue solver requires PETSc/SLEPc >= 3.8"); Message::Error("Nonlinear eigenvalue solver requires PETSc/SLEPc >= 3.8");
#endif #endif
} }
#endif #endif
void EigenSolve_SLEPC(struct DofData * DofData_P, int numEigenValues, void EigenSolve_SLEPC(struct DofData *DofData_P, int numEigenValues,
double shift_r, double shift_i, int FilterExpressionIndex, double shift_r, double shift_i, int FilterExpressionIndex,
List_T *RationalCoefsNum, List_T *RationalCoefsDen, List_T *RationalCoefsNum, List_T *RationalCoefsDen,
List_T *ApplyResolventRealFreqs, struct DofData * DofData_ List_T *ApplyResolventRealFreqs,
P2) struct DofData *DofData_P2)
{ {
// Warn if we are not in harmonic regime (we won't be able to compute/store // Warn if we are not in harmonic regime (we won't be able to compute/store
// complex eigenvectors). // complex eigenvectors).
if(Current.NbrHar != 2){ if(Current.NbrHar != 2) {
Message::Info("EigenSolve will only store the real part of the eigenvectors; Message::Info(
" "EigenSolve will only store the real part of the eigenvectors; "
"Define the system with \"Type Complex\" if this is an issue") "Define the system with \"Type Complex\" if this is an issue");
;
} }
#if !defined(PETSC_USE_COMPLEX) #if !defined(PETSC_USE_COMPLEX)
if(Current.NbrHar == 2){ if(Current.NbrHar == 2) {
Message::Warning("Using PETSc in real arithmetic for complex-simulated-real Message::Warning(
matrices"); "Using PETSc in real arithmetic for complex-simulated-real matrices");
} }
#endif #endif
// GenerateSeparate[] can create up to six matrices M6, M5, M4, M3, M2, M1 suc // GenerateSeparate[] can create up to six matrices M6, M5, M4, M3, M2, M1
h that // such that i*w^5 M6 x + w^4 M5 x + -iw^3 M4 x + -w^2 M3 x + iw M2 x + M1 x =
// i*w^5 M6 x + w^4 M5 x + -iw^3 M4 x + -w^2 M3 x + iw M2 x + M1 x = 0 // 0 check Flag_Init[i] to see which operators exist.
// check Flag_Init[i] to see which operators exist. if(!DofData_P->Flag_Init[7]) {
if(!DofData_P->Flag_Init[7]){ if(!DofData_P->Flag_Init[1] || !DofData_P->Flag_Init[3]) {
if(!DofData_P->Flag_Init[1] || !DofData_P->Flag_Init[3]){ Message::Error("No System available for EigenSolve: check 'DtDt' and "
Message::Error("No System available for EigenSolve: check 'DtDt' and 'Gene "'GenerateSeparate'");
rateSeparate'");
return; return;
} }
if(!DofData_P->Flag_Init[4]&& !DofData_P->Flag_Init[5]&& !DofData_P->Flag_In if(!DofData_P->Flag_Init[4] && !DofData_P->Flag_Init[5] &&
it[6]){ !DofData_P->Flag_Init[6]) {
if(!DofData_P->Flag_Init[2]){ if(!DofData_P->Flag_Init[2]) {
// the shift refers to w^2 // the shift refers to w^2
_linearEVP(DofData_P, numEigenValues, shift_r, shift_i, FilterExpression _linearEVP(DofData_P, numEigenValues, shift_r, shift_i,
Index); FilterExpressionIndex);
} }
else{ else {
// the shift refers to w // the shift refers to w
_quadraticEVP(DofData_P, numEigenValues, shift_r, shift_i, FilterExpress _quadraticEVP(DofData_P, numEigenValues, shift_r, shift_i,
ionIndex); FilterExpressionIndex);
} }
} }
else{ else {
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5)
Message::Error("Please upgrade to slepc >= 3.5.1 for polynomial EVP suppor Message::Error(
t!"); "Please upgrade to slepc >= 3.5.1 for polynomial EVP support!");
return; return;
#else #else
// the shift refers to w // the shift refers to w
_polynomialEVP(DofData_P, numEigenValues, shift_r, shift_i, FilterExpressi _polynomialEVP(DofData_P, numEigenValues, shift_r, shift_i,
onIndex); FilterExpressionIndex);
#endif #endif
} }
} }
else{ else {
#if defined(PETSC_USE_COMPLEX) #if defined(PETSC_USE_COMPLEX)
#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 7) #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 7)
Message::Error("Please upgrade to slepc >= 3.7.3 for non-linear EVP support! Message::Error(
"); "Please upgrade to slepc >= 3.7.3 for non-linear EVP support!");
return; return;
#else #else
_rationalEVP(DofData_P, numEigenValues, shift_r, shift_i, FilterExpressionIn _rationalEVP(DofData_P, numEigenValues, shift_r, shift_i,
dex, FilterExpressionIndex, RationalCoefsNum, RationalCoefsDen,
RationalCoefsNum, RationalCoefsDen, ApplyResolventRealFreqs, Do ApplyResolventRealFreqs, DofData_P2);
fData_P2);
#endif #endif
#else #else
Message::Error("Please compile Petsc/Slepc with complex arithmetic for non l Message::Error("Please compile Petsc/Slepc with complex arithmetic for non "
inear EVP support!"); "linear EVP support!");
#endif #endif
} }
} }
#endif #endif
 End of changes. 151 change blocks. 
377 lines changed or deleted 365 lines changed or added
To the C++ Programs section of the getdp source changes report or the top of this page
Mainly generated by pkgdiff using rfcdiff
Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)