<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Jul 28, 2016 at 3:38 AM, Klaij, Christiaan <span dir="ltr"><<a href="mailto:C.Klaij@marin.nl" target="_blank">C.Klaij@marin.nl</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex">I'm trying to understand how to assemble a block matrix in a<br>
format-independent manner, so that I can switch between types<br>
mpiaij and matnest.<br>
<br>
The manual states that the key to format-independent assembly is<br>
to use MatGetLocalSubMatrix. So, in the code below, I'm using<br>
this to assemble a 3-by-3 block matrix A and setting the diagonal<br>
of block A02. This seems to work for type mpiaij, but not for<br>
type matnest. What am I missing?<br>
<br>
Chris<br>
<br>
<br>
$ cat mattry.F90<br>
program mattry<br>
<br>
use petscksp<br>
implicit none<br>
#include <petsc/finclude/petsckspdef.h><br>
<br>
PetscInt :: n=4 ! setting 4 cells per process<br>
<br>
PetscErrorCode :: ierr<br>
PetscInt :: size,rank,i<br>
Mat :: A,A02<br>
IS :: isg0,isg1,isg2<br>
IS :: isl0,isl1,isl2<br>
ISLocalToGlobalMapping :: map<br>
<br>
integer, allocatable, dimension(:) :: idx<br>
<br>
call PetscInitialize(PETSC_NULL_CHARACTER,ierr); CHKERRQ(ierr)<br>
call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr); CHKERRQ(ierr)<br>
call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRQ(ierr)<br>
<br>
! local index sets for 3 fields<br>
allocate(idx(n))<br>
idx=(/ (i-1, i=1,n) /)<br>
call ISCreateGeneral(PETSC_COMM_WORLD,n,idx,PETSC_COPY_VALUES,isl0,ierr);CHKERRQ(ierr)<br>
call ISCreateGeneral(PETSC_COMM_WORLD,n,idx+n,PETSC_COPY_VALUES,isl1,ierr);CHKERRQ(ierr)<br>
call ISCreateGeneral(PETSC_COMM_WORLD,n,idx+2*n,PETSC_COPY_VALUES,isl2,ierr);CHKERRQ(ierr)<br>
! call ISView(isl3,PETSC_VIEWER_STDOUT_WORLD,ierr); CHKERRQ(ierr)<br>
deallocate(idx)<br>
<br>
! global index sets for 3 fields<br>
allocate(idx(n))<br>
idx=(/ (i-1+rank*3*n, i=1,n) /)<br>
call ISCreateGeneral(PETSC_COMM_WORLD,n,idx,PETSC_COPY_VALUES,isg0,ierr);CHKERRQ(ierr)<br>
call ISCreateGeneral(PETSC_COMM_WORLD,n,idx+n,PETSC_COPY_VALUES,isg1,ierr); CHKERRQ(ierr)<br>
call ISCreateGeneral(PETSC_COMM_WORLD,n,idx+2*n,PETSC_COPY_VALUES,isg2,ierr); CHKERRQ(ierr)<br>
! call ISView(isg3,PETSC_VIEWER_STDOUT_WORLD,ierr); CHKERRQ(ierr)<br>
deallocate(idx)<br>
<br>
! local-to-global mapping<br>
allocate(idx(3*n))<br>
idx=(/ (i-1+rank*3*n, i=1,3*n) /)<br>
call ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,3*n,idx,PETSC_COPY_VALUES,map,ierr); CHKERRQ(ierr)<br>
! call ISLocalToGlobalMappingView(map,PETSC_VIEWER_STDOUT_WORLD,ierr); CHKERRQ(ierr)<br>
deallocate(idx)<br>
<br>
! create the 3-by-3 block matrix<br>
call MatCreate(PETSC_COMM_WORLD,A,ierr); CHKERRQ(ierr)<br>
call MatSetSizes(A,3*n,3*n,PETSC_DECIDE,PETSC_DECIDE,ierr); CHKERRQ(ierr)<br>
! call MatSetType(A,MATNEST,ierr); CHKERRQ(ierr)<br>
call MatSetUp(A,ierr); CHKERRQ(ierr)<br></blockquote><div><br></div><div>I am sorry I have not had time to run this, but I believe you need to insert a call here:</div><div><pre style="color:rgb(0,0,0)"><span style="font-family:arial,sans-serif"> MatNestSetSubMats(A, 3, [isg0, isg1, isg2], 3, [isg0, isg1, isg2], NULL);</span><br></pre></div><div>coming from <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatNestSetSubMats.html#MatNestSetSubMats">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatNestSetSubMats.html#MatNestSetSubMats</a></div><div><br></div><div>Thanks,</div><div><br></div><div> Matt</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex">
call MatSetOptionsPrefix(A,"A_",ierr); CHKERRQ(ierr)<br>
call MatSetLocalToGlobalMapping(A,map,map,ierr); CHKERRQ(ierr)<br>
call MatSetFromOptions(A,ierr); CHKERRQ(ierr)<br>
<br>
! set diagonal of block A02 to 0.65<br>
call MatGetLocalSubmatrix(A,isl0,isl2,A02,ierr); CHKERRQ(ierr)<br>
do i=1,n<br>
call MatSetValuesLocal(A02,1,i-1,1,i-1,0.65d0,INSERT_VALUES,ierr); CHKERRQ(ierr)<br>
end do<br>
call MatRestoreLocalSubMatrix(A,isl0,isl2,A02,ierr); CHKERRQ(ierr)<br>
call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)<br>
call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)<br>
<br>
! verify<br>
call MatGetSubmatrix(A,isg0,isg2,MAT_INITIAL_MATRIX,A02,ierr); CHKERRQ(ierr)<br>
call MatView(A02,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr)<br>
<br>
call PetscFinalize(ierr)<br>
<br>
end program mattry<br>
<br>
$ mpiexec -n 2 ./mattry -A_mat_type mpiaij<br>
Mat Object: 2 MPI processes<br>
type: mpiaij<br>
row 0: (0, 0.65)<br>
row 1: (1, 0.65)<br>
row 2: (2, 0.65)<br>
row 3: (3, 0.65)<br>
row 4: (4, 0.65)<br>
row 5: (5, 0.65)<br>
row 6: (6, 0.65)<br>
row 7: (7, 0.65)<br>
<br>
$ mpiexec -n 2 ./mattry -A_mat_type nest<br>
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>
[0]PETSC ERROR: Null argument, when expecting valid pointer<br>
[0]PETSC ERROR: Null Pointer: Parameter # 3<br>
[0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>
[0]PETSC ERROR: Petsc Release Version 3.7.3, Jul, 24, 2016<br>
[0]PETSC ERROR: ./mattry on a linux_64bit_debug named lin0322.marin.local by cklaij Thu Jul 28 10:31:04 2016<br>
[0]PETSC ERROR: Configure options --with-mpi-dir=/home/cklaij/ReFRESCO/Dev/trunk/Libs/install/openmpi/1.8.7 --with-clanguage=c++ --with-x=1 --with-debugging=1 --with-blas-lapack-dir=/opt/intel/composer_xe_2015.1.133/mkl --with-shared-libraries=0<br>
[0]PETSC ERROR: #1 MatNestFindIS() line 298 in /home/cklaij/ReFRESCO/Dev/trunk/Libs/build/petsc/3.7.3-dbg/src/mat/impls/nest/matnest.c<br>
[0]PETSC ERROR: [1]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>
[1]PETSC ERROR: Null argument, when expecting valid pointer<br>
[1]PETSC ERROR: Null Pointer: Parameter # 3<br>
[1]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>
[1]PETSC ERROR: Petsc Release Version 3.7.3, Jul, 24, 2016<br>
[1]PETSC ERROR: ./mattry on a linux_64bit_debug named lin0322.marin.local by cklaij Thu Jul 28 10:31:04 2016<br>
[1]PETSC ERROR: Configure options --with-mpi-dir=/home/cklaij/ReFRESCO/Dev/trunk/Libs/install/openmpi/1.8.7 --with-clanguage=c++ --with-x=1 --with-debugging=1 --with-blas-lapack-dir=/opt/intel/composer_xe_2015.1.133/mkl --with-shared-libraries=0<br>
[1]PETSC ERROR: #1 MatNestFindIS() line 298 in /home/cklaij/ReFRESCO/Dev/trunk/Libs/build/petsc/3.7.3-dbg/src/mat/impls/nest/matnest.c<br>
[1]PETSC ERROR: #2 MatNestFindSubMat() line 371 in /home/cklaij/ReFRESCO/Dev/trunk/Libs/build/petsc/3.7.3-dbg/src/mat/impls/nest/matnest.c<br>
[1]PETSC ERROR: #3 MatGetLocalSubMatrix_Nest() line 414 in /home/cklaij/ReFRESCO/Dev/trunk/Libs/build/petsc/3.7.3-dbg/src/mat/impls/nest/matnest.c<br>
[1]PETSC ERROR: #4 MatGetLocalSubMatrix() line 10099 in /home/cklaij/ReFRESCO/Dev/trunk/Libs/build/petsc/3.7.3-dbg/src/mat/interface/matrix.c<br>
#2 MatNestFindSubMat() line 371 in /home/cklaij/ReFRESCO/Dev/trunk/Libs/build/petsc/3.7.3-dbg/src/mat/impls/nest/matnest.c<br>
[0]PETSC ERROR: #3 MatGetLocalSubMatrix_Nest() line 414 in /home/cklaij/ReFRESCO/Dev/trunk/Libs/build/petsc/3.7.3-dbg/src/mat/impls/nest/matnest.c<br>
[0]PETSC ERROR: #4 MatGetLocalSubMatrix() line 10099 in /home/cklaij/ReFRESCO/Dev/trunk/Libs/build/petsc/3.7.3-dbg/src/mat/interface/matrix.c<br>
--------------------------------------------------------------------------<br>
MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD<br>
with errorcode 85.<br>
<br>
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.<br>
You may or may not see output from other processes, depending on<br>
exactly when Open MPI kills them.<br>
--------------------------------------------------------------------------<br>
[lin0322.marin.local:11985] 1 more process has sent help message help-mpi-api.txt / mpi-abort<br>
[lin0322.marin.local:11985] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages<br>
$<br>
<br>
<br>
dr. ir. Christiaan Klaij | CFD Researcher | Research & Development<br>
MARIN | T +31 317 49 33 44 | mailto:<a href="mailto:C.Klaij@marin.nl">C.Klaij@marin.nl</a> | <a href="http://www.marin.nl" rel="noreferrer" target="_blank">http://www.marin.nl</a><br>
<br>
MARIN news: <a href="http://www.marin.nl/web/News/News-items/Ship-design-in-EU-project-Holiship.htm" rel="noreferrer" target="_blank">http://www.marin.nl/web/News/News-items/Ship-design-in-EU-project-Holiship.htm</a><br>
<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">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</div>
</div></div>