On Mon, Sep 27, 2010 at 6:10 PM, Alexis Morris <span dir="ltr"><<a href="mailto:amorris@mtroyal.ca">amorris@mtroyal.ca</a>></span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Hi everyone,<br>
<br>
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?<br>
<br>
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?<br>
</blockquote><div><br></div><div>The best thing to do is replace your data structure with the PETSc Mat object. I suppose that you construct your</div><div>CSR matrix one row at a time. You can do this with PETSc using</div>
<div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatSetValues.html">http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatSetValues.html</a></div>
<div> </div><div>so it should be simple translation. To be efficient, there is the one additional step of telling PETSc the</div><div>length of each row, but you must do this to assemble your CSR structure, so again it should be no extra work.</div>
<div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
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.<br>
<br>
subroutine eigen(A)<br>
#include "finclude/petscdef.h"<br>
#include "finclude/slepcepsdef.h"<br>
use slepceps<br>
use kinds_module !Contains data type definitions, like r8=double precision real.<br>
implicit none<br>
real (kind=r8), dimension(:,:), intent(in) :: A<br>
<br>
! local variables<br>
<br>
PetscErrorCode :: ierr<br>
PetscMPIInt :: N ! Size of A<br>
PetscMPIInt :: nconv ! number of converged eigenpairs<br>
Mat :: A_PETSC ! The PETSc equivalent of my input A<br>
<br>
<br>
! Get size of A<br>
<br>
N = size(A, dim=1)<br>
<br>
! Get my matrix into SLEPc.<br>
<br>
call SlepcInitialize(PETSC_NULL_CHARACTER,ierr); CHKERRQ(ierr)<br>
<br>
< This is where I don't know what to do ><br>
< How to I do A_PETSC = A><br>
<br>
call MatAssemblyBegin(A_PETSC,MAT_FINAL_ASSEMBLY,ierr)<br>
call MatAssemblyEnd(A_PETSC,MAT_FINAL_ASSEMBLY,ierr)<br>
<br>
! Find eigenpairs with SLEPc<br>
<br>
EPSCreate( PETSC_COMM_WORLD, eps, ierr)<br>
EPSSetOperators( eps, A_PETSC, PETSC_NULL )<br>
EPSSetProblemType( eps, EPS_HERMITIAN );<br>
EPSSetFromOptions( eps );<br>
EPSSolve( eps );<br>
EPSGetConverged( eps, nconv );<br>
...<br>
EPSDestroy( eps );<br>
<br>
call SlepcFinalize(ierr)<br>
end subroutine eigen<br>
<br>
<br>
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.?<br></blockquote><div><br></div><div>When calling PETSc functions, it is best to use the PETSc types, such as PetscInt.</div>
<div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
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.<br>
<br>
Best regards,<br>
Alexis<br>
<br>
<br>
-- <br>
Dr Alexis Morris<br>
Assistant Professor<br>
Department of Mathematics, Physics and Engineering<br>
Mount Royal University<br>
4825 Mount Royal Gate SW<br>
Calgary, Alberta, Canada<br>
T3E 6K6<br>
Phone: (403) 440-8507<br>
Fax: (403) 440-6505<br>
<br>
</blockquote></div><br><br clear="all"><br>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener<br>