How to get the direct sum of mass matrix and itself

In Python, you can do it with the block_kronecker function in cbc.block (
https://launchpad.net/cbc.block),
which uses the Python interface of Dolfin and Trilinos.

If you want do it in C++ it is probably best to cast the matrix to its PETSc

or Trilinos type and let PETSc or Trilinos perform the sum.

Kent

On 31 August 2011 17:55, Nguyen Van Dang <
<email address hidden>> wrote:

> New question #169703 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/169703
>
> Hello,
> I would like to take the direct sum of two matrices.
> Ex:
> A=
> | 1 2 3 |
> | 2 3 1 |
> | 3 2 1 |
>
> B = A \oplus A=
>
> | 1 2 3 0 0 0 |
> | 2 3 1 0 0 0 |
> | 3 2 1 0 0 0 |
> | 0 0 0 1 2 3 |
> | 0 0 0 2 3 1 |
> | 0 0 0 3 2 1 |
>
> In my problem, I need the direct sum of the mass matrix and itself in C++.
> I hope that this sum is also a sparse matrix. Is there any way to do this?
> Thanks in advance.
> Best regards.
> Dang
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Dear Kent,
I still have problem when I use 'solve' function for direct-sum matrix. I send you my code as well as the errors I got. Would you like to help me again?
Thanks in advance.
Dang

My code:
#include <dolfin.h>
#include "mass_matrix.h"

using namespace dolfin;

int main()
{
  #ifdef HAS_SLEPC

  // Create mesh
  const double L1=1.0;
  const double L2=1.0;
  const double nx=2;
  const double ny=2;
  Rectangle mesh(-L1,-L2,L1,L2,nx,ny);

    // Build stiffness matrix
  PETScMatrix A;
  mass_matrix::FunctionSpace V(mesh);
  mass_matrix::BilinearForm a(V, V);
  assemble(A, a);
  int neqn = A.size(0);
  // variables
    Array<double> x1(neqn),x2(2*neqn);
    PETScVector rhs1(neqn),rhs2(2*neqn);
    PETScVector Y1(neqn),Y2(2*neqn);

    for (int j=0;j<neqn;j++)
  {
 x1[j] = j;
        x2[j] = j;
        x2[j+neqn] = j;
  }
  rhs1.set_local(x1);
  rhs2.set_local(x2);
  solve(A,Y1,rhs1);
  for (int j=0;j<neqn;j++)
        {
             printf(" %f",Y1[j]);
        }
        printf("\n");
  dolfin::uint row=0;
  std::vector<dolfin::uint> columns;
  std::vector<double> values;
  PETScMatrix M=A;
  M.resize(2*neqn,2*neqn);
  for (row=0;row<neqn;row++)
  {
   A.getrow(row, columns, values);
        M.setrow(row, columns, values);
        for (int j=0;j<columns.size();j++)
        {
             columns[j] = columns[j] + neqn;
        }
         M.setrow(row+neqn, columns, values);
  }

  solve(M,Y2,rhs2);

  #else

    cout << "SLEPc must be installed to run this demo." << endl;

  #endif

  return 0;
}

errors:

nguyenvandang@ubuntu:~/Documents/time-dept-prob/direct-sum$ ./direct-sum
 2.571429 -2.571429 19.714286 -2.571429 -2.571429 11.142857 67.714286 11.142857 43.714286
[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Object is in wrong state!
[0]PETSC ERROR: Not for unassembled matrix!
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Petsc Release Version 3.1.0, Patch 5, Mon Sep 27 11:51:54 CDT 2010
[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: Unknown Name on a linux-gnu named ubuntu by nguyenvandang Wed Aug 24 22:15:51 2011
[0]PETSC ERROR: Libraries linked from /build/buildd/petsc-3.1.dfsg/linux-gnu-c-opt/lib
[0]PETSC ERROR: Configure run at Mon Mar 7 18:34:33 2011
[0]PETSC ERROR: Configure options --with-shared --with-debugging=0 --useThreads 0 --with-clanguage=C++ --with-c-support --with-fortran-interfaces=1 --with-mpi-dir=/usr/lib/openmpi --with-mpi-shared=1 --with-blas-lib=-lblas --with-lapack-lib=-llapack --with-blacs=1 --with-blacs-include=/usr/include --with-blacs-lib="[/usr/lib/libblacsCinit-openmpi.so,/usr/lib/libblacs-openmpi.so]" --with-scalapack=1 --with-scalapack-include=/usr/include --with-scalapack-lib=/usr/lib/libscalapack-openmpi.so --with-mumps=1 --with-mumps-include=/usr/include --with-mumps-lib="[/usr/lib/libdmumps.so,/usr/lib/libzmumps.so,/usr/lib/libsmumps.so,/usr/lib/libcmumps.so,/usr/lib/libmumps_common.so,/usr/lib/libpord.so]" --with-umfpack=1 --with-umfpack-include=/usr/include/suitesparse --with-umfpack-lib="[/usr/lib/libumfpack.so,/usr/lib/libamd.so]" --with-spooles=1 --with-spooles-include=/usr/include/spooles --with-spooles-lib=/usr/lib/libspooles.so --with-hypre=1 --with-hypre-dir=/usr --with-scotch=1 --with-scotch-include=/usr/include/scotch --with-scotch-lib=/usr/lib/libscotch.so --with-hdf5=1 --with-hdf5-dir=/usr
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: MatGetOrdering() line 188 in src/mat/order/sorder.c
[0]PETSC ERROR: PCSetUp_LU() line 125 in src/ksp/pc/impls/factor/lu/lu.c
[0]PETSC ERROR: PCSetUp() line 795 in src/ksp/pc/interface/precon.c
[0]PETSC ERROR: KSPSetUp() line 237 in src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: KSPSolve() line 353 in src/ksp/ksp/interface/itfunc.c

I just call M.apply("insert"); before call "solve"

Thank you very much for your help. This solved my problem.