On Fri, Jan 20, 2012 at 6:02 PM, Xuefei (Rebecca) Yuan <xyuan at lbl.gov>wrote:
> 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) > > > > > 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 > > > > > 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. > MatSetOption(MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE); Matt > Thanks very much! > > Have a nice weekend! > > Cheers, > > Rebecca > > > > > > > > > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener -------------- next part -------------- An HTML attachment was scrubbed... URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120120/2dcaefea/attachment.htm>