Hello everyone,
I want to use PCApplyBAorAB symmetrically with a preconditioner from an ICC
factorization.
However, I don't understand why the result is not the same if I use a
natural or a RCM ordering during the factorization.
This problem does not occur if I use this routine with PC_LEFT or PC_RIGHT
I have put in attachment a code to illustrate the problem.
I compute a symmetric matix A in R^{n*n} and a full rank matrix B in R^{m*n}
I compute an ICC factorization of the SPD matrix A+B^TB to precondition A.
Thank you!
Sylvain Mercier
program main
implicit none
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! include files
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! petscsys.h - base PETSc routines petscvec.h - vecteurs
! petscmat.h - matrices petscpc.h - preconditionneurs
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscmat.h>
#include <finclude/petscpc.h>
Vec x,y,work
Mat A,K,C,Ct
PC pc
PetscScalar one,none
PetscScalar value(3)
PetscInt i,i1,i2,i3,intbid
PetscInt m,n,col(3)
PetscMPIInt rank
PetscBool flg
PetscErrorCode ierr
MatFactorInfo info
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
none = -1.D0
one = 1.D0
n = 10
m = 5
i1 = 1
i2 = 2
i3 = 3
call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
! -------------------------------------------------------------------
! Compute a symmetric matrix A (size n,n) and
! a full rank matrix C (size m,n)
! -------------------------------------------------------------------
call MatCreate(PETSC_COMM_WORLD,A,ierr)
call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
call MatSetFromOptions(A,ierr)
value(1) = -1.0
value(2) = 2.0
value(3) = -1.0
do 60 i=1,n-2
col(1) = i-1
col(2) = i
col(3) = i+1
call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr)
60 continue
i = n - 1
col(1) = n - 2
col(2) = n - 1
call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
i = 0
col(1) = 0
col(2) = 1
value(1) = 2.0
value(2) = -1.0
call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
call MatCreate(PETSC_COMM_WORLD,C,ierr)
call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m,n,ierr)
call MatSetFromOptions(C,ierr)
call MatSetValue(C,0,0,one,INSERT_VALUES,ierr)
call MatSetValue(C,0,5,one,INSERT_VALUES,ierr)
call MatSetValue(C,1,1,one,INSERT_VALUES,ierr)
call MatSetValue(C,2,2,one,INSERT_VALUES,ierr)
call MatSetValue(C,3,3,one,INSERT_VALUES,ierr)
call MatSetValue(C,3,4,none,INSERT_VALUES,ierr)
call MatSetValue(C,4,5,one,INSERT_VALUES,ierr)
call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
! call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)
! call MatView(C,PETSC_VIEWER_STDOUT_WORLD,ierr)
! -------------------------------------------------------------------
! Compute K=A+C^TC (symmetric positive definite)
! -------------------------------------------------------------------
call MatTranspose(C,MAT_INITIAL_MATRIX,Ct,ierr)
call MatMatMult(Ct,C,MAT_INITIAL_MATRIX,1.D0,K,ierr)
call MatAXPY(K,1.D0,A,DIFFERENT_NONZERO_PATTERN,ierr)
! call MatView(K,PETSC_VIEWER_STDOUT_WORLD,ierr)
! -------------------------------------------------------------------
! Create the PC object to precondition A with an ICC of K
! and compute an example of PCApplyBAorAB (symmetric)
! -------------------------------------------------------------------
call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n,x,ierr)
call VecDuplicate(x,y,ierr)
call VecDuplicate(y,work,ierr)
call VecSet(x,one,ierr)
! Using MATORDERINGNATURAL
call PCCreate(PETSC_COMM_WORLD,pc,ierr)
call PCSetOperators(pc,A,K,SAME_NONZERO_PATTERN,ierr)
call PCSetType(pc,PCICC,ierr)
call PCFactorSetLevels(pc,5,ierr)
call PCFactorSetMatOrderingType(pc,MATORDERINGNATURAL,ierr)
call PCApplyBAorAB(pc,PC_SYMMETRIC,x,y,work,ierr)
write(*,*)'Natural Ordering:'
call VecView(y,PETSC_VIEWER_STDOUT_WORLD,ierr)
! Using MATORDERING RCM
call PCCreate(PETSC_COMM_WORLD,pc,ierr)
call PCSetOperators(pc,A,K,SAME_NONZERO_PATTERN,ierr)
call PCSetType(pc,PCICC,ierr)
call PCFactorSetLevels(pc,5,ierr)
call PCFactorSetMatOrderingType(pc,MATORDERINGRCM,ierr)
call PCApplyBAorAB(pc,PC_SYMMETRIC,x,y,work,ierr)
write(*,*)'RCM Ordering:'
call VecView(y,PETSC_VIEWER_STDOUT_WORLD,ierr)
! Free work space. All PETSc objects should be destroyed when they
! are no longer needed.
call VecDestroy(x,ierr)
call VecDestroy(y,ierr)
call VecDestroy(work,ierr)
call MatDestroy(A,ierr)
call MatDestroy(C,ierr)
call MatDestroy(Ct,ierr)
call MatDestroy(K,ierr)
! Always call PetscFinalize() before exiting a program.
call PetscFinalize(ierr)
end
