[petsc-users] The multiplication of large matrices

Wed May 29 15:14:59 CDT 2013

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);
      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);
         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 