Hello,
I am new to Slepc and trying to work out the eigenvalues and eigenvectors of my
Jacobian matrix. I am using a shell matrix to work out the matrix-vector
product and I am using the default Krylov-schur method.
My first attempt was not successful and I got the following errors:
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Missing or incorrect user input
[0]PETSC ERROR: The inner product is not well defined: indefinite matrix
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.17.4, unknown
[0]PETSC ERROR: cfdtest on a arch-debug named ming by feng Mon Aug 22 19:21:41
2022
[0]PETSC ERROR: Configure options --with-cc=mpicc --with-cxx=mpicxx --with-fc=0
PETSC_ARCH=arch-debug
[0]PETSC ERROR: #1 BV_SafeSqrt() at
/home/feng/cfd/slepc-3.17.1/include/slepc/private/bvimpl.h:130
[0]PETSC ERROR: #2 BV_SquareRoot_Default() at
/home/feng/cfd/slepc-3.17.1/include/slepc/private/bvimpl.h:365
[0]PETSC ERROR: #3 BVOrthogonalizeCGS1() at
/home/feng/cfd/slepc-3.17.1/src/sys/classes/bv/interface/bvorthog.c:101
[0]PETSC ERROR: #4 BVOrthogonalizeGS() at
/home/feng/cfd/slepc-3.17.1/src/sys/classes/bv/interface/bvorthog.c:177
[0]PETSC ERROR: #5 BVOrthonormalizeColumn() at
/home/feng/cfd/slepc-3.17.1/src/sys/classes/bv/interface/bvorthog.c:402
[0]PETSC ERROR: #6 BVMatArnoldi() at
/home/feng/cfd/slepc-3.17.1/src/sys/classes/bv/interface/bvkrylov.c:91
[0]PETSC ERROR: #7 EPSSolve_KrylovSchur_Default() at
/home/feng/cfd/slepc-3.17.1/src/eps/impls/krylov/krylovschur/krylovschur.c:261
[0]PETSC ERROR: #8 EPSSolve() at
/home/feng/cfd/slepc-3.17.1/src/eps/interface/epssolve.c:147
[0]PETSC ERROR: #9 slepc_eigen_comp() at domain/cfd/slepc_eigen_solve.cpp:77
Could someone please shine some light on this? I have also attached my code.
The code is part of my big library and I cannot attach the whole code, sorry
about this. but I am happy to provide more information. The attached code has
some arrangements for halo exchange, but for the moment, it assumes it is a
serial run.
Many thanks,
Feng
using namespace std;
# include <iomanip> // std::setprecision
# include <domain/cfd/domain.h>
#ifdef PETSC
PetscErrorCode cFdDomain::slepc_eigen_comp()
{
cout << "slepc_eigen_solve\n";
PetscErrorCode ierr;
Mat slepc_A_mf; /* eigenvalue problem matrix */
Int iq, iqs, iqe;
PetscInt blocksize;
MPI_Comm* A_COMM_WORLD;
PetscInt maxneig=8;
PetscInt nev;
EPS eps; /* eigenproblem solver context */
EPSType type;
PetscViewer viewer;
movegrid();
frame();
weights();
cout << "show vtk file\n";
vtk();
A_COMM_WORLD = dev->getcomm();
dof->range( dev->getrank(), &iqs,&iqe );
//number of ghost cells and their global indexes (in the new numbering)
ig_mat = dof->get_igmat();
ighost = new Int [dof->size()];
ighost_local = new Int [dof->size()];
ilocal = new Int [dof->size()];
nghost=0;
nlocal=0;
for(iq=0; iq<dof->size(); iq++)
{
if(iq>=iqs && iq<iqe)
{
//non-halo cells
ilocal[nlocal] = iq;
nlocal++;
}
else
{
ighost[nghost] = ig_mat[iq];
ighost_local[nghost] = iq;
nghost++;
}
}
blocksize = nv;
cout << nlocal << " " << blocksize << "\n";
//create shell matrix
PetscCall(MatCreateShell(*A_COMM_WORLD, nlocal*blocksize, nlocal*blocksize, PETSC_DETERMINE, PETSC_DETERMINE, this, &slepc_A_mf));
PetscCall(MatShellSetOperation(slepc_A_mf,MATOP_MULT,(void(*)(void))mymult_slepc));
// Create eigensolver context
PetscCall(EPSCreate(*A_COMM_WORLD,&eps));
// Set operators. In this case, it is a standard eigenvalue problem
PetscCall(EPSSetOperators(eps,slepc_A_mf,NULL));
PetscCall(EPSSetProblemType(eps,EPS_NHEP));
// Ask for the rightmost eigenvalues
PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
// Set other solver options at runtime
PetscCall(EPSSetFromOptions(eps));
// Solve the eigensystem
PetscCall(EPSSolve(eps));
/*
Optional: Get some information from the solver and display it
*/
PetscCall(EPSGetType(eps,&type));
PetscCall(PetscPrintf(*A_COMM_WORLD," Solution method: %s\n\n",type));
PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
PetscCall(PetscPrintf(*A_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Display solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* show detailed info unless -terse option is given by user */
PetscCall(PetscViewerASCIIGetStdout(*A_COMM_WORLD,&viewer));
PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
PetscCall(EPSConvergedReasonView(eps,viewer));
PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer));
PetscCall(PetscViewerPopFormat(viewer));
PetscCall(EPSDestroy(&eps));
PetscCall(MatDestroy(&slepc_A_mf));
ierr=0;
return ierr;
}
PetscErrorCode cFdDomain::mymult_slepc(Mat m ,Vec x, Vec y)
{
PetscErrorCode ierr;
void *ctx;
cFdDomain *myctx;
Vec localin;
MatShellGetContext(m, &ctx);
myctx = (cFdDomain*)ctx;
// PetscCall(VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD));
// PetscCall(VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD));
// PetscCall(VecGhostGetLocalForm(x,&localin));
ierr = myctx->petsc_set_fd_dq(x); CHKERRQ(ierr);
ierr = myctx->slepc_fd_rhs2(y); CHKERRQ(ierr);
// VecGhostRestoreLocalForm(x,&localin);
return ierr;
}
PetscErrorCode cFdDomain::slepc_fd_rhs2(Vec var)
{
Int iqs, iqe;
Int iv, iq, jq;
Real delta_rhs;
PetscErrorCode ierr;
PetscScalar *array;
Real eps;
Real dq_norm;
dof->range( dev->getrank(), &iqs,&iqe );
dq_norm = 0;
for(iq=iqs; iq<iqe; iq++)
{
for(iv=0; iv<nv; iv++)
{
dq_norm += dq[iv][iq]*dq[iv][iq];
}
}
dev->gsum(1, &dq_norm);
dq_norm = sqrt(dq_norm/dof->gsize());
//eps = sqrt(1+csv_mag_fd)*1.e-4/dq_norm;
eps = 1.e-7/dq_norm;
cout << eps << "\n";
//compute rhs with unperturbed flow states
setv( 0,dof->size(), nv, 0., rhs );
grad( q,dqdx );
gradb( q,dqdx );
gtrhf( rhs );
gtrhm( rhs );
gtrhv( rhs );
//perturb q, halo and non-halo cells
//convert to primitive variable perturbations
fld->dvar( 0,dof->size(), q, aux, dq, daux );
for(iv=0; iv<nv; iv++)
{
for(iq=0; iq<dof->size(); iq++)
{
q[iv][iq] += daux[iv][iq]*eps;
}
}
//compute rhs with perturbed flow states
setv( 0,dof->size(), nv, 0., wrkq );
grad( q,dqdx );
gradb( q,dqdx );
gtrhf( wrkq );
gtrhm( wrkq );
gtrhv( wrkq );
//matrix-vector product
ierr = VecGetArray(var,&array); CHKERRQ(ierr);
for(iv=0; iv<nv; iv++)
{
for(iq=0; iq<nlocal; iq++)
{
jq = ilocal[iq];
array[iv + nv*iq] = (wrkq[iv][jq] - rhs[iv][jq])/eps;
}
}
ierr = VecRestoreArray(var,&array); CHKERRQ(ierr);
///reset q
for(iv=0; iv<nv; iv++)
{
for(iq=0; iq<dof->size(); iq++)
{
q[iv][iq] -= daux[iv][iq]*eps;
}
}
ierr = 0;
return ierr;
}
#endif
