[petsc-users] The multiplication of large matrices

Wed May 29 15:42:16 CDT 2013

Joon hee Choi <choi240 at purdue.edu> writes:

> Dear Matthew,
>
> Thank you for your fast reply.
> I am attaching the actual error output and my code below. The error output is repeated by for-loop. Could you let me know what might be wrong? 
>
> Thank you,
> Joon
>
> [0]PETSC ERROR: --------------------- Error Message ------------------------------------
> [0]PETSC ERROR: No support for this operation for this object type!
> [0]PETSC ERROR: MatMatMult not supported for B of type localref!
> [0]PETSC ERROR: ------------------------------------------------------------------------
> [0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 6, Mon Feb 11 12:26:34 CST 2013
> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
> [0]PETSC ERROR: See docs/index.html for manual pages.
> [0]PETSC ERROR: ------------------------------------------------------------------------
> [0]PETSC ERROR: ./tensor on a linux-sta named rossmann-a009.rcac.purdue.edu by choi240 Wed May 29 05:53:46 2013
> [0]PETSC ERROR: Libraries linked from /apps/rhel5/petsc-3.3-p6/64/impi-4.1.0.024_intel-13.0.1.117_ind64/linux-static/lib
> [0]PETSC ERROR: Configure run at Tue May 21 15:56:45 2013
> [0]PETSC ERROR: Configure options --with-cc=mpiicc --with-cxx=mpiicpc --with-fc=mpiifort --with-scalar-type=real --with-shared-libraries=0 --with-pic=1 --with-clanguage=C++ --with-fortran --with-fortran-kernels=1 --with-64-bit-indices=1 --with-debugging=0 --with-blas-lapack-dir=/opt/intel/composer_xe_2013.1.117/mkl/lib/intel64 --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --download-hdf5=no --download-metis=no --download-parmetis=no --download-superlu_dist=no --download-mumps=no --download-scalapack=yes --download-blacs=yes --download-hypre=no --download-spooles=no
> [0]PETSC ERROR: ------------------------------------------------------------------------
> [0]PETSC ERROR: MatMatMult() line 8603 in /apps/rhel5/petsc-3.3-p6/64/impi-4.1.0.024_intel-13.0.1.117_ind64/src/mat/interface/matrix.c
> [0]PETSC ERROR: --------------------- Error Message ------------------------------------
> [0]PETSC ERROR: No support for this operation for this object type!
> [0]PETSC ERROR: MatMatMult not supported for B of type localref!
>
>    for (i=0; i<I; i++) idxXrow[i] = i;
>    for (i=0; i<I; i++) idxBrow[i] = i;
>    for (j=0; j<J; j++) idxCrow[j] = j;
>    ierr = ISCreateGeneral(PETSC_COMM_SELF, I, idxXrow, PETSC_COPY_VALUES, &isXrow); CHKERRQ(ierr);
>    ierr = ISCreateGeneral(PETSC_COMM_SELF, I, idxBrow, PETSC_COPY_VALUES, &isBrow); CHKERRQ(ierr);
>    ierr = ISCreateGeneral(PETSC_COMM_SELF, J, idxCrow, PETSC_COPY_VALUES, &isCrow); CHKERRQ(ierr);
>
>    for (r=0; r<R; r++) {
>       ierr = ISCreateGeneral(PETSC_COMM_SELF, 1, &r, PETSC_COPY_VALUES, &isBCcol); CHKERRQ(ierr);
>       ierr = MatGetLocalSubMatrix(B, isBrow, isBCcol, &tempB); CHKERRQ(ierr);
>       ierr = MatGetLocalSubMatrix(C, isCrow, isBCcol, &tempC); CHKERRQ(ierr);

>From the man page:

   Depending on the format of mat, the returned submat may not implement MatMult().  Its communicator may be
   the same as mat, it may be PETSC_COMM_SELF, or some other subcomm of mat's.

   The submat always implements MatSetValuesLocal().  If isrow and iscol have the same block size, then
   MatSetValuesBlockedLocal() will also be implemented.

>       ierr = MatAssemblyBegin(tempB, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>       ierr = MatAssemblyEnd(tempB, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>       ierr = MatAssemblyBegin(tempC, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>       ierr = MatAssemblyEnd(tempC, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
>       for (j=0; j<J; j++)
>          for (i=0; i<I; i++) idxXcol[i] = i + j*I;
>          ierr = ISCreateGeneral(PETSC_COMM_SELF, I, idxXcol, PETSC_COPY_VALUES, &isXcol); CHKERRQ(ierr);
>          ierr = MatGetLocalSubMatrix(X1, isXrow, isXcol, &tempX); CHKERRQ(ierr);
>          ierr = MatAssemblyBegin(tempX, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>          ierr = MatAssemblyEnd(tempX, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>          ierr = MatMatMult(tempX, tempB, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tempM);

It does not implement MatMult.  I don't see enough context in this email
about what you are trying to do, but MatGetLocalSubMatrix is not
intended for what you are doing.  MatGetSubMatrices() might do what you
want.

>          ierr = MatRestoreLocalSubMatrix(X1, isXrow, isXcol, &tempX); CHKERRQ(ierr);
>       }
>       ierr = MatRestoreLocalSubMatrix(B, isBrow, isBCcol, &tempB); CHKERRQ(ierr);
>       ierr = MatRestoreLocalSubMatrix(C, isCrow, isBCcol, &tempC); CHKERRQ(ierr);
>    }
>
> ----- Original Message -----
> From: "Matthew Knepley" <knepley at gmail.com>
> To: "Joon hee Choi" <choi240 at purdue.edu>
> Cc: petsc-users at mcs.anl.gov
> Sent: Wednesday, May 29, 2013 7:01:26 AM
> Subject: Re: [petsc-users] The multiplication of large matrices
>
>
> On Wed, May 29, 2013 at 6:43 AM, Joon hee Choi < choi240 at purdue.edu > wrote: 
>
>
>
>
> 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. 
>
>
>
> Look, we have rules on this list. Without them, we cannot help you. You MUST send the ENTIRE error output. 
> Without that, we are just guessing. 
>
>
>
>
> 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. 
>
>
>
> If you get an Out Of Memory error (which we need to see), it means that you are out of memory. 
>
>
> Matt 
>
>
> 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 
>
>
>
>
> -- 
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. 
> -- Norbert Wiener 