dolfin team mailing list archive

[Patricia Brunner <question149003@xxxxxxxxxxxxxxxxxxxxx>]

Question #149003 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/149003

    Status: Needs information => Open

Patricia Brunner gave more information on the question:
Hi,

I already found out how to implement this in principle (using a class and redefining mult). Nevertheless it does not work until now.
I have to solve a system with PETScKrylovSolver (in the mult part of my class) and I  get the error:
TypeError: unsupported operand type(s) for *: 'Matrix' and 'SwigPyObject'

I don't know if the problem could be solved using an already transposed
matrix instead of the class TransposeMatrix. I'd like to try this, but I
don't know how to set the AT-Matrix-entities to the ones of At =
transpose(A)?

Thanks for your help,
Patricia

I'm attaching my code below.

***************************
from dolfin import *
import numpy
from numpy import array
from numpy import transpose


class TransposeMatrix(PETScKrylovMatrix) : 
    def __init__(self, A) :
        PETScKrylovMatrix.__init__(self, A.size(0), A.size(1));
        self.A = A;

    def mult(self, x, b) :
        print "TransposeMatrix.mult";
        print self.A;
        #self.A.transpmult(*args);
        #self.A.mult(*args);
        b = self.A * x;

# Implements @(p) omCon .* A^{-T}(Mo * p)
class SaMatrix(PETScKrylovMatrix) :
    def __init__(self, A, Mo, V) :
        PETScKrylovMatrix.__init__(self, A.size(0), A.size(1));
        self.At = TransposeMatrix(A);
        self.Mo = Mo;
        self.solver = PETScKrylovSolver("gmres", "none");
        self.V = V;
        self.yt = Function(V);
        self.yt_petsc = PETScVector(self.yt.vector().size());

    def mult(self, *args) :
        Mop = self.Mo * args[0];
        self.yt_petsc.zero();
        self.solver.solve(self.At, self.yt_petsc, Mop);
        yf = project(self.yt, self.V);
        args[1].set_local(yf.vector().array());


parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

mesh = UnitSquare(64,64)
V0 = FunctionSpace(mesh, "DG", 0)
V2 = FunctionSpace(mesh, "CG", 2)

w = TrialFunction(V2)
v = TestFunction(V2)

z = Function(V2)
z.vector()[:] = 1.0

p = Function(V2)

a0 = (dot(grad(w),grad(v)) + w*v)*dx + 0.25*(w*v)*ds
f0 = v*dx

Mo = assemble(v*w*dx)

A = assemble(a0)
At = transpose(A)

AT = Matrix()

Sa = SaMatrix(A, Mo, V2)
y = PETScVector(z.vector().size());
Sa.mult(z.vector() - y, p);

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.