[petsc-users] The matrix allocated by calling MatMPIAIJSetProallocation() was flushed out after calling DMMGSetSNESLocal()
Matthew Knepley
knepley at gmail.com
Sat Jan 21 12:27:30 CST 2012
On Sat, Jan 21, 2012 at 12:19 PM, Xuefei (Rebecca) Yuan <xyuan at lbl.gov>wrote:
> Hello,
>
> After I was not able to get rid of the zeros in the matrix by
>
> ierr = MatSetOption(jacobian, MAT_IGNORE_ZERO_ENTRIES,
>>> PETSC_TRUE);CHKERRQ(ierr);
>>>
>>
>
> I would like to create the matrix by the following routine:
>
> ierr = MatCreate(PETSC_COMM_WORLD, &jacobian);CHKERRQ(ierr);
> ierr = MatSetType(jacobian, MATMPIAIJ);CHKERRQ(ierr);
> ierr = MatSetSizes(jacobian, PETSC_DECIDE, PETSC_DECIDE, (PetscInt)(
> info.mx*info.my*4), (PetscInt)(info.mx*info.my*4));CHKERRQ(ierr);
> ierr = MatMPIAIJSetPreallocation(jacobian, 11, PETSC_NULL, 18,
> PETSC_NULL);CHKERRQ(ierr);
>
> instead of using
>
> ierr = DMGetMatrix(DMMGGetDM(dmmg), MATAIJ, &jacobian);CHKERRQ(ierr);
>
> to save the memory. However, I found that the routine DMGetMatrix() was
> still called inside DMMGSetSNESLocal(), therefore, the matrix I created was
> useless:
>
What exactly are you getting out of using the DM? If it does not express
the pattern of your problem, why use it?
Matt
> *********************************************
>
> Breakpoint 2, main (argc=3, argv=0x7fff5fbff770) at twcartffxmhd.c:255
> 255 ierr = MatCreate(PETSC_COMM_WORLD, &jacobian);CHKERRQ(ierr);
> (gdb) n
> 256 ierr = MatSetType(jacobian, MATMPIAIJ);CHKERRQ(ierr);
> (gdb)
> 257 ierr = MatSetSizes(jacobian, PETSC_DECIDE, PETSC_DECIDE, (PetscInt)(
> info.mx*info.my*4), (PetscInt)(info.mx*info.my*4));CHKERRQ(ierr);
> (gdb)
> 258 ierr = MatMPIAIJSetPreallocation(jacobian, 11, PETSC_NULL, 18,
> PETSC_NULL);CHKERRQ(ierr);
> (gdb)
>
> Breakpoint 3, main (argc=3, argv=0x7fff5fbff770) at twcartffxmhd.c:260
> 260 ierr = DMMGSetSNESLocal(dmmg, FormFunctionLocal,
> FormJacobianLocal,0,0);CHKERRQ(ierr);
> (gdb) s
> DMMGSetSNESLocal_Private (dmmg=0x102004480, function=0x100016ba9
> <FormFunctionLocal>, jacobian=0x100050e83 <FormJacobianLocal>,
> ad_function=0, admf_function=0) at damgsnes.c:932
> 932 PetscInt i,nlevels = dmmg[0]->nlevels;
> (gdb) n
> 934 PetscErrorCode
> (*computejacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*) = 0;
> (gdb) s
> 937 PetscFunctionBegin;
> (gdb) n
> 938 if (jacobian) computejacobian = DMMGComputeJacobian;
> (gdb)
> 942 CHKMEMQ;
> (gdb)
> 943 ierr = PetscObjectGetCookie((PetscObject)
> dmmg[0]->dm,&cookie);CHKERRQ(ierr);
> (gdb)
> 944 if (cookie == DM_COOKIE) {
> (gdb)
> 948 ierr = PetscOptionsHasName(PETSC_NULL,
> "-dmmg_form_function_ghost", &flag);CHKERRQ(ierr);
> (gdb)
> 949 if (flag) {
> (gdb)
> 952 ierr =
> DMMGSetSNES(dmmg,DMMGFormFunction,computejacobian);CHKERRQ(ierr);
> (gdb)
> Matrix Object:
> type=seqaij, rows=100, cols=100
> total: nonzeros=5776, allocated nonzeros=5776
> using I-node routines: found 25 nodes, limit used is 5
> 954 for (i=0; i<nlevels; i++) {
> (gdb) b twcartffxmhd.c:280
> Breakpoint 4 at 0x100004032: file twcartffxmhd.c, line 282.
> (gdb) c
> Continuing.
> Matrix Object:
> type=mpiaij, rows=100, cols=100
> total: nonzeros=0, allocated nonzeros=2900
> using I-node (on process 0) routines: found 20 nodes, limit used is 5
>
> *********************************************
>
> The final matrix after FormJacobianlLocal() call still have 5776 nonzeros:
>
> % Size = 100 100
> 2 % Nonzeros = 5776
> 3 zzz = zeros(5776,3);
>
> Is there a way that I can pass the matrix created via line 255-258 to
> FormJacobianLocal()?
>
> Thanks very much!
>
> Cheers,
>
> Rebecca
>
>
>
>
>
>
>
> On Jan 20, 2012, at 5:23 PM, Matthew Knepley wrote:
>
> On Fri, Jan 20, 2012 at 7:07 PM, Xuefei (Rebecca) Yuan <xyuan at lbl.gov>wrote:
>
>> Here is the output:
>>
>>
>>
>>
>>
>> On Jan 20, 2012, at 5:01 PM, Matthew Knepley wrote:
>>
>> On Fri, Jan 20, 2012 at 6:55 PM, Xuefei (Rebecca) Yuan <xyuan at lbl.gov>wrote:
>>
>>> Hello Matt,
>>>
>>> I tried several times for 3.1-p8 and dev version by putting MatSetOption
>>>
>>
>> Are you sure your entries are exactly 0.0?
>>
>>
> Are you using ADD_VALUES?
>
>
> http://petsc.cs.iit.edu/petsc/petsc-dev/file/783e93230143/src/mat/impls/aij/seq/aij.c#l310
>
> Matt
>
>
>> Matt
>>
>>
>>> 1) right after creation of the matrix:
>>>
>>> #ifdef petscDev
>>> ierr = DMCreateMatrix(DMMGGetDM(dmmg), MATAIJ,
>>> &jacobian);CHKERRQ(ierr);
>>> #else
>>> ierr = DAGetMatrix(DMMGGetDA(dmmg), MATAIJ,
>>> &jacobian);CHKERRQ(ierr);
>>> #endif
>>> ierr = MatSetOption(jacobian, MAT_IGNORE_ZERO_ENTRIES,
>>> PETSC_TRUE);CHKERRQ(ierr);
>>>
>>> 2) at the beginning of the FormJacobianLocal() routine:
>>>
>>> PetscFunctionBegin;
>>> ierr = MatSetOption(jacobian, MAT_IGNORE_ZERO_ENTRIES,
>>> PETSC_TRUE);CHKERRQ(ierr);
>>>
>>> 3) before calling MatAssemblyBegin() in FormJacobianLocal() routine:
>>>
>>> ierr = MatSetOption(jacobian, MAT_IGNORE_ZERO_ENTRIES,
>>> PETSC_TRUE);CHKERRQ(ierr);
>>> ierr = MatAssemblyBegin(jacobian,
>>> MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>>> ierr = MatAssemblyEnd(jacobian,
>>> MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>>>
>>> None of those works. What is wrong here?
>>>
>>> Thanks,
>>>
>>> R
>>>
>>>
>>> On Jan 20, 2012, at 4:32 PM, Matthew Knepley wrote:
>>>
>>> On Fri, Jan 20, 2012 at 6:28 PM, Xuefei (Rebecca) Yuan <xyuan at lbl.gov>wrote:
>>>
>>>> Hello Matt,
>>>>
>>>> I have changed the code as
>>>>
>>>> ierr = MatSetOption(jacobian, MAT_IGNORE_ZERO_ENTRIES,
>>>> PETSC_TRUE);CHKERRQ(ierr);
>>>> ierr = MatAssemblyBegin(jacobian,
>>>> MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>>>> ierr = MatAssemblyEnd(jacobian,
>>>> MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>>>>
>>>
>>> You have to set it before you start setting values, so we know to ignore
>>> them.
>>>
>>> Matt
>>>
>>>
>>>> but still get the same result as before, the matrix still has 5776
>>>> nonzeros:
>>>>
>>>> % Size = 100 100
>>>> 2 % Nonzeros = 5776
>>>> 3 zzz = zeros(5776,3);
>>>>
>>>> Then I switch the order as
>>>>
>>>> ierr = MatAssemblyBegin(jacobian,
>>>> MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>>>> ierr = MatAssemblyEnd(jacobian,
>>>> MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>>>> ierr = MatSetOption(jacobian, MAT_IGNORE_ZERO_ENTRIES,
>>>> PETSC_TRUE);CHKERRQ(ierr);
>>>>
>>>> nothing changed.
>>>>
>>>> The version is 3.1-p8.
>>>>
>>>> Thanks very much!
>>>>
>>>> Best regards,
>>>>
>>>> Rebecca
>>>>
>>>>
>>>>
>>>>
>>>> On Jan 20, 2012, at 4:09 PM, Matthew Knepley wrote:
>>>>
>>>> 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
>>>>
>>>>
>>>>
>>>
>>>
>>> --
>>> 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
>>>
>>>
>>>
>>
>>
>> --
>> 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
>>
>>
>>
>>
>
>
> --
> 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
>
>
>
--
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/20120121/491d9476/attachment-0001.htm>
More information about the petsc-users
mailing list