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

[Shri]

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.

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?

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