[petsc-users] The multiplication of large matrices

[Joon hee Choi]

Hello all,

I am trying to compute the multiplication of large matrices using petsc. First,

I=2.6*10^7, J=4.8*10^7.
The matrix X is as follows: size of I*(I*J), block size of I, and non-zeros of 1.4*10^8.
-> X=[X1 X2 ...]
The matrix B is as follows: size of I*1, and dense.
The matrix M is as follows: size of I*J. M is the multiplication of each block of the matrix X and the matrix B. 
-> M=[X1*B X2*B ...]

I have to get the matrix M, given X and B. I successfully set up very large aij matrix X. However, I don't know what I have to do next. I tried two ways, but I failed.

First, I tried to get each block X1, X2,... using ISCreate and MatGetLocalSubMatrix, compute X1*B, X2*B,... using MatMatMult, and then set up the matrix M using ISLocalToGlobalMappingCreate and MatSetLocalToGlobalMapping. However, I got "No support for this operation for this object type" error. Also, this code was very slow because of 48 million loops.

Second, I tried to set up the new matrix BB from B. The matrix BB has the size of (I*J)*J and the block size of I*1. And every diagonal block of BB is B and all other blocks are 0 matrices. That is,

     | B 0 0 .. 0 |
BB = | 0 B 0 .. 0 |
     | ...        |
     | 0 0 0 .. B |

This is not a block diagonal matrix because B is not square. Anyway, if I set up BB, I can get M easily because M=X*BB. However, I got "out of memory" error from MatCreateSeqAIJ. I think this is because BB has non-zeros of I*J=10^15.

I never have any other ideas. If someone can fix my wrong ways correctly or has new ideas, then please please let me know. Thank you very much.

Joon