<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Thu, Jul 4, 2013 at 10:34 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
The j should be tracking the y variable and the i the x variable. See for example src/ksp/ksp/examples/tutorials/ex29.c<br>
<br>
for (j=ys; j<ys+ym; j++) {<br>
for (i=xs; i<xs+xm; i++) {<br>
<div class="im"> row.i = i; row.j = j;<br>
<br></div></blockquote><div>Thanks Barry, tracking y variable with i, and x with j was the problem! Interchanging them worked. <br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div class="im">
</div> After you fix this run with a 3 by 4 mesh and check the line number printed indicating the problem. Then look at the row and column values for that problem and see in the code how it could get out of bounds.<br>
</blockquote><div><br></div><div>I'm afraid I did not understand exactly which line number that gets printed are you talking about ? Is it this one in the error message ?<br> [0]PETSC ERROR: main() line 75 in "unknowndirectory/"/user/bkhanal/home/works/cmake_tuts/petsc_test/examples/FDdmda_test.cxx <br>
</div><div>where there is a call to the function MatSetValuesStencil ? <br></div><div>Could you please clarify a bit more ?<br><br></div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<span class=""><font color="#888888"><br>
Barry<br>
</font></span><div class=""><div class="h5"><br>
On Jul 4, 2013, at 8:11 AM, Bishesh Khanal <<a href="mailto:bisheshkh@gmail.com">bisheshkh@gmail.com</a>> wrote:<br>
<br>
><br>
> On Thu, Jul 4, 2013 at 2:59 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
> On Thu, Jul 4, 2013 at 7:32 AM, Bishesh Khanal <<a href="mailto:bisheshkh@gmail.com">bisheshkh@gmail.com</a>> wrote:<br>
> Hi all,<br>
> I'm trying to use DMCreateMatrix and MatStencil to fill up a matrix that results from a finite difference discretization of a PDE (objective is to solve the resulting linear system in the form of Ax=b).<br>
><br>
> However there are some errors and I'm not sure about the issue!<br>
> A short demo:<br>
> A 2D mXn grid, with two variables at each node (dof=2),<br>
> so the resulting A would be 2mn X 2mn.<br>
> Let's say the two variables are vx and vy, and the two associated<br>
> equations discretized are x-eq and y-eq.<br>
><br>
> Here is the relevant part of the code: (mat1 variable for the A matrix):<br>
><br>
> Just sending the whole code is better. I suspect the logic is wrong for selecting the boundary.<br>
><br>
> Matt<br>
><br>
> Thanks, here is the complete code:<br>
><br>
> #include <petscsys.h><br>
> #include <petscdmda.h><br>
><br>
> #undef __FUNCT__<br>
> #define __FUNCT__ "main"<br>
> int main(int argc,char **argv)<br>
> {<br>
> PetscErrorCode ierr;<br>
> DM da;<br>
> PetscInt m=10,n=10;<br>
><br>
> ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);<br>
><br>
> /* Create a DMDA and an associated vector */<br>
> ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,<br>
> PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr);<br>
><br>
> Mat mat1;<br>
> ierr = DMCreateMatrix(da,MATMPIAIJ,&mat1); CHKERRQ(ierr);<br>
> MatStencil row, col[4]; //let's test 4 non-zeros in one row.<br>
> PetscScalar val[4];<br>
> PetscScalar coeff = 2.; //just a constant coefficient for testing.<br>
> PetscInt i,j;<br>
> DMDALocalInfo info;<br>
><br>
> ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);<br>
><br>
> for(i = info.ys; i < info.ys+info.ym; ++i){<br>
> for(j = info.xs; j < info.xs+info.xm; ++j){<br>
> row.i = i; row.j = j; //one node at a time<br>
> if (j == 0 || j == info.mx-1){ //left and right borders<br>
> //vx(i,j) = 0;<br>
> row.c = 0;<br>
> col[0].c = 0; val[0] = 1;<br>
> col[0].i = i; col[0].j = j;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,1,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
><br>
> //vy: //vy(i,j) - c*vy(i,j+-1) = 0;<br>
> row.c = 1;<br>
> col[0].c = 1; val[1] = coeff;<br>
> col[1].c = 1;<br>
> col[1].i = i;<br>
> if(j == 0) //vy(i,j) - c*vy(i,j+1) = 0;<br>
> col[1].j = j+1;<br>
> else //vy(i,j) - c*vy(i,j-1) = 0;<br>
> col[1].j = j-1;<br>
><br>
> ierr = MatSetValuesStencil(mat1,1,&row,2,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
> }<br>
> else if (i == 0 || i == info.my-1){ //top and bottom borders<br>
> //vx: vx(i,j) - c* vx(i+-1,j) = 0;<br>
> row.c = 0;<br>
> col[0].c = 0; val[0] = 1;<br>
> col[0].i = i; col[0].j = j;<br>
> col[1].c = 0; val[1] = coeff;<br>
> if (i == 0) //vx(i,j) - c*vx(i+1,j) = 0;<br>
> col[1].i = i+1;<br>
> else //vx(i,j) - c*vx(i-1,j) = 0;<br>
> col[1].i = i-1;<br>
> col[1].j = j;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,2,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
><br>
> //vy(i,j) = 0;<br>
> row.c = 1;<br>
> col[0].c = 1;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,1,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
> }<br>
> else { //Interior points:<br>
> row.c = 0;//x-eq<br>
> col[0].c = 0; val[0] = 2*coeff;<br>
> col[0].i = i; col[0].j = j+1;<br>
> col[1].c = 0; val[1] = -val[0] - coeff;<br>
> col[1].i = i; col[1].j = j;<br>
> col[2].c = 1; val[2] = 4*coeff;<br>
> col[2].i = i+1; col[2].j = j;<br>
> col[3].c = 1; val[3] = 4*coeff;<br>
> col[3].i = i; col[3].j = j;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,4,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
><br>
> row.c = 1; //y-eq<br>
> col[0].c = 1; val[0] = 2*coeff;<br>
> col[0].i = i; col[0].j = j+1;<br>
> col[1].c = 1; val[1] = -val[0] - coeff;<br>
> col[1].i = i; col[1].j = j;<br>
> col[2].c = 0; val[2] = 4*coeff;<br>
> col[2].i = i+1; col[2].j = j;<br>
> col[3].c = 0; val[3] = 4*coeff;<br>
> col[3].i = i; col[3].j = j;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,4,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
> }<br>
> }<br>
> }<br>
><br>
> MatAssemblyBegin(mat1,MAT_FINAL_ASSEMBLY);<br>
> MatAssemblyEnd(mat1,MAT_FINAL_ASSEMBLY);<br>
> ierr = MatSetOption(mat1,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);<br>
> ierr = PetscObjectSetName((PetscObject)mat1,"mat1");CHKERRQ(ierr);<br>
><br>
> /* clean up and exit */<br>
> ierr = DMDestroy(&da);CHKERRQ(ierr);<br>
> ierr = MatDestroy(&mat1);<br>
> ierr = PetscFinalize();<br>
> return 0;<br>
> }<br>
><br>
><br>
><br>
><br>
> PetscInt m = 10, n=10;<br>
> ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,<br>
> PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr);<br>
><br>
> Mat mat1;<br>
> ierr = DMCreateMatrix(da,MATMPIAIJ,&mat1); CHKERRQ(ierr);<br>
> MatStencil row, col[4]; //let's test 4 non-zeros in one row.<br>
> PetscScalar val[4];<br>
> PetscScalar coeff = 2.; //just a constant coefficient for testing.<br>
> PetscInt i,j;<br>
> DMDALocalInfo info;<br>
> ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);<br>
> //Now fill up the matrix:<br>
> for(i = info.ys; i < info.ys+info.ym; ++i){<br>
> for(j = info.xs; j < info.xs+info.xm; ++j){<br>
> row.i = i; row.j = j; //one node at a time<br>
> if (j == 0 || j == info.mx-1){ //left and right borders<br>
> //vx(i,j) = 0;<br>
> row.c = 0;<br>
> col[0].c = 0; val[0] = 1;<br>
> col[0].i = i; col[0].j = j;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,1,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
><br>
> //vy: //vy(i,j) - c*vy(i,j+-1) = 0;<br>
> row.c = 1;<br>
> col[0].c = 1; val[1] = coeff;<br>
> col[1].c = 1;<br>
> col[1].i = i;<br>
> if(j == 0) //vy(i,j) - c*vy(i,j+1) = 0;<br>
> col[1].j = j+1;<br>
> else //vy(i,j) - c*vy(i,j-1) = 0;<br>
> col[1].j = j-1;<br>
><br>
> ierr = MatSetValuesStencil(mat1,1,&row,2,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
> }<br>
> else if (i == 0 || i == info.my-1){ //top and bottom borders<br>
> //vx: vx(i,j) - c* vx(i+-1,j) = 0;<br>
> row.c = 0;<br>
> col[0].c = 0; val[0] = 1;<br>
> col[0].i = i; col[0].j = j;<br>
> col[1].c = 0; val[1] = coeff;<br>
> if (i == 0) //vx(i,j) - c*vx(i+1,j) = 0;<br>
> col[1].i = i+1;<br>
> else //vx(i,j) - c*vx(i-1,j) = 0;<br>
> col[1].i = i-1;<br>
> col[1].j = j;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,2,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
><br>
> //vy(i,j) = 0;<br>
> row.c = 1;<br>
> col[0].c = 1;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,1,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
> }<br>
> else { //Interior points:<br>
> row.c = 0;//x-eq<br>
> col[0].c = 0; val[0] = 2*coeff;<br>
> col[0].i = i; col[0].j = j+1;<br>
> col[1].c = 0; val[1] = -val[0] - coeff;<br>
> col[1].i = i; col[1].j = j;<br>
> col[2].c = 1; val[2] = 4*coeff;<br>
> col[2].i = i+1; col[2].j = j;<br>
> col[3].c = 1; val[3] = 4*coeff;<br>
> col[3].i = i; col[3].j = j;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,4,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
><br>
> row.c = 1; //y-eq<br>
> col[0].c = 1; val[0] = 2*coeff;<br>
> col[0].i = i; col[0].j = j+1;<br>
> col[1].c = 1; val[1] = -val[0] - coeff;<br>
> col[1].i = i; col[1].j = j;<br>
> col[2].c = 0; val[2] = 4*coeff;<br>
> col[2].i = i+1; col[2].j = j;<br>
> col[3].c = 0; val[3] = 4*coeff;<br>
> col[3].i = i; col[3].j = j;<br>
> ierr = MatSetValuesStencil(mat1,1,&row,4,col,val,INSERT_VALUES);CHKERRQ(ierr);<br>
> }<br>
> }<br>
> }<br>
><br>
> MatAssemblyBegin(mat1,MAT_FINAL_ASSEMBLY);<br>
> MatAssemblyEnd(mat1,MAT_FINAL_ASSEMBLY);<br>
><br>
> However I get the following error when run with 2 processors.<br>
> [0]PETSC ERROR: --------------------- Error Message ------------------------------------<br>
> [0]PETSC ERROR: Argument out of range!<br>
> [0]PETSC ERROR: Local index 120 too large 120 (max) at 0!<br>
> [0]PETSC ERROR: ------------------------------------------------------------------------<br>
> [0]PETSC ERROR: Petsc Release Version 3.4.1, Jun, 10, 2013<br>
> [0]PETSC ERROR: See docs/changes/index.html for recent updates.<br>
> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.<br>
> [0]PETSC ERROR: See docs/index.html for manual pages.<br>
> [0]PETSC ERROR: ------------------------------------------------------------------------<br>
> [0]PETSC ERROR: examples/FDdmda_test on a arch-linux2-cxx-debug named edwards by bkhanal Thu Jul 4 14:20:11 2013<br>
> [0]PETSC ERROR: Libraries linked from /home/bkhanal/Documents/softwares/petsc-3.4.1/arch-linux2-cxx-debug/lib<br>
> [0]PETSC ERROR: Configure run at Wed Jun 19 11:04:51 2013<br>
> [0]PETSC ERROR: Configure options --with-cc=gcc --with-fc=g77 --with-cxx=g++ --download-f-blas-lapack=1 --download-mpich=1 -with-clanguage=cxx --download-hypre=1<br>
> [0]PETSC ERROR: ------------------------------------------------------------------------<br>
> [0]PETSC ERROR: ISLocalToGlobalMappingApply() line 444 in /home/bkhanal/Documents/softwares/petsc-3.4.1/src/vec/is/utils/isltog.c<br>
> [0]PETSC ERROR: MatSetValuesLocal() line 1967 in /home/bkhanal/Documents/softwares/petsc-3.4.1/src/mat/interface/matrix.c<br>
> [0]PETSC ERROR: MatSetValuesStencil() line 1339 in /home/bkhanal/Documents/softwares/petsc-3.4.1/src/mat/interface/matrix.c<br>
> [0]PETSC ERROR: main() line 75 in "unknowndirectory/"/user/bkhanal/home/works/cmake_tuts/petsc_test/examples/FDdmda_test.cxx<br>
> application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0<br>
> [cli_0]: aborting job:<br>
> application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0<br>
><br>
> ===================================================================================<br>
> = BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES<br>
> = EXIT CODE: 63<br>
> = CLEANING UP REMAINING PROCESSES<br>
> = YOU CAN IGNORE THE BELOW CLEANUP MESSAGES<br>
> ===================================================================================<br>
><br>
> I do not understand how the index (i,j) are out of range when set to corresponding fields in row and col variables. What could be the possible problem ? And any suggestions on the way to debug this sort of issues ?<br>
><br>
> Thanks,<br>
> Bishesh<br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
> --<br>
> 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<br>
><br>
<br>
</div></div></blockquote></div><br></div></div>