[petsc-users] Getting a matrix into PETSc

[Barry Smith]

On Sep 27, 2010, at 5:10 PM, Alexis Morris wrote:

> Hi everyone,
> 
>   I am new to PETSc, and am really excited about using it in combination with SLEPc to find the eigenvalues of a large hermitian sparse matrix.  I have read through the user manual, and one thing is really confusing me: how do I get my matrix into PETSc?
> 
> My code is organized into separate modules (I am using Fortran).  One of the modules calculates my matrix H and stores it in a compressed sparse row format.  

   What do you mean calculates and stores in compressed sparse row format? Are you calculating it in parallel? Do you want to calculate in parallel? Do you want to solve the eigenvalue problem in parallel? If you just want a sequential code you can use MatCreateSeqAIJWithArrays() to directly use the CSR matrix you created this will be simplier than reorganizing your code to use MatSetValues(). If your code is not  in parallel but you want it to be in parallel then you will need to do a lot of reorganization and use MatSetValues.


   Barry

> I would like to then call from my main program an eigenvalue subroutine and pass the matrix H to it.   In this eigenvalue subroutine, I would use PETSc/SLEPc.  How do I get it to use my pre-calculated matrix?
> 
> For example, here is some code for my eigenvalue subroutine.  In this particular problem, my input matrix is just stored in the regular full matrix format so that I could more easily test the routine.
> 
> subroutine eigen(A)
> #include "finclude/petscdef.h"
> #include "finclude/slepcepsdef.h"
>   use slepceps
>   use kinds_module  !Contains data type definitions, like r8=double precision real.
>   implicit none
>   real (kind=r8), dimension(:,:), intent(in) :: A
> 
> ! local variables
> 
>    PetscErrorCode :: ierr
>    PetscMPIInt :: N            ! Size of A
>    PetscMPIInt :: nconv    ! number of converged eigenpairs
>    Mat :: A_PETSC            ! The PETSc equivalent of my input A
> 
> 
> ! Get size of A
> 
>    N = size(A, dim=1)
> 
> ! Get my matrix into SLEPc.
> 
>    call SlepcInitialize(PETSC_NULL_CHARACTER,ierr); CHKERRQ(ierr)
> 
> < This is where I don't know what to do >
> < How to I do A_PETSC = A>
> 
>    call MatAssemblyBegin(A_PETSC,MAT_FINAL_ASSEMBLY,ierr)
>    call MatAssemblyEnd(A_PETSC,MAT_FINAL_ASSEMBLY,ierr)
> 
> ! Find eigenpairs with SLEPc
> 
>    EPSCreate( PETSC_COMM_WORLD, eps, ierr)
>    EPSSetOperators( eps, A_PETSC, PETSC_NULL )
>    EPSSetProblemType( eps, EPS_HERMITIAN );
>    EPSSetFromOptions( eps );
>    EPSSolve( eps );
>    EPSGetConverged( eps, nconv );
>    ...
>    EPSDestroy( eps );
> 
>    call SlepcFinalize(ierr)
>    end subroutine eigen
> 
> 
> Another thing I am confused about: my main program uses the usual fortran data types.  Should I worry about these clashing with the PETSc ones like PetscMPIInt, Vec, Mat, etc.?
> 
> I apologize if these things have already been discussed at length.  I have been searching the internet and this mailing list, and can't seem to find anything about getting a matrix into PETSc.  I appreciate any help.
> 
> Best regards,
> Alexis
> 
> 
> -- 
> Dr Alexis Morris
> Assistant Professor
> Department of Mathematics, Physics and Engineering
> Mount Royal University
> 4825 Mount Royal Gate SW
> Calgary, Alberta, Canada
> T3E 6K6
> Phone: (403) 440-8507
> Fax: (403) 440-6505
>