A^TA + B?

[Zhifeng Sheng]

Hong Zhang wrote:

>
> Zhifeng,
>
>>
>> I would like to implement B = A^TA + B.
>>
>> And the transpose multiplication can be done with the following 
>> function,
>>
>> PetscErrorCode PETSCMAT_DLLEXPORT MatMatMultTranspose(Mat A,Mat 
>> B,MatReuse scall,PetscReal fill,Mat *C) since I would like to save 
>> some memory, I don't want to create an intermediate matrix. In the 
>> above function, the discription of
>>
>> the argument "scall" is not clear, I know that MAT_INITIAL_MATRIX 
>> will create a new empty matrix and the result will be stored.
>
>
> Yes.
>
>>
>> what about MAT_REUSE_MATRIX ? does it mean that the C matrix will be 
>> added with the result?
>
>
> When the product matrix C=A^T*B is previously computed,
> and you want repeating C=A^T*B, in which A and B maintain
> the same non-zero sparse pattern but with different numerical values, 
> using
> MAT_REUSE_MATRIX will skip the symbolic computation of C=A^T*B
> and reuse the exiting matrix data structure and memory of C.
>
> Why do you need B = A^TA + B?
> If your matrices are dense, you should use lapack or scalapack.
> Sparse matrix product should be avoided because
> the product matrix C usually is much denser than A and B,
> and C=A^T*B cannot be implemented efficiently in general.
> We provide MatMatMultTranspose() mainly for multigrid computation
> in which the matrices have special data structure.
>
> Best,
>
> Hong
>
Thanks for your reply.

I need to implement a least-squares method, in which A^T*A will be 
computed, however, since I use finite elements, local matrices will be 
computed and then added to the global matrix. The A matrix is very 
sparse, e.g 9 non-zeros in a row.

Actually, in this case, I only need to do b = a^T*a + b on local 
matrices, which are very small e.g (12 by 12). The elementary matrices 
will be assembled in the end. Should I use lapack in combination with Petsc?

Best regards
Zhifeng