[petsc-users] Matrix free generalized eigenvalue problem
Миша Сальников
rolekst095 at gmail.com
Sun May 15 07:39:08 CDT 2022
Hi, I am developing a program for solving generalized eigenvalue problem Ax
= k Bx when both operators A and B are matrix-free. I found some examples
on SLEPc documentation but and there was and example where problem Ax = k x
are solving matrix free. So I tried to change this example a little bit and
added operator B as identity matrix-free operator. But for now I am getting
an error. My source code is attached. And the error message is below
Can you please explain to me what I am doing wrong?
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: See https://petsc.org/release/overview/linear_solve_table/
for possible LU and Cholesky solvers
[0]PETSC ERROR: Could not locate a solver type for factorization type LU
and matrix type shell.
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.15.5, Sep 29, 2021
[0]PETSC ERROR: ./exe_eps_impl_matr on a named LAPTOP-7DN2DH7N by misha
Sun May 15 15:21:59 2022
[0]PETSC ERROR: Configure options --build=x86_64-linux-gnu --prefix=/usr
--includedir=${prefix}/include --mandir=${prefix}/share/man
--infodir=${prefix}/share/info --sysconfdir=/etc --localstatedir=/var
--with-option-checking=0 --with-silent-rules=0
--libdir=${prefix}/lib/x86_64-linux-gnu --runstatedir=/run
--with-maintainer-mode=0 --with-dependency-tracking=0 --with-64-bit-indices
--with-debugging=0 --shared-library-extension=64_real
--with-shared-libraries --with-pic=1 --with-cc=mpicc --with-cxx=mpicxx
--with-fc=mpif90 --with-cxx-dialect=C++11 --with-opencl=1
--with-blas-lib=-lblas --with-lapack-lib=-llapack --with-scalapack=1
--with-scalapack-lib=-lscalapack-openmpi --with-ptscotch=1
--with-ptscotch-include=/usr/include/scotch
--with-ptscotch-lib="-lptesmumps -lptscotch -lptscotcherr" --with-fftw=1
--with-fftw-include="[]" --with-fftw-lib="-lfftw3 -lfftw3_mpi"
--with-yaml=1 --with-hdf5-include=/usr/include/hdf5/openmpi
--with-hdf5-lib="-L/usr/lib/x86_64-linux-gnu/hdf5/openmpi
-L/usr/lib/x86_64-linux-gnu/openmpi/lib -lhdf5 -lmpi"
--CXX_LINKER_FLAGS=-Wl,--no-as-needed --with-hypre=1
--with-hypre-include=/usr/include/hypre64m --with-hypre-lib=-lHYPRE64m
--with-mumps=1 --with-mumps-include="[]" --with-mumps-lib="-ldmumps_64
-lzmumps_64 -lsmumps_64 -lcmumps_64 -lmumps_common_64 -lpord_64"
--with-suitesparse=1 --with-suitesparse-include=/usr/include/suitesparse
--with-suitesparse-lib="-lumfpack -lamd -lcholmod -lklu" --with-zoltan=1
--with-zoltan-include=/usr/include/trilinos
--with-zoltan-lib=-ltrilinos_zoltan
--prefix=/usr/lib/petscdir/petsc64-3.15/x86_64-linux-gnu-real
--PETSC_ARCH=x86_64-linux-gnu-real-64 CFLAGS="-g -O2
-ffile-prefix-map=/build/petsc-tFq6Xk/petsc-3.15.5+dfsg1=. -flto=auto
-ffat-lto-objects -flto=auto -ffat-lto-objects -fstack-protector-strong
-Wformat -Werror=format-security -fPIC" CXXFLAGS="-g -O2
-ffile-prefix-map=/build/petsc-tFq6Xk/petsc-3.15.5+dfsg1=. -flto=auto
-ffat-lto-objects -flto=auto -ffat-lto-objects -fstack-protector-strong
-Wformat -Werror=format-security -fPIC" FCFLAGS="-g -O2
-ffile-prefix-map=/build/petsc-tFq6Xk/petsc-3.15.5+dfsg1=. -flto=auto
-ffat-lto-objects -flto=auto -ffat-lto-objects -fstack-protector-strong
-fPIC -ffree-line-length-0" FFLAGS="-g -O2
-ffile-prefix-map=/build/petsc-tFq6Xk/petsc-3.15.5+dfsg1=. -flto=auto
-ffat-lto-objects -flto=auto -ffat-lto-objects -fstack-protector-strong
-fPIC -ffree-line-length-0" CPPFLAGS="-Wdate-time -D_FORTIFY_SOURCE=2"
LDFLAGS="-Wl,-Bsymbolic-functions -flto=auto -ffat-lto-objects -flto=auto
-Wl,-z,relro -fPIC" MAKEFLAGS=w
[0]PETSC ERROR: #1 MatGetFactor() at ./src/mat/interface/matrix.c:4756
[0]PETSC ERROR: #2 PCSetUp_LU() at ./src/ksp/pc/impls/factor/lu/lu.c:82
[0]PETSC ERROR: #3 PCSetUp() at ./src/ksp/pc/interface/precon.c:1017
[0]PETSC ERROR: #4 KSPSetUp() at ./src/ksp/ksp/interface/itfunc.c:406
[0]PETSC ERROR: #5 STSetUp_Shift() at
./src/sys/classes/st/impls/shift/shift.c:107
[0]PETSC ERROR: #6 STSetUp() at ./src/sys/classes/st/interface/stsolve.c:582
[0]PETSC ERROR: #7 EPSSetUp() at ./src/eps/interface/epssetup.c:350
[0]PETSC ERROR: #8 EPSSolve() at ./src/eps/interface/epssolve.c:136
[0]PETSC ERROR: #9 main() at ex_implicit_matr.cpp:95
[0]PETSC ERROR: No PETSc Option Table entries
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to petsc-maint at mcs.anl.gov----------
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 92.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220515/c5f7e545/attachment.html>
-------------- next part --------------
/*
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SLEPc - Scalable Library for Eigenvalue Problem Computations
Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
This file is part of SLEPc.
SLEPc is distributed under a 2-clause BSD license (see LICENSE).
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
// 2-D Laplacian Eigenproblem (matrix-free version), N=100 (10x10 grid)
// Solution method: krylovschur
// Number of requested eigenvalues: 1
// ---------------------- --------------------
// k ||Ax-kx||/||kx||
// ---------------------- --------------------
// 7.837972 1.84054e-09
// 7.601493 6.49209e-09
// ---------------------- --------------------
static char help[] = "Solves the same eigenproblem as in example ex2, but using a shell matrix. "
"The problem is a standard symmetric eigenproblem corresponding to the 2-D Laplacian operator.\n\n"
"The command line options are:\n"
" -n <n>, where <n> = number of grid subdivisions in both x and y dimensions.\n\n";
#include <slepceps.h>
/*
User-defined routines
*/
PetscErrorCode MatMult_Laplacian2D(Mat A, Vec x, Vec y);
PetscErrorCode MatMult_Laplacian2D_right(Mat A, Vec x, Vec y);
PetscErrorCode MatGetDiagonal_Laplacian2D(Mat A, Vec diag);
int main(int argc, char **argv)
{
Mat A; /* operator matrix */
Mat B; /* operator matrix */
EPS eps; /* eigenproblem solver context */
EPSType type;
PetscMPIInt size;
PetscInt N, n = 10, nev;
PetscBool terse;
CHKERRQ(SlepcInitialize(&argc, &argv, (char *)0, help));
CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
CHKERRQ(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
N = n * n;
CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "\n2-D Laplacian Eigenproblem (matrix-free version), N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n", N, n, n));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the operator matrix that defines the eigensystem, Ax=kx
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
CHKERRQ(MatCreateShell(PETSC_COMM_WORLD, N, N, N, N, &n, &A));
CHKERRQ(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))MatMult_Laplacian2D));
CHKERRQ(MatCreateShell(PETSC_COMM_WORLD, N, N, N, N, &n, &B));
CHKERRQ(MatShellSetOperation(B, MATOP_MULT, (void (*)(void))MatMult_Laplacian2D_right));
CHKERRQ(MatSetFromOptions(A));
CHKERRQ(MatSetFromOptions(B));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the eigensolver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create eigensolver context
*/
CHKERRQ(EPSCreate(PETSC_COMM_WORLD, &eps));
/*
Set operators. In this case, it is a standard eigenvalue problem
*/
CHKERRQ(EPSSetOperators(eps, A, B));
// CHKERRQ(EPSSetOperators(eps, A, NULL));
/*
Set solver parameters at runtime
*/
CHKERRQ(EPSSetFromOptions(eps));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the eigensystem
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
CHKERRQ(EPSSolve(eps));
/*
Optional: Get some information from the solver and display it
*/
CHKERRQ(EPSGetType(eps, &type));
CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, " Solution method: %s\n\n", type));
CHKERRQ(EPSGetDimensions(eps, &nev, NULL, NULL));
CHKERRQ(PetscPrintf(PETSC_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 */
CHKERRQ(PetscOptionsHasName(NULL, NULL, "-terse", &terse));
if (terse)
CHKERRQ(EPSErrorView(eps, EPS_ERROR_RELATIVE, NULL));
else
{
CHKERRQ(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
CHKERRQ(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD));
CHKERRQ(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
}
CHKERRQ(EPSDestroy(&eps));
CHKERRQ(MatDestroy(&A));
CHKERRQ(SlepcFinalize());
return 0;
}
/*
Compute the matrix vector multiplication y<---T*x where T is a nx by nx
tridiagonal matrix with DD on the diagonal, DL on the subdiagonal, and
DU on the superdiagonal.
*/
static void tv(int nx, const PetscScalar *x, PetscScalar *y)
{
PetscScalar dd, dl, du;
int j;
dd = 4.0;
dl = -1.0;
du = -1.0;
y[0] = dd * x[0] + du * x[1];
for (j = 1; j < nx - 1; j++)
y[j] = dl * x[j - 1] + dd * x[j] + du * x[j + 1];
y[nx - 1] = dl * x[nx - 2] + dd * x[nx - 1];
}
/*
Matrix-vector product subroutine for the 2D Laplacian.
The matrix used is the 2 dimensional discrete Laplacian on unit square with
zero Dirichlet boundary condition.
Computes y <-- A*x, where A is the block tridiagonal matrix
| T -I |
|-I T -I |
A = | -I T |
| ... -I|
| -I T|
The subroutine TV is called to compute y<--T*x.
*/
PetscErrorCode MatMult_Laplacian2D(Mat A, Vec x, Vec y)
{
void *ctx;
int nx, lo, i, j;
const PetscScalar *px;
PetscScalar *py;
PetscFunctionBeginUser;
CHKERRQ(MatShellGetContext(A, &ctx));
nx = *(int *)ctx;
CHKERRQ(VecGetArrayRead(x, &px));
CHKERRQ(VecGetArray(y, &py));
tv(nx, &px[0], &py[0]);
for (i = 0; i < nx; i++)
py[i] -= px[nx + i];
for (j = 2; j < nx; j++)
{
lo = (j - 1) * nx;
tv(nx, &px[lo], &py[lo]);
for (i = 0; i < nx; i++)
py[lo + i] -= px[lo - nx + i] + px[lo + nx + i];
}
lo = (nx - 1) * nx;
tv(nx, &px[lo], &py[lo]);
for (i = 0; i < nx; i++)
py[lo + i] -= px[lo - nx + i];
CHKERRQ(VecRestoreArrayRead(x, &px));
CHKERRQ(VecRestoreArray(y, &py));
PetscFunctionReturn(0);
}
PetscErrorCode MatMult_Laplacian2D_right(Mat A, Vec x, Vec y)
{
void *ctx;
int nx, lo, i, j;
const PetscScalar *px;
PetscScalar *py;
PetscFunctionBeginUser;
CHKERRQ(MatShellGetContext(A, &ctx));
nx = *(int *)ctx;
CHKERRQ(VecGetArrayRead(x, &px));
CHKERRQ(VecGetArray(y, &py));
for (i = 0; i < nx*nx; i++)
py[i] = px[i];
CHKERRQ(VecRestoreArrayRead(x, &px));
CHKERRQ(VecRestoreArray(y, &py));
PetscFunctionReturn(0);
}
/*TEST
test:
suffix: 1
args: -n 72 -eps_nev 4 -eps_ncv 20 -terse
requires: !single
TEST*/
More information about the petsc-users
mailing list