Hello everyone,
I have a question on setting matrix entries which are not in the stencil
width. Take ksp ex45.c as an example,
http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex45.c.html
Instead of setting up the standard 7-point stencil, now for each cell, the
matrix also has a additional dependency on the cell (Mx, My, Mz). Namely,
for each row, the col corresponding to (Mx, My, Mz) is always nonzero. I
modify the example code to add this entries:
+ MatSetOption(B,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
+ MatSetOption(jac,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
+ v[7] = 100; col[7].i = mx-1; col[7].j = my-1; col[7].k = mz-1;
+ ierr = MatSetValuesStencil(B,1,&row,8,col,v,INSERT_VALUES);CHKERRQ(ierr);
It is okay to for np=1, but crash for np>=2 with the error message:
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Local index 342 too large 244 (max) at 7
[0]PETSC ERROR: #1 ISLocalToGlobalMappingApply() line 423 in
petsc-3.7.4/src/vec/is/utils/isltog.c
[0]PETSC ERROR: #2 MatSetValuesLocal() line 2052 in
petsc-3.7.4/src/mat/interface/matrix.c
[0]PETSC ERROR: #3 MatSetValuesStencil() line 1447 in
petsc-3.7.4/src/mat/interface/matrix.c
[0]PETSC ERROR: #4 ComputeMatrix() line 151 in extest.c
Can I add new entries to the cols not in the stencil width into the dmda
matrix or Jacobian?
Attached please find the modifed ex45 example, the diff file as well as the
run log.
Thanks for your help.
Best,
Xiangdong
diff --git a/ex45.c b/ex45.c
index b4bb565..6bcb074 100644
--- a/ex45.c
+++ b/ex45.c
@@ -117,8 +117,8 @@ PetscErrorCode ComputeMatrix(KSP ksp,Mat jac,Mat B,void
*ctx)
DM da;
PetscErrorCode ierr;
PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
- PetscScalar v[7],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
- MatStencil row,col[7];
+ PetscScalar v[8],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
+ MatStencil row,col[8];
PetscFunctionBeginUser;
ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
@@ -127,6 +127,9 @@ PetscErrorCode ComputeMatrix(KSP ksp,Mat jac,Mat B,void
*ctx)
HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
+ MatSetOption(B,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
+ MatSetOption(jac,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
+
for (k=zs; k<zs+zm; k++) {
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
@@ -142,7 +145,10 @@ PetscErrorCode ComputeMatrix(KSP ksp,Mat jac,Mat B,void
*ctx)
v[4] = -HyHzdHx;col[4].i = i+1; col[4].j = j; col[4].k = k;
v[5] = -HxHzdHy;col[5].i = i; col[5].j = j+1; col[5].k = k;
v[6] = -HxHydHz;col[6].i = i; col[6].j = j; col[6].k = k+1;
- ierr =
MatSetValuesStencil(B,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr);
+ //ierr =
MatSetValuesStencil(B,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr);
+
+ v[7] = 100; col[7].i = mx-1; col[7].j = my-1; col[7].k = mz-1;
+ ierr =
MatSetValuesStencil(B,1,&row,8,col,v,INSERT_VALUES);CHKERRQ(ierr);
}
}
}
/*
Laplacian in 3D. Modeled by the partial differential equation
- Laplacian u = 1,0 < x,y,z < 1,
with boundary conditions
u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
This uses multigrid to solve the linear system
See src/snes/examples/tutorials/ex50.c
Can also be run with -pc_type exotic -ksp_type fgmres
*/
static char help[] = "Solves 3D Laplacian using multigrid.\n\n";
#include <petscksp.h>
#include <petscdm.h>
#include <petscdmda.h>
extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
extern PetscErrorCode ComputeInitialGuess(KSP,Vec,void*);
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
PetscErrorCode ierr;
KSP ksp;
PetscReal norm;
DM da;
Vec x,b,r;
Mat A;
PetscInitialize(&argc,&argv,(char*)0,help);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-7,-7,-7,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);CHKERRQ(ierr);
ierr = KSPSetDM(ksp,da);CHKERRQ(ierr);
ierr = KSPSetComputeInitialGuess(ksp,ComputeInitialGuess,NULL);CHKERRQ(ierr);
ierr = KSPSetComputeRHS(ksp,ComputeRHS,NULL);CHKERRQ(ierr);
ierr = KSPSetComputeOperators(ksp,ComputeMatrix,NULL);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
ierr = KSPGetRhs(ksp,&b);CHKERRQ(ierr);
ierr = VecDuplicate(b,&r);CHKERRQ(ierr);
ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr);
ierr = MatMult(A,x,r);CHKERRQ(ierr);
ierr = VecAXPY(r,-1.0,b);CHKERRQ(ierr);
ierr = VecNorm(r,NORM_2,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);CHKERRQ(ierr);
ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
#undef __FUNCT__
#define __FUNCT__ "ComputeRHS"
PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
{
PetscErrorCode ierr;
PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
DM dm;
PetscScalar Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
PetscScalar ***barray;
PetscFunctionBeginUser;
ierr = KSPGetDM(ksp,&dm);CHKERRQ(ierr);
ierr = DMDAGetInfo(dm,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
ierr = DMDAVecGetArray(dm,b,&barray);CHKERRQ(ierr);
for (k=zs; k<zs+zm; k++) {
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
barray[k][j][i] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
} else {
barray[k][j][i] = Hx*Hy*Hz;
}
}
}
}
ierr = DMDAVecRestoreArray(dm,b,&barray);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "ComputeInitialGuess"
PetscErrorCode ComputeInitialGuess(KSP ksp,Vec b,void *ctx)
{
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = VecSet(b,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "ComputeMatrix"
PetscErrorCode ComputeMatrix(KSP ksp,Mat jac,Mat B,void *ctx)
{
DM da;
PetscErrorCode ierr;
PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
PetscScalar v[8],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
MatStencil row,col[8];
PetscFunctionBeginUser;
ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
ierr = DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
MatSetOption(B,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
MatSetOption(jac,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);
for (k=zs; k<zs+zm; k++) {
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
row.i = i; row.j = j; row.k = k;
if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
v[0] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
ierr = MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
} else {
v[0] = -HxHydHz;col[0].i = i; col[0].j = j; col[0].k = k-1;
v[1] = -HxHzdHy;col[1].i = i; col[1].j = j-1; col[1].k = k;
v[2] = -HyHzdHx;col[2].i = i-1; col[2].j = j; col[2].k = k;
v[3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);col[3].i = row.i; col[3].j = row.j; col[3].k = row.k;
v[4] = -HyHzdHx;col[4].i = i+1; col[4].j = j; col[4].k = k;
v[5] = -HxHzdHy;col[5].i = i; col[5].j = j+1; col[5].k = k;
v[6] = -HxHydHz;col[6].i = i; col[6].j = j; col[6].k = k+1;
//ierr = MatSetValuesStencil(B,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr);
v[7] = 100; col[7].i = mx-1; col[7].j = my-1; col[7].k = mz-1;
ierr = MatSetValuesStencil(B,1,&row,8,col,v,INSERT_VALUES);CHKERRQ(ierr);
}
}
}
}
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
DM Object: 2 MPI processes
type: da
Processor [0] M 7 N 7 P 7 m 1 n 1 p 2 w 1 s 1
X range of indices: 0 7, Y range of indices: 0 7, Z range of indices: 0 4
Processor [1] M 7 N 7 P 7 m 1 n 1 p 2 w 1 s 1
X range of indices: 0 7, Y range of indices: 0 7, Z range of indices: 4 7
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Local index 342 too large 244 (max) at 7
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.7.4, Oct, 02, 2016
[0]PETSC ERROR: #1 ISLocalToGlobalMappingApply() line 423 in
petsc-3.7.4/src/vec/is/utils/isltog.c
[0]PETSC ERROR: #2 MatSetValuesLocal() line 2052 in
petsc-3.7.4/src/mat/interface/matrix.c
[0]PETSC ERROR: #3 MatSetValuesStencil() line 1447 in
petsc-3.7.4/src/mat/interface/matrix.c
[0]PETSC ERROR: #4 ComputeMatrix() line 151 in extest.c
[0]PETSC ERROR: #5 KSPSetUp() line 341 in
petsc-3.7.4/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #6 KSPSolve() line 599 in
petsc-3.7.4/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #7 main() line 51 in extest.c
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -dm_view
[0]PETSC ERROR: ----------------End of Error Message -------send entire error
message to [email protected]
