[petsc-users] Passing array of PETSc Vec's to a function and returning it
Rubem Mondaini
rmondaini at csrc.ac.cn
Wed Jul 22 21:48:54 CDT 2020
Dear Jose,
thank you very much for your prompt reply. I did check before I was not
requesting more eigenvectors than the ones converged.
In fact, to make things more practical, I made a short version of my
application, where I reproduce the problem.
test.c:
################################
#include <slepceps.h>
#include <complex.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
// NOTE: Petsc was compiled with "--with-scalar-type=complex"
static char help[] = "Hermitian Eigenproblem\n\n";
PetscErrorCode foo(Vec **y, PetscInt n, MPI_Comm mpi_comm)
{
PetscErrorCode ierr;
EPS eps; // eigenproblem solver context
Mat A;
PetscInt nnz, nrows, nconv;
PetscScalar *val_sp_A;
PetscInt *i_sp_A, *j_sp_A;
PetscReal norm;
Vec x;
nnz = 6; nrows = 4;
// The matrix looks like:
// [1 0 0 1+i]
// [0 1 0 0 ]
// [A] = [0 0 1 0 ]
// [1-i 0 0 1 ]
// I am doing everything with *one* MPI task for the sake of simplicity
val_sp_A = malloc(nnz*sizeof(PetscScalar));
i_sp_A = malloc( (nrows+1)*sizeof(PetscInt));
j_sp_A = malloc(nnz*sizeof(PetscInt));
// Building the matrix (CSR first)
val_sp_A[0] = 1.0; val_sp_A[1] = 1.0 + 1.0*PETSC_i; val_sp_A[2] =
1.0; val_sp_A[3] = 1.0; val_sp_A[4] = 1.0 -1.0*PETSC_i; val_sp_A[5] = 1.0;
j_sp_A[0] = 0; j_sp_A[1] = 3; j_sp_A[2] = 1; j_sp_A[3] =
2; j_sp_A[4] = 0; j_sp_A[5] = 3;
// zero-based index
i_sp_A[0] = 0; i_sp_A[1] = 2; i_sp_A[2] = 3; i_sp_A[3] = 4;
i_sp_A[4] = 6;
// Building the matrix and vecs
ierr = MatCreate(mpi_comm, &A);
CHKERRQ(ierr);
ierr = MatCreateMPIAIJWithArrays(mpi_comm, nrows, nrows, nrows,
nrows, i_sp_A, j_sp_A, val_sp_A, &A);
ierr = MatSetUp(A); CHKERRQ(ierr);
ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);
ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);
ierr = MatCreateVecs(A, &x, NULL);
CHKERRQ(ierr);
// Solve eigenvec problem
ierr = EPSCreate(mpi_comm, &eps);
CHKERRQ(ierr);
ierr = EPSSetOperators( eps, A, NULL);
CHKERRQ(ierr);
ierr = EPSSetProblemType( eps, EPS_HEP );
CHKERRQ(ierr);
ierr = EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL);
CHKERRQ(ierr);
ierr = EPSSetDimensions(eps, n, 2*n, 2*n);
CHKERRQ(ierr);
ierr = EPSSolve(eps); CHKERRQ(ierr);
ierr = EPSGetConverged(eps, &nconv);
CHKERRQ(ierr);
printf("\nnconv = %lu\n", nconv); // I am getting all four
eigenvalues converged, you can test it.
assert(n < nconv); // This guarantees the number
of backed up eigenvectors *will be smaller than nconv*
// Now doing the backup of the eigenvectors
ierr = VecDuplicateVecs(x, n, y); CHKERRQ(ierr);
for (PetscInt i = 0; i < n ; i++) { // I can guarantee that n <
nconv
ierr = EPSGetEigenvector(eps, i, *y[i], NULL); CHKERRQ(ierr);
ierr = VecNorm(*y[i],NORM_2,&norm); CHKERRQ(ierr); //
this prints out fine for i = 0 (norm = 1)
printf("i = %lu\tnorm = %f\n", i, norm);
}
// Deallocating
ierr = EPSDestroy(&eps); CHKERRQ(ierr);
ierr = VecDestroy(&x); CHKERRQ(ierr);
ierr = MatDestroy(&A); CHKERRQ(ierr);
free(val_sp_A);
free(i_sp_A);
free(j_sp_A);
return ierr;
}
int main(int argc,char **argv)
{
PetscErrorCode ierr;
PetscScalar norm;
Vec *y;
PetscInt n = 2; // Backing up *two* eigenvectors
ierr = SlepcInitialize(&argc,&argv,(char*)0,help); CHKERRQ(ierr);
foo(&y, n, PETSC_COMM_WORLD);
ierr = VecDestroyVecs(n, &y); CHKERRQ(ierr);
return 0;
}
#################################################
Makefile:
#################################################
#Source File name
SRC = test
CC = mpiicc
default: $(SRC)
include ${SLEPC_DIR}/lib/slepc/conf/slepc_common
# include $(SLEPC_DIR)/lib/slepc/conf/slepc_variables
CFLAGS_FOR_ICC = -qopt-report=4 -qopt-report-phase ipo -O3 -g -w2
-std=c99 -qopenmp -DMKL_ILP64 -I$(MKLROOT)/include
INTEL_LIB = -Wl,--start-group
${MKLROOT}/lib/intel64/libmkl_intel_ilp64.a
${MKLROOT}/lib/intel64/libmkl_intel_thread.a
${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -liomp5 -lpthread
-lm -ldl
CFLAGS=$(CFLAGS_FOR_ICC)
$(SRC): $(SRC).o
-${CLINKER} -o $(SRC) $(SRC).o ${SLEPC_EPS_LIB} ${INTEL_LIB}
# # ${RM} *.o
$(SRC).o: $(SRC).c
-${CLINKER} -I${PETSC_DIR}/include
-I${PETSC_DIR}/linux-intel/include -I${SLEPC_DIR}/include
-I${SLEPC_DIR}/linux-intel/include -c $(SRC).c
#################################################
Execution
#################################################
[rmondaini at manager slepc_test]$ ./test
nconv = 4
i = 0 norm = 1.000000
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Wrong type of object: Parameter # 3
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.13.3, Jul 01, 2020
[0]PETSC ERROR: ./test on a linux-intel named manager by rmondaini Thu
Jul 23 10:43:15 2020
[0]PETSC ERROR: Configure options PETSC_ARCH=linux-intel
--with-cc=mpiicc --with-cxx=mpiicpc --with-fc=mpiifort
--with-blaslapack-dir=/opt/intel/compilers_and_libraries_2016/linux/mkl/lib/intel64
--with-mpi-include=/opt/intel/compilers_and_libraries_2016/linux/mpi/intel64/include
--with-mpi-lib=/opt/intel/compilers_and_libraries_2016/linux/mpi/intel64/lib/libmpicxx.a
--with-mpiexec=/opt/intel/compilers_and_libraries_2016/linux/mpi/intel64/bin/mpiexec
--with-scalar-type=complex --with-64-bit-blas-indices
--with-64-bit-indices --download-make --force
[0]PETSC ERROR: #1 EPSGetEigenvector() line 504 in
/home/rmondaini/libraries/slepc-3.13.3/src/eps/interface/epssolve.c
[0]PETSC ERROR: #2 foo() line 66 in test.c
#################################################
So the offending statement is precisely the one to retrieve
eigenvectors, and occurs when i = 1 (as in my previous much larger
application).
Do you have suggestions in how to proceed?
Thank you *very* much!
Best,
Rubem
On 7/22/20 11:10 PM, Jose E. Roman wrote:
> Probably you are requesting more eigenvectors than actually computed. Argument i should be smaller than nconv, see https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSGetEigenvector.html
> Jose
>
>
>> El 22 jul 2020, a las 10:25, rmondaini at csrc.ac.cn escribió:
>>
>> I am trying to pass an array of Vec's in PETSc to a function, modify it internally and retrieve the results. The idea is to copy a handful of eigenvectors from the EPS solver to the main routine. A pseudocode is as follows:
>>
>> #include <slepceps.h>
>>
>> PetscErrorCode foo(Vec **y, int n) {
>>
>> EPS eps; // eigenproblem solver context
>>
>> // ...
>>
>> ierr = MatCreateVecs(H, &x, NULL); CHKERRQ(ierr);
>>
>> ierr = EPSSolve(eps); CHKERRQ(ierr);
>>
>> // ...
>>
>> ierr = VecDuplicateVecs(x, n, y); CHKERRQ(ierr);
>>
>> for (int i = 0; i < n ; i++) { // I can guarantee that n < nconv
>>
>> ierr = EPSGetEigenvector(eps, i, *y[i], NULL); CHKERRQ(ierr); // this breaks for i = 1
>>
>> ierr = VecNorm(*y[i],NORM_2,&norm); CHKERRQ(ierr); // this prints out fine for i = 0 (norm = 1)
>>
>> printf("norm = %f\n", norm);
>>
>> }
>>
>> ierr = EPSDestroy(&eps); CHKERRQ(ierr);
>> ierr = VecDestroy(&x); CHKERRQ(ierr);
>>
>> return ierr;
>>
>> }
>>
>> int main(int argc,char **argv)
>> {
>> PetscErrorCode ierr;
>> PetscScalar norm;
>> Vec *y;
>>
>> foo(&y, 3);
>>
>> ierr = VecDestroyVecs(3, &y); CHKERRQ(ierr);
>>
>> return 0;
>> }
>>
>> Am I making a naive mistake here?
>>
More information about the petsc-users
mailing list