[petsc-users] Passing array of PETSc Vec's to a function and returning it
Rubem Mondaini
rmondaini at csrc.ac.cn
Thu Jul 23 02:31:34 CDT 2020
That's is precisely the issue.
Thank you!
On 7/23/20 3:08 PM, Jose E. Roman wrote:
> 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