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