Petcs Matrix Multiplication and Transposition

Question

Given two matrices A and P I want to compute the PtAP product B = P^T A P.
To this aim, I have created a cpp extension to wrap the routine MatPtAP in PETSc.

Even though the results of the matrix multiplication is correct, something is still missing in the dolfin::PETScMatrix object. In particular, even though I can perform matvecs with B, the method B.init_vector(x,dim) does not seem to do its job correctly.
More specifically, if I try to call x.set_local( np.array(....) )
I get an error saying that "Local to global never set with VecSetLocalToGlobalMapping()".

I am attaching a basic example to reproduce this issue, so that I can get some feedback.

File "testLinalg/myla.h"

#include <dolfin/la/GenericMatrix.h>
#include <dolfin/la/Matrix.h>
#include <dolfin/la/PETScMatrix.h>

namespace dolfin
{
class cpp_linalg
{
public:
    Matrix MatPtAP(const GenericMatrix & A, const GenericMatrix & P);
};
}

File "testLinalg/myla.cpp

#include "myla.h"
namespace dolfin
{

 Matrix cpp_linalg::MatPtAP(const GenericMatrix & A, const GenericMatrix & P)
{
    const PETScMatrix* Ap = &as_type<const PETScMatrix>(A);
    const PETScMatrix* Pp = &as_type<const PETScMatrix>(P);
    Mat CC;\
    PetscErrorCode ierr = ::MatPtAP(Ap->mat(),Pp->mat(),MAT_INITIAL_MATRIX, 1.0,&CC);\
    PETScMatrix CCC = PETScMatrix(CC);\
    return Matrix(CCC);
  }

}

Python driver

import dolfin as dl
import numpy as np
import os

#Compile the cpp module:
abspath = os.path.dirname( os.path.abspath(__file__) )
sdir = os.path.join(abspath,"testLinalg")
header_file = open(os.path.join(sdir,"myla.h"), "r")
code = header_file.read()
header_file.close()
cpp_sources = ["myla.cpp"]  
cpp_module = dl.compile_extension_module(
code=code, source_directory=sdir, sources=cpp_sources,
include_dirs=[".", sdir])

#Generate a mesh and some matrices
nx = 32
ny = 32
mesh = dl.UnitSquareMesh(nx, ny)
Vh = dl.FunctionSpace(mesh, 'Lagrange', 1)

ndof = Vh.dim()

uh, vh = dl.TrialFunction(Vh), dl.TestFunction(Vh)

A = dl.assemble( dl.inner( dl.nabla_grad(uh), dl.nabla_grad(vh) ) *dl.dx )
M = dl.assemble( dl.inner( uh, vh ) *dl.dx )
ones = dl.interpolate(dl.Constant(1.), Vh).vector()
Mones = M*ones
s = Mones
s.set_local( np.ones(ndof) / Mones.array() )
M.zero()
M.set_diagonal(s)

myla = cpp_module.cpp_linalg()
B = myla.MatPtAP(M, A)

x = dl.Vector()
B.init_vector(x,1)

#This works fine
B.mult(s,x)

# Here I get an error
x.set_local(s.array())

In particular, here it is the error:

[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Object is in wrong state!
[0]PETSC ERROR: Local to global never set with VecSetLocalToGlobalMapping()!
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 3, Wed Aug 29 11:26:24 CDT 2012 
[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 darwin13. named dhcp-67-195.ices.utexas.edu by uvilla Tue May 26 16:33:30 2015
[0]PETSC ERROR: Libraries linked from /Users/johannr/fenics-1.5.0/local/lib
[0]PETSC ERROR: Configure run at Mon Jan 12 22:05:37 2015
[0]PETSC ERROR: Configure options --prefix=/Users/johannr/fenics-1.5.0/local COPTFLAGS=-O2 --with-debugging=0 --with-clanguage=cxx --with-c-support=1 --with-blas-lapack-dir=/usr --with-umfpack=1 --with-umfpack-include=/Users/johannr/fenics-1.5.0/local/include/suitesparse --with-umfpack-lib="[/Users/johannr/fenics-1.5.0/local/lib/libumfpack.a,/Users/johannr/fenics-1.5.0/local/lib/libamd.a]" --with-spooles=1 --with-spooles-include=/Users/johannr/fenics-1.5.0/local/include --with-spooles-lib=/Users/johannr/fenics-1.5.0/local/lib/libspooles.a --with-ptscotch=1 --with-ptscotch-dir=/Users/johannr/fenics-1.5.0/local --with-ml=1 --with-ml-include=/Users/johannr/fenics-1.5.0/local/include/trilinos --with-ml-lib=/Users/johannr/fenics-1.5.0/local/lib/libml.dylib --with-hdf5=1 --with-hdf5-dir=/Users/johannr/fenics-1.5.0/local --with-x=0 -with-x11=0 --with-fortran=0 --with-shared-libraries=1 PETSC_DIR=/Users/johannr/fenics-1.5.0/fenics-superbuild/build-fenics/CMakeExternals/src/PETSc PETSC_ARCH=darwin13.4.0-cxx-opt
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: VecSetValuesLocal() line 964 in /Users/johannr/fenics-1.5.0/fenics-superbuild/build-fenics/CMakeExternals/src/PETSc/src/vec/vec/interface/rvector.c
Traceback (most recent call last):
  File "/Users/uvilla/python_workspace/siamreviewed/codes/unit_test/testLinalg.py", line 41, in <module>
    x.set_local(s.array())

Any suggestion, it is greatly appreciated

Answer 1

Ok, I was able to find the solution.

Apparently in Petsc when matrix-matrix multiplication are performed the Local to Global mapping is not propagated from the inputs to the outputs. Therefore my solution was to use MatGetLocalToGlobalMapping/MatSetLocalToGlobalMapping to manually set the ISLocalToGlobalMapping.

Below, is how the code looks like:

Matrix cpp_linalg::MatPtAP(const GenericMatrix & A, const GenericMatrix & P)
{
    const PETScMatrix* Ap = &as_type<const PETScMatrix>(A);
    const PETScMatrix* Pp = &as_type<const PETScMatrix>(P);
    Mat CC;
    PetscErrorCode ierr = ::MatPtAP(Ap->mat(),Pp->mat(),MAT_INITIAL_MATRIX, 1.0,&CC);

    //Manually set the LocalToGlobalMapping
    ISLocalToGlobalMapping mapping;
    MatGetLocalToGlobalMapping(Pp->mat(),NULL, &mapping);
    MatSetLocalToGlobalMapping(CC, mapping, mapping);


    PETScMatrix CCC = PETScMatrix(CC);

    return Matrix(CCC);
}

Answer 2

Maybe this will help : http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/
in particular : http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMultTranspose.html

Comments

Hi Claire,

Thanks for your answer, but it does not really help my case. What I need is to build a NEW matrix B = P^t A P, and not just to be able to apply that operator to a vector.

I can create the matrix B in Fenics using the instant compiler, but the result is still missing some information regarding the Local to Global Mapping.

Best,

Umbe