PETSc operators

How to perform simple operations in PETSc like matrix addition, multiplication, constant times a matrix? I’ve tried MatMatMult and MatScale both did not work.

MatMatMult should work if you have a recent enough FreeFEM installation, cf. https://github.com/FreeFem/FreeFem-sources/blob/develop/plugin/mpi/PETSc-code.hpp#L4460. MatScale is simply A *= alpha. Addition is A += alpha * B where both A and B are Mat. You can also do, assuming your Mat follow the same distributions:

matrix a = alpha * b + gamma * c; // sequential FreeFEM matrix
Mat B(A, a); // B and A follow the same distribution

This is usually more efficient since the sequential matrix operations do not require any communication.

1 Like

Is there a direct way to obtain the number of rows and columns like in sequential? ex: B.n or B.m

.n will give you the number of local rows of a Mat. .m does not exist, but since you need to Mat do define a rectangular one, i.e., Mat C(A, B), instead of doing C.m (which does not exist), you can do B.n.

Ah I see, thank you for your reply. I would like to add something about my first comment regarding matrix operations, whenever I try to do for instance A+=alpha I get the following error:
any idea on how to proceed?

[0]PETSC ERROR: ------------------------------------------------------------------------

[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range

[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger

[0]PETSC ERROR: or see https://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind

[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors

[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run

[0]PETSC ERROR: to get more information on the crash.

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------

[0]PETSC ERROR: Signal received

[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.

[0]PETSC ERROR: Petsc Release Version 3.15.0, Mar 30, 2021

[0]PETSC ERROR: /softs/local/freefem/4.9/bin/FreeFem+±mpi on a named node001.cluster by mnoun Mon Sep 6 10:39:22 2021

[0]PETSC ERROR: Configure options --prefix=/softs/local/freefem/4.9/ff-petsc/r MAKEFLAGS= --with-debugging=0 COPTFLAGS="-O3 -mtune=generic" CXXOPTFLAGS="-O3 -mtune=generic" FOPTFLAGS="-O3 -mtune=generic" --with-cxx-dialect=C++11 --with-ssl=0 --with-x=0 --with-fortran-bindings=0 --with-cudac=0 --with-scalar-type=real --with-blaslapack-include= --with-blaslapack-lib="-llapack -lblas" --download-metis --download-ptscotch --download-hypre --download-parmetis --download-mmg --download-parmmg --download-superlu --download-suitesparse --download-tetgen --download-slepc --download-hpddm --with-cc=gcc --with-cxx=g++ --download-mpich --with-fc=gfortran --download-scalapack --download-mumps --download-slepc-configure-arguments=–download-arpack=https://github.com/prj-/arpack-ng/archive/b64dccb.tar.gz PETSC_ARCH=fr

[0]PETSC ERROR: #1 User provided function() at unknown file:0

[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.

That’s a very generic error. It’s nearly impossible to help without a MWE to reproduce the error.

This is a part of the code : testpetsc.edp (3.0 KB) I defined three matrices K, C and M and I am trying to combine them in the matrix a1.

Neither a1 nor a2 are defined, so it doesn’t make sense to do a1 += something with an undefined a1. Just put everything in a matrix linearComb, and then do Mat a1(K, linearComb).

I defined a1 and a2 on line 81 as Mat a1, a2. Could you elaborate a little bit more your method please?
Edit: K, C and M are Mat type also

Sorry, I meant that none of the matrices are initialized, so you cannot do unintializedMat += something, it’s undefined behavior (and results in a segfault).

Instead of

K = stiffness(Vh, Vh);
C = damping(Vh,Vh);
M = mass(Vh,Vh);

Do

matrix locK = stiffness(Vh, Vh);
K = locK;
matrix locC = damping(Vh,Vh);
C = locC;
matrix locM = damping(Vh,Vh);
M = locM;
matrix locLinearComb = alpha * locK + beta * locC + gamma * locM;
Mat linearComb(M, locLinearComb);