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
