[petsc-users] Getting a matrix into PETSc

[Matthew Knepley]

On Mon, Sep 27, 2010 at 6:10 PM, Alexis Morris <amorris at mtroyal.ca> 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.  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?
>

The best thing to do is replace your data structure with the PETSc Mat
object. I suppose that you construct your
CSR matrix one row at a time. You can do this with PETSc using


http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatSetValues.html

so it should be simple translation. To be efficient, there is the one
additional step of telling PETSc the
length of each row, but you must do this to assemble your CSR structure, so
again it should be no extra work.

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.?
>

When calling PETSc functions, it is best to use the PETSc types, such as
PetscInt.

  Thanks,

     Matt


> 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
>
>


-- 
What most experimenters take for granted before they begin their experiments
is infinitely more interesting than any results to which their experiments
lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20100927/0d261dd0/attachment-0001.htm>