[petsc-users] Reuse Matrix

[Baron Law]

Thanks Barry,

I originally plan to use the matrix for KSP with LSQR solver. Now what I try
to do is implement my own sparse matrix using Vec as I just need to store
the column index. Instead of using KSP, I try to implement my own regression
solver as the result beta[i], is simply the mean of all RHS y where the
non-zero column = i.

Best,
Baron

-----Original Message-----
From: Barry Smith [mailto:bsmith at mcs.anl.gov] 
Sent: Wednesday, April 08, 2015 3:23 PM
To: Baron Law
Cc: petsc-users at mcs.anl.gov
Subject: Re: [petsc-users] Reuse Matrix


   Baron,

     This is tricky. when the 1 moves between the "diagonal part" and the
"off-diagonal" part it will trigger a new allocation and thus slow things
down a great deal. I think the traditional PETSc matrix formats are not good
for your situation.

     I think it is likely much better to use a MATSHELL for this matrix.
What do you use the matrix for? You do matrix-vector product with it,
anything else? You can write your own matrix vector product with MATSHELL
customized for your problem.

  Barry


I cheat would be to have a single nonzero entry always in the diagonal part
AND the off-diagonal part of each row and just set one or the other to 0.0
each time so you always use the exact same matrix layout. But this will
induce unneeded communication since PETSc will need to set up the scatter
for that off diagonal entry even though it is zero.

> On Apr 8, 2015, at 1:39 PM, Baron Law <baronlaw at purdue.edu> wrote:
> 
> Hi,
>  
> I need to compute the solution of least square regression for many times
(> 10^6) for large matrix A (>10^9x10^3). I use a hypercube basis, so each
row of the matrix A has exactly one non-zero entry with value 1. On
different runs, the positions of the 1 are different. I would like to reuse
the sparse structure of A, what is the best way to do in this scenario?
>  
> I have heard of MatCreateMPIAIJWithSplitArrays, but it is kind of
difficult to use. I try it a bit and seems like after creation if I need to
change the column index, I need to skip the diagonal block. Also, I am not
100% sure, how to find which columns belong to the diagonal block for each
process. But as mentioned, I need to run the regression million times, even
it is difficult to use, I will try if it is fastest method.
>  
> Best,
> Baron