This was educational. One of the problems with depending on external
software is that there is a need for understanding everything
correctly or else it will come back to bite you. A lesson I wont
forget now.

The issue was neither one of the cases I specified earlier.

I was originally creating a block diagonal matrix by using the nnz
structure of a single variable. The problem occurs because I was
calling close() on the matrix every time I updated the block matrix
(variable block). I was looking specifically on nnz allocation but
failed to realize that when I call AssemblyEnd() on a matrix, it also
deallocates matrix extra storage. And the next time I add the block
diagonal elements for the next variable, I was forcing a reallocation
of non-zero elements. And hence in total, there was nvar-1
deallocations and nvar-1 allocations.

The way my code interacts with libMesh made it a little hard to get to
the problem. I can actually use the nnz/n_vars strategy to get the
correct allocation now. I apologize for all of your time and thanks
for your help in getting to the core of the issue :)

Vijay

On Fri, Dec 4, 2009 at 8:43 PM, Vijay S. Mahadevan <[email protected]> wrote:
> Yes, the runs were in serial. hmm. I just remembered that I also
> supply my CouplingMatrix for the system matrix and so the nnz should
> probably be divided by sum(CouplingMatrix(ivar,:)) where ivar would be
> the variable number. I hope that makes sense. I'll check this again
> tomorrow or later tonight and see if that fixes it.
>
> Or I could also be creating an expanded stencil in my physics based
> preconditioning that is triggering the extra mallocs. This should be
> evident if I spy the matrix in matlab and compare with my sparsity for
> matrix in the code.
>
> On Fri, Dec 4, 2009 at 8:14 PM, Kirk, Benjamin (JSC-EG311)
> <[email protected]> wrote:
>> That will be interesting.  Is this in serial?  In parallel you have 
>> 'on-processor' and 'off-processor' nonzeros, the total number in a given row 
>> is the sum of the two...
>>
>>
>>
>> ----- Original Message -----
>> From: Vijay S. Mahadevan <[email protected]>
>> To: Kirk, Benjamin (JSC-EG311)
>> Cc: [email protected] <[email protected]>; 
>> [email protected] <[email protected]>
>> Sent: Fri Dec 04 20:06:29 2009
>> Subject: Re: [Libmesh-users] Block matrix nnz's
>>
>> Yes they are. I did divide the total number of nonzeros in the row by
>> the number of variables but I still get lots of malloc hits during the
>> first assembly. And that's why I posed this question. I'll have to
>> print out the matrix and load it in matlab to see how my nnz
>> allocation differs from the actual matrix sparsity. I'll do this
>> tomorrow and let you know if there was something wrong in my
>> assumption.
>>
>> On Fri, Dec 4, 2009 at 7:24 PM, Kirk, Benjamin (JSC-EG311)
>> <[email protected]> wrote:
>>> Are all your variables of the same type?  If not, the concept of a block is 
>>> more tedious...  If so, can't you just divide the total number of nonzeros 
>>> in the row by the number of variables in the system?
>>>
>>> -Ben
>>>
>>>
>>> ----- Original Message -----
>>> From: Vijay S. Mahadevan <[email protected]>
>>> To: Derek Gaston <[email protected]>
>>> Cc: [email protected] 
>>> <[email protected]>
>>> Sent: Fri Dec 04 14:32:04 2009
>>> Subject: Re: [Libmesh-users] Block matrix nnz's
>>>
>>>> Believe it or not... Dr. Ben Kirk has thought of everything... and there is
>>>> no reason for us lesser beings to burn brain cycles over this kind of thing
>>>> ;-)
>>>
>>> I have come to the same conclusion several times before too :-)
>>>
>>> Derek, I have used the CouplingMatrix option before but was thinking
>>> along the lines of say calling DofMap with CouplingMatrix twice. One
>>> with the right coupling and one with a reduced coupling strategy
>>> (store only one variable block). This way, all true matvecs happen
>>> using the true coupled matrix while the reduced size matrix is used
>>> for performing preconditioning operations. There are some neat physics
>>> based preconditioners that can use these single blocks effectively
>>> (reduced storage and sometimes symmetric even if system matrix is
>>> block unsymmetric) and depending on the problem, can work better than
>>> ILU or some other numerical preconditioner.
>>>
>>> I got a reply from Rahul Sampath and his suggestion was to add a new
>>> ImplicitSystem with one variable to get the right nnz distribution for
>>> the single block. That seems reasonable and clean for now. Are there
>>> other alternatives ?
>>>
>>> Vijay
>>>
>>> On Fri, Dec 4, 2009 at 1:32 PM, Derek Gaston <[email protected]> wrote:
>>>> Believe it or not... Dr. Ben Kirk has thought of everything... and there is
>>>> no reason for us lesser beings to burn brain cycles over this kind of thing
>>>> ;-)
>>>>
>>>> In this case what you need to do is create a CouplingMatrix
>>>> (http://libmesh.sourceforge.net/doxygen/classCouplingMatrix.php) and attach
>>>> it to the DofMap _before_ compute_sparsity() gets called...
>>>>
>>>> You attach it to the DofMap by setting a public member:
>>>>
>>>> http://libmesh.sourceforge.net/doxygen/classDofMap.php#4d50b4c013a93b8036c126ced15ecbe6
>>>>
>>>> Derek
>>>>
>>>> On Dec 4, 2009, at 12:27 AM, Vijay S. Mahadevan wrote:
>>>>
>>>>> Hi,
>>>>>
>>>>> I'm trying to figure out a clean way to find the nnz sparsity
>>>>> structure for a given variable. That is given an EquationSystem and
>>>>> its corresponding Dofmap, can I find the correct nnz structure to
>>>>> create a matrix ignoring the coupling between variables. It is useful
>>>>> for creating block preconditioners where the entire matrix (with all
>>>>> coupling between variables) is not needed but just individual blocks
>>>>> would be enough to perform Block-Jacobi type operations for the
>>>>> preconditioner.
>>>>>
>>>>> I can always sweep through all the elements, compute the dofs for the
>>>>> variable of interest and create this nnz but that seems hackish. This
>>>>> information should definitely be derived from DofMap but looking at
>>>>> the DofMap code, I could not figure out how this can be done correctly
>>>>> either. Is there a simpler and efficient way to do this ?
>>>>>
>>>>> Vijay
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>> Join us December 9, 2009 for the Red Hat Virtual Experience,
>>>>> a free event focused on virtualization and cloud computing.
>>>>> Attend in-depth sessions from your desk. Your couch. Anywhere.
>>>>> http://p.sf.net/sfu/redhat-sfdev2dev
>>>>> _______________________________________________
>>>>> Libmesh-users mailing list
>>>>> [email protected]
>>>>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>>
>>>>
>>>
>>> ------------------------------------------------------------------------------
>>> Join us December 9, 2009 for the Red Hat Virtual Experience,
>>> a free event focused on virtualization and cloud computing.
>>> Attend in-depth sessions from your desk. Your couch. Anywhere.
>>> http://p.sf.net/sfu/redhat-sfdev2dev
>>> _______________________________________________
>>> Libmesh-users mailing list
>>> [email protected]
>>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>>>
>>
>

------------------------------------------------------------------------------
Join us December 9, 2009 for the Red Hat Virtual Experience,
a free event focused on virtualization and cloud computing. 
Attend in-depth sessions from your desk. Your couch. Anywhere.
http://p.sf.net/sfu/redhat-sfdev2dev
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users

Reply via email to