[petsc-users] Ordering with MatCreateNest and MPI

[FRANCAVILLA MATTEO ALESSANDRO]

Hi,

I'm trying to use a matrix composed of different sub-matrix blocks; I create 
the matrix with MatCreateNest(). I'm finding it hard to understand the 
ordering (in MPI applications) of the result of MatMult operations. 
Apparently the result of the multiplication is correct, except that both 
input and output vectors go through some sort of permutation. How should I 
use the MatMult routine, i.e. how to pass the input vector with the correct 
ordering and how to get the correct ordering of the result?

Below a simple Fortran example to demonstrate what I'm trying to do (running 
the example with 1 process works fine, with  2 or more processes I don't 
understand how to get the correct ordering). When using standard matrices 
(i.e., not MatNest) I've never had to care about orderings.

Thanks,
Alessandro



PROGRAM Main

!.....================================

implicit none

#include <petsc/finclude/petscsys.h>
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>

!*********************************************************************
!                       PETSc variables
!*********************************************************************

Vec						:: x, b
Mat						:: matArray(4), A
PetscInt				:: II, JJ, FromRow, ToRow, IMAT
PetscErrorCode    		:: IERR
PetscInt, ALLOCATABLE     :: INDICES(:)
PetscScalar, ALLOCATABLE  :: VALUES(:)
PetscInt, PARAMETER		:: N = 2			! size of sub-matrices

CALL PetscInitialize(PETSC_NULL_CHARACTER,IERR)


ALLOCATE(INDICES(N),VALUES(N))
INDICES = (/ (ii, ii=0, N-1) /)

! Form 4 blocks of size NxN

DO IMAT = 1, 4
	CALL 
MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,PETSC_NULL_SCALAR,matArray(IMAT),IERR)
	CALL MatGetOwnershipRange(matArray(IMAT),FromRow,ToRow,IERR)
	DO II = FromRow, ToRow-1
		DO JJ = 1, N
	    	VALUES(JJ) = (IMAT-1)*100 + N*II + JJ
		END DO
	  	CALL 
MatSetValues(matArray(IMAT),1,II,N,INDICES,VALUES,INSERT_VALUES,IERR)
	END DO
	CALL MatAssemblyBegin(matArray(IMAT),MAT_FINAL_ASSEMBLY,IERR)
	CALL MatAssemblyEnd(matArray(IMAT),MAT_FINAL_ASSEMBLY,IERR)
END DO

! Form a 2x2 block matrix made as [ matArray(1) matArray(2) ; matArray(3) 
matArray(4) ]
! If N = 2 the matrix is:
! A =            [   1   2 | 101 102 ;
!          	  3   4 | 102 103 ;
!		___________________
!               201 202 | 401 402 ;
!	     203 204 | 403 404 ]

CALL 
MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matArray,A,IERR)

CALL MatCreateVecs(A,x,b,IERR)
CALL VecSet(x,0.+0.*PETSC_i,IERR)
CALL VecSetValues(x,1,2,(1.+0.*PETSC_i),INSERT_VALUES,IERR)
CALL VecAssemblyBegin(x,IERR)
CALL VecAssemblyEnd(x,IERR)
CALL VecView(x,PETSC_VIEWER_STDOUT_WORLD,IERR)

CALL MatMult(A,x,b,IERR)
! A*x extracts the 3rd column of A

CALL VecView(b,PETSC_VIEWER_STDOUT_WORLD,IERR)

DEALLOCATE(INDICES,VALUES)

CALL PetscFinalize(IERR)

END PROGRAM