<div dir="ltr"><div dir="ltr">
<div><div><div><div><div><br clear="all"></div>Hi Jed/PETSc-developers, <br><br></div>My goal is to invert a set of these PDE's to obtain a series of parameters F_t (with TSSolve and TSAdjoint for function/gradient computation). I was planning to use TAO for setting up the inverse problem but given that TAO doesn't support complex scalars, I'm re-thinking about converting this to a real formulation.<br><br></div>In <a href="https://doi.org/10.1007/s00466-006-0047-8" target="_blank">https://doi.org/10.1007/s00466-006-0047-8</a>, Mark Adams explains that the K1 formulation of Day/Heroux is well suited to block matrices in PETSc but the PDE described there has an large SPD operator arising out of the underlying elliptic PDE.
( Day/Heroux in their paper claim that K2/K3 formulations are problematic for Krylov solvers due to non ideal eigenvalue spectra created by conversion to real formulation).
<br><br></div>Now, the system I have is :
u_t = A*(u_xx + u_yy) + F_t*u; (A is purely imaginary and F_t is complex, with abs(A/F) ~ 1e-16. The parabolic PDE is converted to a series of TS solves each being elliptic). I implemented the K1/K4 approaches by using DMDA to manage the 2-dof grid instead of setting it up as one large vector of [real,imag] and those didn't converge well either (at least with simple preconditioners). <br><br>Any pointers as to what I could do to make the real formulation well conditioned ? Or should I not bother with this for now and implement a first order gradient descent method in PETSc (while approximating the regularizer as a cost integrand) ?<br></div><div><div><div><div><div><div><div><div><div><div><div><br><br><div>Thank You,<br></div><div><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div style="font-size:12.8px">Sajid Ali | PhD Candidate<br></div><div style="font-size:12.8px">Applied Physics<br></div><div style="font-size:12.8px">Northwestern University</div><div style="font-size:12.8px"><a href="http://s-sajid-ali.github.io" target="_blank">s-sajid-ali.github.io</a></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div>
</div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Mar 27, 2019 at 9:36 PM Jed Brown <<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">When you roll your own equivalent real formulation, PETSc has no way of<br>
knowing what conjugate transpose might mean, thus symmetry is lost. I<br>
would suggest just using the AVX2 implementation for now and putting in<br>
a request (or contributing a patch) for AVX-512 complex optimizations.<br>
<br>
Sajid Ali via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> writes:<br>
<br>
> Hi,<br>
><br>
> I'm able to solve the following equation using complex numbers (with<br>
> ts_type cn and pc_type gamg) :<br>
> u_t = A*u'' + F_t*u;<br>
> (where A = -1j/(2k) amd u'' refers to u_xx+u_yy implemented with the<br>
> familiar 5-point stencil)<br>
><br>
> Now, I want to solve the same problem using real numbers. The equivalent<br>
> equations are:<br>
> u_t_real = 1/(2k) * u''_imag + F_real*u_real - F_imag*u_imag<br>
> u_t_imag = -1/(2k) * u''_real + F_imag*u_real - F_real*u_imag<br>
><br>
> Thus, if we now take our new u vector to have twice the length of the<br>
> problem we're solving, keeping the first half as real and the second half<br>
> as imaginary, we'd get a matrix that had matrices computing the laplacian<br>
> via the 5-point stencil in the top-right and bottom-left corners and a<br>
> diagonal [F_real+F_imag, F_real-F_imag] term.<br>
><br>
> I tried doing this and the gamg preconditioner complains about an<br>
> unsymmetric matrix. If i use the default preconditioner, I get<br>
> DIVERGED_NONLINEAR_SOLVE.<br>
><br>
> Is there a way to better organize the matrix ?<br>
><br>
> PS: I'm trying to do this using only real numbers because I realized that<br>
> the optimized avx-512 kernels for KNL are not implemented for complex<br>
> numbers. Would that be implemented soon ?<br>
><br>
> Thank You,<br>
> Sajid Ali<br>
> Applied Physics<br>
> Northwestern University<br>
</blockquote></div><br clear="all"><br>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div style="font-size:12.8px">Sajid Ali | PhD Candidate<br></div><div style="font-size:12.8px">Applied Physics<br></div><div style="font-size:12.8px">Northwestern University</div><div style="font-size:12.8px"><a href="http://s-sajid-ali.github.io" target="_blank">s-sajid-ali.github.io</a></div></div></div></div></div></div></div></div>