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
