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