# Creating Mat<complex> from Mat<real>

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

You canâ€™t. PETSc is written in C, to this date, it is not possible to mix both type of `PetscScalar`, see this 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.

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.

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

1 Like

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);
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);
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;
}``````

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).

1 Like

Good point! This probably only works simply for homogeneous BCâ€™s.

@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);``````

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.

``````mesh Th = square(2, 2);
include "macro_ddm.idp"
fespace Vh(Th, P1);
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.

``````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.

``````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`.

1 Like

Fantastic! Thank you

@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>
``````

``````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.

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.

1 Like

Thank you for the tip. Iâ€™ll give that a shot!