[petsc-users] Passing array of PETSc Vec's to a function and returning it
Jose E. Roman
jroman at dsic.upv.es
Thu Jul 23 02:08:50 CDT 2020
It's a problem with double pointers. You should write (*y)[i] instead of *y[i].
Jose
> El 23 jul 2020, a las 4:48, Rubem Mondaini <rmondaini at csrc.ac.cn> escribió:
>
> 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