Cong Li <[email protected]> writes: > Thanks for the answer. > > Do you mean get each process's local data and save them in the C array > locally, then let the program call MPI_Gather to combine the arrrays on one > process? > > Actuall I am trying to implement block cg method for my research.
Okay, then you don't want the P^T A P matrix to ever be MPIAIJ. You can compute it by doing the local product and MPI_Allreduce'ing so that everyone has the k*k matrix. > It is like this AX=B, X and B are mxk matrices, A is mxm large sparse marix. > > Given initial guess X0, R0=B-AX > P=R0 > for i=1,2,...do > gamma(i)=inv(Pi(transpose)APi)Pi(transpose)Ri This should be a solve, not "inv". > X(i+1) = Xi +P(i) gamma(i) > ..... > .....psi(i) = - inv(Pi(transpose)APi)PiA(transpose)R(i+1) > ..... > > Since the program store Pi(transpose)APi in MATMPIAJI (beacase I am using > MatPtAP call and A is MATMPIAIJ, I have to save the result in MATMPIAIJ for > the compatibabilty concern), I need to convert Pi(transpose)APi to > MATSEQDENSE to get its inverse. I'm afraid PETSc does not have a native dense format that stores redundantly, but that's what you want. (In Elemental notation, this is a [*,*] distribution.) This could be added to PETSc and we can advise if you are interested in doing that, but if you want to do the minimal work to implement your method, you should do the product (A P), then get out the arrays and do the local part followed by an MPI_Allreduce. > By the way, I also think elemental is a good choice, but I don't how to > call its functions from PETSc, could you show me an example (data type > conversion from MATMPIAIJ to elemental type, and how can I call elemental > function in PETSc code) ? Your problem doesn't really benefit from Elemental since you really just need local operations. The 2d distributions that PETSc's elemental supports are not what you want for your purposes.
pgpEpP7dlOhzj.pgp
Description: PGP signature
