[petsc-dev] various dot products messed up for complex numbers?

[Barry Smith]

On Apr 17, 2012, at 9:57 PM, Jed Brown wrote:

> Where is the update to every KSP (at least the CGs, but I think all of them) to change the order?

   The assumption is that all the KSPs were coded assuming the complex conjugation of the second argument.  If they were actually coded using the complex conjugate of the first argument then they have been WRONG for many years because PETSc has always implemented the inner products complex conjugating the second argument. (PETSc never used the zdotc() routine in BLAS because of the returning a complex double business so PETSc always had handcoded dot() within a #if def(complex) using the complex conjugate of the second argument.

   Since I do not know if they have always been coded wrong (and am assuming they are coded correctly, at least cg) changing them all by brute force now would make them all wrong now.  Certainly if some are wrong now (which means they have always been wrong) they should be fixed, but that needs to be done by inspection/testing to determine which ones are wrong.

   The truth is likely that except for CG most of the methods were coded without regard to complex and are likely wrong but not necessarily wrong in a consistent way or even self-consistently wrong.


    Barry


> I don't see it in the repository. Note that this notation is different from that used by most textbook and paper descriptions of the methods.
> 
> On Apr 17, 2012 7:09 PM, "Barry Smith" <bsmith at mcs.anl.gov> wrote:
> >>
> >>
> >>   Since PETSc was written and always implemented with mathematician style* inner products, I have decided to keep it that way and call the BLASdot_() with the arguments reversed. Note that until Matt added the support for zdotc() we never called the BLAS for complex inner products and thus never used the engineers style inner products.
> >>
> >>    Please report any problems,
> >>
> >>    Barry
> >>
> >> * please do not be upset by my joking reference to mathematician and engineering style inner products.
> >>
> >>
> >> On Apr 15, 2012, at 1:09 PM, Jed Brown wrote:
> >>
> >> > On Sun, Apr 15, 2012 at 12:55, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >> >   Real mathematicians conjugate the second one :-)
> >> >
> >> > I never unrderstood that choice. I like vectors to be column vectors, one forms to be row vectors, and operations to have standard fixity. How some mathematicians ended up with infix operators and postfix inner products is beyond me.
> >> >
> >> >
> >> >   The code to check is PETSc's conjugate gradient, there the conjugate does mater :-).
> >> >
> >> > Looks right to me.
> >> >
> >> >        ierr = KSP_MatMult(ksp,Amat,P,W);CHKERRQ(ierr);          /*     w <- Ap         */
> >> >        ierr = VecXDot(P,W,&dpi);CHKERRQ(ierr);                  /*     dpi <- p'w     */
> >> >
> >>