[petsc-users] forming a matrix from a set of vectors

[Barry Smith]

   Nicolas,

    I believe the best way to do this is to write a code that specifically knows the SeqAIJ matrix structure  and does the product directly with that data structure instead of trying to "cook something up" using the standard higher level PETSc routines. So you need to have 

#include <../src/mat/impls/aij/seq/aij.h>  

directly in your code. See for example src/mat/impls/aij/seq/aij.c 

Next you need to store c = transpose(d) since storing the matrix d in PETSc sparse format (which is row based) is terrible for the operation you want to perform. So ideally just create the matrix c and put its values in; if that is too difficult you can use MatTranspose() to get c from d.

   Now to the algorithm 

    T =     c *(A*transpose(c)) 

note that the columns of transpose(c) (i.e. the columns of d)  are the rows of c, let c(i,*) represent a row of c. The diagonals of T are 

    T(i,i) = c(i,*) * (A * c(i,*))

   Let s = A * c(i,*) now to implement this just 
loop over i. Compute s (an array) using directly the data in the sparse matrix c (use MatDenseGetArray()) to access the values in A) then compute 
c(i,*) * s; then move on to the next i.   This will be very efficient (computes only exactly what is needed) and requires only the array s which is small (size of number of rows of A, not c).

  Barry





> On Aug 26, 2015, at 4:01 PM, Nicolas Pozin <nicolas.pozin at inria.fr> wrote:
> 
> Yes, this is sequentially and yes again, d is stored as a AIJ matrix
> 
> ----- Mail original -----
>> De: "Barry Smith" <bsmith at mcs.anl.gov>
>> À: "Nicolas Pozin" <nicolas.pozin at inria.fr>
>> Cc: "Jed Brown" <jed at jedbrown.org>, petsc-users at mcs.anl.gov
>> Envoyé: Mercredi 26 Août 2015 22:58:12
>> Objet: Re: [petsc-users] forming a matrix from a set of vectors
>> 
>> 
>>   Since A is tiny I am assuming you are doing this sequentially only?
>> 
>>   Do you have d stored as a AIJ matrix or is transpose(d) stored as a AIJ
>>   matrix?
>> 
>> 
>> 
>> 
>>> On Aug 26, 2015, at 3:41 PM, Nicolas Pozin <nicolas.pozin at inria.fr> wrote:
>>> 
>>> Actually I want to get the diagonal of the matrix : transpose(d)*A*d where
>>> -d is a sparse matrix of size (n1,m1)
>>> -A is a dense symetric matrix of size size (n1,n1)
>>> with m1 very big compared to n1 (1 million against a few dozens).
>>> 
>>> The problem is too big to allow the use of MatMatMult.
>>> What I planned to do :
>>> -compute the vectors Vi defined by transpose(d)*Ai where Ai is the i-th
>>> column of A : quick since d is sparse and n1 is small
>>> -deduce the matrix transpose(d)*A = [V1 ... Vn]
>>> and then get the diagonal of transpose(d)*transpose([V1 ...Vn]) through
>>> -transpose([V1 ...Vn]) and get its columns C1 ... Cn
>>> -conclude on the i-th diagonal value which is the i-th component of
>>> tranpose(d)*Ci
>>> 
>>> 
>>> 
>>> ----- Mail original -----
>>>> De: "Barry Smith" <bsmith at mcs.anl.gov>
>>>> À: "Nicolas Pozin" <nicolas.pozin at inria.fr>
>>>> Cc: "Jed Brown" <jed at jedbrown.org>, petsc-users at mcs.anl.gov
>>>> Envoyé: Mercredi 26 Août 2015 22:21:04
>>>> Objet: Re: [petsc-users] forming a matrix from a set of vectors
>>>> 
>>>> 
>>>>> On Aug 26, 2015, at 3:06 PM, Nicolas Pozin <nicolas.pozin at inria.fr>
>>>>> wrote:
>>>>> 
>>>>> Thank you for this answer.
>>>>> 
>>>>> What I want to do is to get the lines of this matrix and store them in
>>>>> vectors.
>>>> 
>>>> If you want to treat the columns of the dense matrix as vectors then use
>>>> MatDenseGetArray() and call VecCreateMPIWithArray() with a pointer to the
>>>> first row of each column of the obtained array (PETSc dense matrices are
>>>> stored by column; same as for example LAPACK).
>>>> 
>>>> But if you explained more why you want to treat something sometimes as a
>>>> Mat (which is a linear operator on vectors) and sometimes as vectors we
>>>> might be able to suggest how to organize your code.
>>>> 
>>>>  Barry
>>>> 
>>>>> 
>>>>> 
>>>>> ----- Mail original -----
>>>>>> De: "Jed Brown" <jed at jedbrown.org>
>>>>>> À: "Nicolas Pozin" <nicolas.pozin at inria.fr>, petsc-users at mcs.anl.gov
>>>>>> Envoyé: Mercredi 26 Août 2015 20:38:37
>>>>>> Objet: Re: [petsc-users]  forming a matrix from a set of vectors
>>>>>> 
>>>>>> Nicolas Pozin <nicolas.pozin at inria.fr> writes:
>>>>>>> Given a set of vectors V1, V2,...,Vn, is there an efficient way to form
>>>>>>> the
>>>>>>> dense matrix [V1 V2 ... Vn]?
>>>>>> 
>>>>>> What do you want to do with that matrix?  The vector representation is
>>>>>> pretty flexible and the memory semantics are similar unless you store
>>>>>> the dense matrix row-aligned (not the default).
>>>>>> 
>>>> 
>>>> 
>> 
>>