Creating Mat<complex> from Mat<real>

cmd · July 8, 2020, 7:37pm

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

prj · July 8, 2020, 7:44pm

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:

casting your real-valued problem into a complex-valued problem, sure there will be some overhead;
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.

cmd · July 8, 2020, 7:53pm

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.

prj · July 8, 2020, 7:58pm

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

cmd · July 9, 2020, 12:52pm

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

prj · July 9, 2020, 3:32pm

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

cmd · July 9, 2020, 4:15pm

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

cmd · July 10, 2020, 6:14pm

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

prj · July 11, 2020, 9:10am

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

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.

cmd · July 12, 2020, 5:51pm

Fantastic! Thank you

cmd · July 13, 2020, 12:42pm

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

prj · July 13, 2020, 12:46pm

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.

cmd · July 13, 2020, 12:48pm

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

Topic		Replies	Views
Complex matrix for EPSSolve General Discussion	4	207	January 17, 2024
X+a*Y with PETSc matrix (Mat) General Discussion	2	420	April 26, 2021
Complex eigenvalues with PETSc/SLEPc General Discussion	2	272	June 18, 2022
Computing numerical solution with parallel version General Discussion	3	303	March 4, 2022
Eigenmode issue for same matrix General Discussion	2	152	January 19, 2024

Creating Mat<complex> from Mat<real>

Related topics