Hello all,
This is a test for np=1 case of the nonzero structure of the jacobian matrix.
The jacobian matrix is created via
ierr =
DMDACreate2d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
parameters.mxfield, parameters.myfield, PETSC_DECIDE, PETSC_DECIDE, 4, 2, 0, 0,
&da);CHKERRQ(ierr);
ierr = DMCreateMatrix(DMMGGetDM(dmmg), MATAIJ, &jacobian);CHKERRQ(ierr);
After creation of the jacobian matrix,
ierr = MatAssemblyBegin(jacobian,
MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jacobian,
MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
PetscViewer viewer;
char fileName[120];
sprintf(fileName, "jacobian_after_creation.m");CHKERRQ(ierr);
FILE * fp;
ierr =
PetscViewerASCIIOpen(PETSC_COMM_WORLD,fileName,&viewer);CHKERRQ(ierr);
ierr =
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
ierr = MatView (jacobian, viewer); CHKERRQ (ierr);
ierr = PetscFOpen(PETSC_COMM_WORLD,fileName,"a",&fp);
CHKERRQ(ierr);
ierr =
PetscViewerASCIIPrintf(viewer,"spy((spconvert(zzz)));\n");CHKERRQ(ierr);
ierr = PetscFClose(PETSC_COMM_WORLD,fp);CHKERRQ(ierr);
PetscViewerDestroy(&viewer);
I took a look at the structure of the jacobian by storing it in the matlab
format, the matrix has 5776 nonzeros entries, however, those values are all
zeros at the moment as I have not insert or add any values into it yet, the
structure shows: (the following figure shows a global replacement of 0.0 by 1.0
for those 5776 numbers)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: jacobian_after_creation.eps
Type: image/eps
Size: 112025 bytes
Desc: not available
URL:
<http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120120/026bc2ca/attachment-0002.bin>
-------------- next part --------------
Inside the FormJacobianLocal() function, I have selected the index to pass to
the nonzero values to jacobian, for example,
ierr = MatSetValuesStencil(jacobian, 1, &row, 6, col, v,
INSERT_VALUES);CHKERRQ(ierr);
where
col[0].i = column[4].i;
col[1].i = column[5].i;
col[2].i = column[6].i;
col[3].i = column[9].i;
col[4].i = column[10].i;
col[5].i = column[12].i;
col[0].j = column[4].j;
col[1].j = column[5].j;
col[2].j = column[6].j;
col[3].j = column[9].j;
col[4].j = column[10].j;
col[5].j = column[12].j;
col[0].c = column[4].c;
col[1].c = column[5].c;
col[2].c = column[6].c;
col[3].c = column[9].c;
col[4].c = column[10].c;
col[5].c = column[12].c;
v[0] = value[4];
v[1] = value[5];
v[2] = value[6];
v[3] = value[9];
v[4] = value[10];
v[5] = value[12];
and did not pass the zero entries into the jacobian matrix. However, after
inserting or adding all values to the matrix, by the same routine above to take
a look at the jacobian matrix in matlab format, the matrix still has 5776
nonzeros, in which 1075 numbers are nonzeros, and the other 4701 numbers are
all zeros. The spy() gives
-------------- next part --------------
A non-text attachment was scrubbed...
Name: jacobian_after_assembly.eps
Type: image/eps
Size: 28032 bytes
Desc: not available
URL:
<http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120120/026bc2ca/attachment-0003.bin>
-------------- next part --------------
for the true nonzero structures.
But the ksp_view will give the nonzeros number as 5776, instead of 1075:
linear system matrix = precond matrix:
Matrix Object: Mat_0x84000000_1 1 MPI processes
type: seqaij
rows=100, cols=100
total: nonzeros=5776, allocated nonzeros=5776
It is a waste of memory to have all those values of zeros been stored in the
jacobian.
Is there anyway to get rid of those zero values in jacobian and has the only
nonzero numbers stored in jacobian? In such a case, the ksp_view will tell that
total: nonzeros=1075.
Thanks very much!
Have a nice weekend!
Cheers,
Rebecca
