On Jan 21, 2012, at 10:27 AM, Matthew Knepley wrote:

> 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?

I do not want to use DMGetMatrix, so I use the MatCreate() to get my own 
matrix. But inside the call DMMGSetSNESLocal(), this DMGetMatrix was called and 
the matrix is set to 5776 nonzeros...


> 
>   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/29fe6855/attachment.htm>

Reply via email to