I am using Petsc Ksp routines. I construct an operator using MatSetValuesStencil, where in each call of this function I specify one row matrix values of length 5. There is a case where I sometimes need to completely replace a row from a 5 length stencil to a 3 length one. Will INSERT_VALUES mode leave the two values on non changed positions or it will discard them to zero?
The elements of the matrix that are not specified in the arguments idxm and idxn of the function MatSetValuesStencil(...) are left unchanged, even if INSERT_VALUES is used.
Here is a little code starting from ksp_ex29 to test it :
static char help[] = "Does INSERT_VALUES changes the whole row ? No.\n\n";
#include <petscdm.h>
#include <petscdmda.h>
#include <petscksp.h>
extern PetscErrorCode ComputeMatrix42(DM da,Mat jac);
extern PetscErrorCode ComputeMatrix(DM da,Mat jac);
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
DM da;
PetscErrorCode ierr;
Mat matrix;
PetscInitialize(&argc,&argv,(char*)0,help);
ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);CHKERRQ(ierr);
DMCreateMatrix(da,&matrix);
ComputeMatrix(da,matrix);
PetscPrintf(PETSC_COMM_WORLD,"A matrix of negative terms : \n");
MatView(matrix, PETSC_VIEWER_STDOUT_WORLD );
ComputeMatrix42(da,matrix);
PetscPrintf(PETSC_COMM_WORLD,"The diagonal, i-1 and i+1 are set to 42 : \n");
MatView(matrix, PETSC_VIEWER_STDOUT_WORLD );
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = MatDestroy(&matrix);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
#undef __FUNCT__
#define __FUNCT__ "ComputeMatrix"
PetscErrorCode ComputeMatrix(DM da,Mat jac)
{
PetscErrorCode ierr;
PetscInt i,j,mx,my,xm,ym,xs,ys;
PetscScalar v[5];
MatStencil row, col[5];
PetscFunctionBeginUser;
ierr = DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
row.i = i; row.j = j;
v[0] = -1; col[0].i = i; col[0].j = j-1;
v[1] = -1; col[1].i = i-1; col[1].j = j;
v[2] = -13; col[2].i = i; col[2].j = j;
v[3] = -1; col[3].i = i+1; col[3].j = j;
v[4] = -1; col[4].i = i; col[4].j = j+1;
ierr = MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "ComputeMatrix42"
PetscErrorCode ComputeMatrix42(DM da,Mat jac)
{
PetscErrorCode ierr;
PetscInt i,j,mx,my,xm,ym,xs,ys;
PetscScalar v[3];
MatStencil row, col[3];
PetscFunctionBeginUser;
ierr = DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
row.i = i; row.j = j;
v[0] = 42; col[0].i = i-1; col[0].j = j;
v[1] = 42; col[1].i = i; col[1].j = j;
v[2] = 42; col[2].i = i+1; col[2].j = j;
ierr = MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
Compile it with the following makefile :
include ${PETSC_DIR}/conf/variables
include ${PETSC_DIR}/conf/rules
main: main.o chkopts
-${CLINKER} -o main main.o ${PETSC_LIB}
${RM} main.o
Output :
A matrix of negative terms :
Mat Object: 1 MPI processes
type: seqaij
row 0: (0, -13) (1, -1) (2, -1) (3, -1) (6, -1)
row 1: (0, -1) (1, -13) (2, -1) (4, -1) (7, -1)
row 2: (0, -1) (1, -1) (2, -13) (5, -1) (8, -1)
row 3: (0, -1) (3, -13) (4, -1) (5, -1) (6, -1)
row 4: (1, -1) (3, -1) (4, -13) (5, -1) (7, -1)
row 5: (2, -1) (3, -1) (4, -1) (5, -13) (8, -1)
row 6: (0, -1) (3, -1) (6, -13) (7, -1) (8, -1)
row 7: (1, -1) (4, -1) (6, -1) (7, -13) (8, -1)
row 8: (2, -1) (5, -1) (6, -1) (7, -1) (8, -13)
The diagonal, i-1 and i+1 are set to 42 :
Mat Object: 1 MPI processes
type: seqaij
row 0: (0, 42) (1, 42) (2, 42) (3, -1) (6, -1)
row 1: (0, 42) (1, 42) (2, 42) (4, -1) (7, -1)
row 2: (0, 42) (1, 42) (2, 42) (5, -1) (8, -1)
row 3: (0, -1) (3, 42) (4, 42) (5, 42) (6, -1)
row 4: (1, -1) (3, 42) (4, 42) (5, 42) (7, -1)
row 5: (2, -1) (3, 42) (4, 42) (5, 42) (8, -1)
row 6: (0, -1) (3, -1) (6, 42) (7, 42) (8, 42)
row 7: (1, -1) (4, -1) (6, 42) (7, 42) (8, 42)
row 8: (2, -1) (5, -1) (6, 42) (7, 42) (8, 42)
I am using PETSC 3.5.2.