# Creating Mat\<complex\> from Mat\<real\>

**URL:** <https://community.freefem.org/t/creating-mat-complex-from-mat-real/502>\
**Category:** General Discussion\
**Created:** [July 8, 2020, 7:37pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502 "2020-07-08T19:37:58Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [July 8, 2020, 7:37pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/1 "2020-07-08T19:37:58Z")

</div>

Hopefully this is an easy question!

I have a problem involving both real and complex arithmetic. Once I have done:

```
Mat<real> Areal;
Mat<complex> Acplx;
buildDmesh(Th);
createMat(Th,Areal,Pk);

```

How can I then apply the same decomposition to `Acplx`?

Best,  
Chris

---

<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:** [July 8, 2020, 7:44pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/2 "2020-07-08T19:44:34Z")

</div>

You can’t. PETSc is written in C, to this date, it is not possible to mix both type of `PetscScalar`, see [this](https://gitlab.com/petsc/petsc/-/issues/32) 7-year old issue. What do you need to do with both `Areal` and `Acplx`?  
Some alternatives are:

1. casting your real-valued problem into a complex-valued problem, sure there will be some overhead;
2. decouple the real-valued part of your script from the complex-valued part of your script. E.g., if you do LSA, first compute your baseflow with `load "PETSc"`, then do the eigensolve in another script with `load "SLEPc-complex"`. You can use `loadDmesh` and `saveDmesh` to save everything PETSc-related.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [July 8, 2020, 7:53pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/3 "2020-07-08T19:53:53Z")

</div>

Bummer! I am solving the Moore-Spence system for a Hopf point where the real baseflow problem and complex eigenvalue problem are coupled. I will try using a 2x2 MatNest to get the real and imaginary parts of the complex problem in real arithmetic.

---

<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:** [July 8, 2020, 7:58pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/4 "2020-07-08T19:58:35Z")

</div>

You can do that as well. But still, you may want to make sure everything goes according to plan with everything in complex arithmetic first (with a 0 imaginary part for your baseflow). Probably easier to debug than a 2x2 MatNest solver, IMHO 😀

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [July 9, 2020, 12:52pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/5 "2020-07-09T12:52:51Z")

</div>

In case anyone is facing a similar problem in the future, this issue is solved pretty painlessly using the MatNest approach. If you wanted to solve a complex problem of the form:

```
Mat<complex> J;
Vh<complex> def(u);
func complex[int] funcRHS(complex[int]& inPETSc) { // func to assemble RHS
  changeNumbering(J,u[], inPETSc, inverse = true, exchange = true);
  complex[int] out(Vh.ndof);
  out = vR(0, Vh, tgv = -1); // evaluate RHS
  changeNumbering(J, out, inPETSc);
  return inPETSc;
}
func int funcLHS(complex[int]& inPETSc) { // func to assemble complex Mat
  changeNumbering(J, u[], inPETSc, inverse = true, exchange = true);
  J = vJ(Vh, Vh, tgv = -1); // evaluate LHS
  return 0;
}

```

Without changing your `varf`, you can simply switch to `Mat<real>` and solve a real problem of twice the size for the real and imaginary parts of the complex solution:

```
Mat<real> J;
Mat<real> J2;
Vh<complex> def(u);
func real[int] funcRHS(real[int]& inPETSc) { // func to assemble RHS
  // NOTE: the following line does not work!
  // changeNumbering([J,J], [u[].re,u[].im], inPETSc, inverse = true, exchange = true);
  { // Instead use:
   real[int] ur(Vh.ndof),ui(Vh.ndof);
   changeNumbering([J, J], [ur, ui], inPETSc, inverse = true, exchange = true);
   u[].re = ur;
   u[].im = ui;
  }
  complex[int] out(Vh.ndof);
  out = vR(0, Vh, tgv = -1); // evaluate RHS
  real[int] outPETSc;
  changeNumbering([J, J], [out.re, out.im], outPETSc);
  return outPETSc;
}
func int funcLHS(real[int]& inPETSc) {
  // NOTE: the following line does not work!
  // changeNumbering([J,J], [u[].re,u[].im], inPETSc, inverse = true, exchange = true);
  { // Instead use:
   real[int] ur(Vh.ndof),ui(Vh.ndof);
   changeNumbering([J, J], [ur, ui], inPETSc, inverse = true, exchange = true);
   u[].re = ur;
   u[].im = ui;
  }
  matrix<complex> mJ = vJ(Vh, Vh, tgv = -1); // evaluate LHS
  matrix mJ11 = mJ.re;
  matrix mJ21 = mJ.im;
  matrix mJ12 = -mJ21;
  Mat J11(J,mJ11);
  Mat J12(J,mJ12);
  Mat J21(J,mJ21);
  Mat Jblock = [[J11,J12],[J21,J11]];
  MatConvert(Jblock,J2);
  return 0;
}
```

---

<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:** [July 9, 2020, 3:32pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/6 "2020-07-09T15:32:04Z")

</div>

Careful about those boundary conditions in mJ21! 🙂  
You may have to do some postprocessing or `tgv = -10` magic (if you have a `varf` with only the imaginary part that you can then assemble into mJ21).

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [July 9, 2020, 4:15pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/7 "2020-07-09T16:15:52Z")

</div>

Good point! This probably only works simply for homogeneous BC’s.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [July 10, 2020, 6:14pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/8 "2020-07-10T18:14:08Z")

</div>

@prj I’m not sure this should count as a bug, but (per my edit above) I found some nuance in the way `changeNumbering` casts to the real and imag parts of complex arrays. In particular,

```
changeNumbering([J,J], [u[].re,u[].im], inPETSc, inverse = true, exchange = true);

```

does not update the complex fe array `u`. To produce the desired result, I did:

```
real[int] ur(Vh.ndof),ui(Vh.ndof);
changeNumbering([J, J], [ur, ui], inPETSc, inverse = true, exchange = true);
u[].re = ur;
u[].im = ui;

```

Note this was not necessary in the other direction, where the following lines achieved the desired effect.

```
real[int] u2;
changeNumbering([J,J], [u[].re,u[].im], u2);
```

---

<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:** [July 11, 2020, 9:10am UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/9 "2020-07-11T09:10:39Z")

</div>

Not really a bug per say, but more like an undocumented feature 😉  
It was not possible to pass “views” of vectors, e.g., `u[].re` is a view of the real part of `u[]`, `u[](0:20)` is another kind of a view, but `u[]` is the real deal.  
Anyway, ask, and it shall be given [you](https://github.com/FreeFem/FreeFem-sources/commit/26df3ea59be636ab1305ad8b6472d1fa89cfc502).

```auto
mesh Th = square(2, 2);
include "macro_ddm.idp"
fespace Vh(Th, P1);
load "PETSc"
Mat A;
createMat(Th, A, P1)
Vh<complex> u;
real[int] in(2 * A.n);
in(0:A.n-1) = 1;
in(A.n:2*A.n-1) = 10;
u[] = 0;
changeNumbering([A, A], [u[].re, u[].im], in, inverse = true, exchange = true);
if(mpirank == 0)
    cout << u[] << endl;
real[int] out(2 * Vh.ndof);
out = 0;
changeNumbering([A, A], [out(0:Vh.ndof-1), out(Vh.ndof:2*Vh.ndof-1)], in, inverse = true, exchange = true);
if(mpirank == 0)
    cout << out << endl;
out = 0;
changeNumbering([A, A], [out(0:2*Vh.ndof:2), out(1:2*Vh.ndof:2)], in, inverse = true, exchange = true);
if(mpirank == 0)
    cout << out << endl;

```

With the previous version of the plugin, you would get, e.g., with two processes.

```auto
8
	(0,0)	(0,0)	(0,0)	(0,0)	(0,0)
	(0,0)	(0,0)	(0,0)
16
	  0 0 0 0 0
	  0 0 0 0 0
	  0 0 0 0 0
	  0
16
	  0 0 0 0 0
	  0 0 0 0 0
	  0 0 0 0 0
	  0

```

And now.

```auto
8
	(1,10)	(1,10)	(1,10)	(1,10)	(1,10)
	(1,10)	(1,10)	(1,10)
16
	  1 1 1 1 1
	  1 1 1 10 10
	 10 10 10 10 10
	 10
16
	  1 10 1 10 1
	 10 1 10 1 10
	  1 10 1 10 1
	 10

```

You just need to recompile the plugin, so `git pull && cd plugin/mpi && ../seq/ff-c++ -auto PETSc.cpp`.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [July 12, 2020, 5:51pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/10 "2020-07-12T17:51:34Z")

</div>

Fantastic! Thank you 😃

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [July 13, 2020, 12:42pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/11 "2020-07-13T12:42:38Z")

</div>

@prj I have not been able to make this fix work yet (though my old way still does). I wonder if it might be related to something I did in my compilation process. Also, as an aside, I had to modify your example code to get it to work. First, I inserted at line 2:

```
mesh Th = square(2, 2);
macro dimension()2//
include "macro_ddm.idp"

```

to get the macros to work, but I still got the following error in the line `real[int] in(2*A.n);`:

```
No operator .n for type <N5PETSc14DistributedCSRIN5HPDDM7SchwarzIdEEEE>

```

So I am instead doing:

```
real[int] in;
changeNumbering([A, A], [u[].re, u[].im], in);
in(0:in.n/2-1) = 1;
in(in.n/2:in.n-1) = 10;

```

Having worked around these things, I have done `git pull` and recompiled using `ff-c++ -auto PETSc.cpp` in both the `./plugin/mpi/` and `./plugin/seq/` directories. However, my output still looks like what you got with the “previous version of the plugin.” I am running FreeFEM version 4.6 on Ubuntu. Does any potential cause for this jump out at you? If not I can always keep doing it using the fix I showed earlier.

---

<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:** [July 13, 2020, 12:46pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/12 "2020-07-13T12:46:44Z")

</div>

You probably made a `make install` with 4.6 which causes the latest version of `PETSc.so` to not be used.  
That’s why I highly recommend not to `make install`, as it just pollutes your system with FreeFEM files, but nothing else. You need to figure out where is the faulty old `.so`, or alternatively, fix your `~/.freefem++.pref` to load _first_ the `.so` from your custom build.

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [July 13, 2020, 12:48pm UTC](https://community.freefem.org/t/creating-mat-complex-from-mat-real/502/13 "2020-07-13T12:48:24Z")

</div>

Thank you for the tip. I’ll give that a shot!