# Error with matrix addition in PETSc

**URL:** <https://community.freefem.org/t/error-with-matrix-addition-in-petsc/4200>\
**Category:** General Discussion\
**Created:** [January 23, 2026, 3:17pm UTC](https://community.freefem.org/t/error-with-matrix-addition-in-petsc/4200 "2026-01-23T15:17:11Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![mayeulchvne](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/mayeulchvne/32/3585_2.png) [@mayeulchvne](https://community.freefem.org/u/mayeulchvne)\
**Post date:** [January 23, 2026, 3:17pm UTC](https://community.freefem.org/t/error-with-matrix-addition-in-petsc/4200/1 "2026-01-23T15:17:11Z")

</div>

Dear Freefem community,

I have an issue with the addition of two matrices of type Mat. I read in the Forum that it should work since v4.9, however I receive the following error message :

Bmat = Amat + Tmat2; error operator + \<N5PETSc14DistributedCSRIN5HPDDM7SchwarzINSt3\_\_17complexIdEEEEEE\>, \<N5PETSc14DistributedCSRIN5HPDDM7SchwarzINSt3\_\_17complexIdEEEEEE\>  
List of choices

Can you help me understand ?

Best,

Mayeul

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [January 25, 2026, 4:17am UTC](https://community.freefem.org/t/error-with-matrix-addition-in-petsc/4200/2 "2026-01-25T04:17:35Z")

</div>

Without the code, no.

---

<div class="post-metadata">

**Author:** ![mayeulchvne](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/mayeulchvne/32/3585_2.png) [@mayeulchvne](https://community.freefem.org/u/mayeulchvne)\
**Post date:** [January 26, 2026, 1:38pm UTC](https://community.freefem.org/t/error-with-matrix-addition-in-petsc/4200/3 "2026-01-26T13:38:53Z")

</div>

Here is a relevant portion of my code, where I define the matrices I try to subtract, namely a variational formulation matrix and a spherical DtN matrix.

// ---------- DtN ----------

complex[int, int] SHvalues = SphericalHarmonics(maxDegree, xvals, yvals, zvals);

cout \<\< “Spherical Harmonics computed” \<\< endl;

int nbHarmonics = numberOfHarmonics(maxDegree);

int nbLegendre = numberOfLegendre(maxDegree);

complex[int,int] Earray(Vh.ndof, nbHarmonics);

complex[int] tmp(Vh.ndof);

matrix E;

VhSext sphericalHarmonic;

for (int l=0; l\<=maxDegree; l++) {

sphericalHarmonic[] = SHvalues(:, numberLegendre(l,0));

cout \<\< "Computing DtN coefficients for l = " \<\< l \<\< “, m = 0” \<\< endl;

varf FourierCoef(u,v) = int2d(Th,10)(v\*sphericalHarmonic);

Earray(:,numberHarmonic(l,0)) = FourierCoef(0,Vh);

for (int m=1; m\<=l; m++) {

```
sphericalHarmonic\[\] = SHvalues(:, numberLegendre(l,m));

cout << "Computing DtN coefficients for l = " << l << ", m = +/-" << m << endl;

```

varf FourierCoef(u,v) = int2d(Th,10)(v\*sphericalHarmonic);

```
tmp = FourierCoef(0,Vh);

Earray(:,numberHarmonic(l,-m)) = tmp;

Earray(:,numberHarmonic(l,m)) = conj(tmp);

```

}

}

E = Earray;

cout \<\< “Matrix of spherical harmonics projections defined” \<\< endl;

// Initialization of the diagonal matrix of the DtN coefficients

matrix D;

complex[int] Darray(nbHarmonics);

int n;

complex[int] valuesHankel = Hankel(maxDegree, k\*Rext);

for (int l = 0; l \<= maxDegree; l++) {

for (int m = -l; m \<= l; m++) {

```
n = numberHarmonic(l,m);

Darray\[n\] = (l/Rext - k\*valuesHankel\[l + 1\]/valuesHankel\[l\])/(Rext\*Rext);

```

}

}

D = [Darray];

// cout \<\< " D = " \<\< D \<\< endl;

cout \<\< “Diagonal matrix of DtN coefficients defined” \<\< endl;

// Final assembly

matrix T;

Mat Tmat;

T = E\*D;

Tmat = T;

cout \<\< “DtN operator defined” \<\< endl;

Mat Eprime;

matrix Estar = E’;

Eprime = Estar;

MatMatMult(Tmat,Eprime,Tmat);

cout \<\< “DtN operator computed” \<\< endl;

// ---------- Verification of the DtN: the fundamental solution Phi(z), on a hard sphere ----------

matrix A;

Mat Amat, Bmat;

// ---------- Variationnal formulation ----------

varf a(u,v) =

int3d(Th)( dx(u) \* dx(v)’ + dy(u) \* dy(v)’ + dz(u) \* dz(v)’ )

+ int3d(Th)(- k^2 \* u \* v’);

A = a(Vh, Vh);

Amat = A;

Bmat = Amat - Tmat;

cout \<\< “Matrix assembled” \<\< endl;

As I wrote above, the error I get is :

Bmat = Amat - Tmat; error operator - \<N5PETSc14DistributedCSRIN5HPDDM7SchwarzINSt3\_\_17complexIdEEEEEE\>, \<N5PETSc14DistributedCSRIN5HPDDM7SchwarzINSt3\_\_17complexIdEEEEEE\>

Thank you in advance for your help,

Mayeul

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [January 26, 2026, 2:22pm UTC](https://community.freefem.org/t/error-with-matrix-addition-in-petsc/4200/4 "2026-01-26T14:22:02Z")

</div>

Your code is unreadable, make sure it’s inside a code block. On top of that, you can see in [+= and \*= for PETSc Mat. · FreeFem/FreeFem-sources@63a93ac · GitHub](https://github.com/FreeFem/FreeFem-sources/commit/63a93acc8a3c978a0a97b9e495b9cd9cee9558cb) that there is support for `+=`, not `+`. You should probably just do operations with FreeFEM `matrix` and then feed the result into a `Mat`.

---

<div class="post-metadata">

**Author:** ![mayeulchvne](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/mayeulchvne/32/3585_2.png) [@mayeulchvne](https://community.freefem.org/u/mayeulchvne)\
**Post date:** [January 26, 2026, 2:40pm UTC](https://community.freefem.org/t/error-with-matrix-addition-in-petsc/4200/6 "2026-01-26T14:40:37Z")

</div>

(sorry I deleted my post by mistake)

Do you want to compile my code ? Because if so you might need the whole code with the definition of the functions that compute the Spherical harmonics and the spherical Hankel functions. If not, here is the code in a code box.

```cpp
// ---------- DtN ----------

complex[int, int] SHvalues = SphericalHarmonics(maxDegree, xvals, yvals, zvals);
cout << “Spherical Harmonics computed” << endl;

int nbHarmonics = numberOfHarmonics(maxDegree);
int nbLegendre = numberOfLegendre(maxDegree);
complex[int,int] Earray(Vh.ndof, nbHarmonics);
complex[int] tmp(Vh.ndof);
matrix<complex> E;
VhSext sphericalHarmonic;
for (int l=0; l<=maxDegree; l++) {
 sphericalHarmonic[] = SHvalues(:, numberLegendre(l,0));
 cout << "Computing DtN coefficients for l = " << l << “, m = 0” << endl;
 varf FourierCoef(u,v) = int2d(Th,10)(v*sphericalHarmonic);
 Earray(:,numberHarmonic(l,0)) = FourierCoef(0,Vh);
 for (int m=1; m<=l; m++) {
  sphericalHarmonic\[\] = SHvalues(:, numberLegendre(l,m));
  cout << "Computing DtN coefficients for l = " << l << ", m = +/-" << m << endl;
  varf FourierCoef(u,v) = int2d(Th,10)(v*sphericalHarmonic);
  tmp = FourierCoef(0,Vh);
  Earray(:,numberHarmonic(l,-m)) = tmp;
  Earray(:,numberHarmonic(l,m)) = conj(tmp);
  }
 }
E = Earray;
cout << “Matrix of spherical harmonics projections defined” << endl;

// Initialization of the diagonal matrix of the DtN coefficients
matrix<complex> D;
complex[int] Darray(nbHarmonics);
int n;
complex[int] valuesHankel = Hankel(maxDegree, k*Rext);
 for (int l = 0; l <= maxDegree; l++) {
  for (int m = -l; m <= l; m++) {
  n = numberHarmonic(l,m);
  Darray[n] = (l/Rext - k*valuesHankel[l + 1]/valuesHankel[l])/(Rext*Rext);
  }
 }
D = [Darray];

// cout << " D = " << D << endl;
cout << “Diagonal matrix of DtN coefficients defined” << endl;

// Final assembly
matrix<complex> T;
Mat<complex> Tmat;
T = E*D;
Tmat = T;
cout << “DtN operator defined” << endl;

Mat<complex> Eprime;
matrix<complex> Estar = E’;
Eprime = Estar;
MatMatMult(Tmat,Eprime,Tmat);
cout << “DtN operator computed” << endl;

// ---------- Verification of the DtN: the fundamental solution Phi(z), on a hard sphere ----------
matrix<complex> A;
Mat<complex> Amat, Bmat;

// ---------- Variationnal formulation ----------
varf a(u,v) =
int3d(Th)( dx(u) * dx(v)’ + dy(u) * dy(v)’ + dz(u) * dz(v)’ )
+ int3d(Th)(- k^2 * u * v’);
A = a(Vh, Vh);
Amat = A;
Bmat = Amat - Tmat;

```

Thank you for your help.

Sincerely,

Mayeul

---

<div class="post-metadata">

**Author:** ![mayeulchvne](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/mayeulchvne/32/3585_2.png) [@mayeulchvne](https://community.freefem.org/u/mayeulchvne)\
**Post date:** [January 26, 2026, 2:53pm UTC](https://community.freefem.org/t/error-with-matrix-addition-in-petsc/4200/7 "2026-01-26T14:53:06Z")

</div>

> [@prj](#):
>
> [+= and \*= for PETSc Mat. · FreeFem/FreeFem-sources@63a93ac · GitHub](https://github.com/FreeFem/FreeFem-sources/commit/63a93acc8a3c978a0a97b9e495b9cd9cee9558cb)

Thank you very much. Indeed, I tried to do the operations with Freefem `Matrix` but the product `MatMatMult(Tmat,Eprime,Tmat)` is becoming very long when the meshsize is decreasing, and it is significantly shorter in `Mat` format. If you know a way to cast back the result of this product from `Mat<complex>` to `matrix<complex>` it could solve my issues.

---

<div class="post-metadata">

**Author:** ![mayeulchvne](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/mayeulchvne/32/3585_2.png) [@mayeulchvne](https://community.freefem.org/u/mayeulchvne)\
**Post date:** [January 26, 2026, 3:29pm UTC](https://community.freefem.org/t/error-with-matrix-addition-in-petsc/4200/8 "2026-01-26T15:29:41Z")

</div>

I solved my issue with the `+=` operator and the `MatMatMult` function thank you for your help.