[petsc-dev] [petsc-users] matrix-matrix addition

[Barry Smith]

On Sep 27, 2010, at 4:28 PM, Shri wrote:

> Barry,
>      I've modified MatAXPY_SeqAIJ as follows when DIFFERENT_NONZERO_STRUCTURE flag is true (i) create a new matrix and do preallocation for it. For determining the number of nonzeros per row in the new matrix, i've used the number of nonzeros in the original two matrices. This is the piece of code which computes the nnz per row in the new matrix
> for(i = 0;i < m;i++) {
>    nzx = xi[i+1] - xi[i]; nzy = yi[i+1] - yi[i];
>    if (nzx > nzy) nnz[i] = nzy + (nzx - nzy);
>    else nnz[i] = nzx + (nzy - nzx);
> }
> 
>  (ii) Set the values in this new matrix by adding values of matrix Y first followed by matrix X (iii) Replace the old matrix by new one using MatHeaderReplace. 
> 
> Let me know if its ok and i'll add the mpiaij version.

   Note that the mpiaij one is much simpler; you only need to call the MatAXPY() on each of the two submatrices A and B.

> 
> Btw, i'm getting memory leak in the code.
> [0]Total space allocated 16 bytes
> [ 0]16 bytes PetscCommDuplicate() line 152 in src/sys/objects/tagm.c
>      [0]  PetscHeaderCreate_Private() line 25 in src/sys/objects/inherit.c
>      [0]  MatCreate() line 67 in src/mat/utils/gcreate.c
> 
> I suspect it might be due to MatHeaderReplace(). Any cues?

   Hmmm, does it it lead memory with the old code? Or only after the change.

  Barry

> 
> Boyce,can you please check if this working for your application.
> 
> Shri
> 
> ----- Boyce Griffith <griffith at cims.nyu.edu> wrote:
>> Hi, Barry et al. --
>> 
>> I'm moving this over to petsc-dev, since that seems like the more 
>> appropriate place.
>> 
>> I'm attaching my stab at implementing MatAXPY for pairs of sequantial 
>> AIJ matrices.  I've done only very minimal testing, but it appears to 
>> work for my application.
>> 
>> In the implementation, I create a new matrix W to store the result of 
>> computing a*X+Y.  Instead of swapping out the innards of W and Y, my 
>> implementation "returns" the newly allocated matrix W.  (I didn't see a 
>> function to swap out matrices.  Is there one that could be used here?)
>> 
>> -- Boyce
>> 
>> On 9/27/10 1:06 PM, Barry Smith wrote:
>>> 
>>> On Sep 27, 2010, at 12:04 PM, Jed Brown wrote:
>>> 
>>>> On Mon, Sep 27, 2010 at 19:01, Barry Smith<bsmith at mcs.anl.gov>  wrote:
>>>>> If you do not merge the column indices from A and B but simply count them all you will get over allocation, at worst a factor of two, though usually it would just be a little.
>>>> 
>>>> You count the number of *unique* entries.  You have two indices that
>>>> run through the arrays aj and bj, counting duplicates only once.  This
>>>> is easy to do since they are both sorted.
>>> 
>>>    Yes, sounds like it could work and definitely better than my suggestion.
>>> 
>>>   Shri, can you try this?
>>> 
>>>    Barry
>>> 
>>>> 
>>>> Jed
>>> 
>>> 
>>> 
>