415 Assert(J.TMV_rowsize() == 5);
416 Assert(J.TMV_colsize() == 5);
418 std::complex<double> zc(
x[0],
x[1]);
419 std::complex<double> g(
x[2],
x[3]);
440 J.col(0) =
dbdE.TMV_subVector(1,6) -
dbdE(0) * f;
445 J.col(1) =
dbdE.TMV_subVector(1,6) -
dbdE(0) * f;
466 J.col(2) =
dbdE.TMV_subVector(1,6) -
dbdE(0) * f;
473 J.col(3) =
dbdE.TMV_subVector(1,6) -
dbdE(0) * f;
490 J.col(4) =
dbdE.TMV_subVector(1,6) -
dbdE(0) * f;
497 if (!
zerob11) J.row(4).setZero();
502 dbg<<
"dE = "<<dE<<std::endl;
523 DMatrix dDdmu_num = (D2-D1)/(2.*dE);
524 DMatrix d2D_num = (D2+D1-2.*D0)/(dE*dE);
526 dbg<<
"Gmu = "<<Gmubig.TMV_subMatrix(0,6,0,6)<<std::endl;
527 dbg<<
"D = "<<Dbig.TMV_subMatrix(0,6,0,6)<<std::endl;
528 dbg<<
"dDdmu = "<<dDdmu.TMV_subMatrix(0,6,0,6)<<std::endl;
529 dbg<<
"dDdmu_num = "<<dDdmu_num.TMV_subMatrix(0,6,0,6)<<std::endl;
530 dbg<<
"Norm(dD/dmu - numeric dD/dmu) = "<<Norm(dDdmu-dDdmu_num)<<std::endl;
531 dbg<<
"dE*Norm(d2D_num) = "<<dE*Norm(d2D_num)<<std::endl;
532 test(Norm(dDdmu-dDdmu_num) < 30.*dE*Norm(d2D_num),
"dDdmu");
550 DMatrix dSdg1_num = (S2-S1)/(2.*dE);
551 DMatrix d2S_num = (S2+S1-2.*S0)/(dE*dE);
553 dbg<<
"Gg1 = "<<Gg1big.TMV_subMatrix(0,6,0,6)<<std::endl;
554 dbg<<
"Gth = "<<Gthbig.TMV_subMatrix(0,6,0,6)<<std::endl;
555 dbg<<
"S = "<<Sbig.TMV_subMatrix(0,6,0,6)<<std::endl;
556 dbg<<
"dSdg1 = "<<dSdg1.TMV_subMatrix(0,6,0,6)<<std::endl;
557 dbg<<
"dSdg1_num = "<<dSdg1_num.TMV_subMatrix(0,6,0,6)<<std::endl;
558 dbg<<
"Norm(dS/dg1 - numeric dS/dg1) = "<<Norm(dSdg1-dSdg1_num)<<std::endl;
559 dbg<<
"dE*Norm(d2S_num) = "<<dE*Norm(d2S_num)<<std::endl;
560 test(Norm(dSdg1-dSdg1_num) < 30.*dE*Norm(d2S_num),
"dSdg1");
570 DMatrix dSdg2_num = (S2-S1)/(2.*dE);
571 d2S_num = (S2+S1-2.*S0)/(dE*dE);
573 dbg<<
"Gg2 = "<<Gg2big.TMV_subMatrix(0,6,0,6)<<std::endl;
574 dbg<<
"Gth = "<<Gthbig.TMV_subMatrix(0,6,0,6)<<std::endl;
575 dbg<<
"S = "<<Sbig.TMV_subMatrix(0,6,0,6)<<std::endl;
576 dbg<<
"dSdg2 = "<<dSdg2.TMV_subMatrix(0,6,0,6)<<std::endl;
577 dbg<<
"dSdg2_num = "<<dSdg2_num.TMV_subMatrix(0,6,0,6)<<std::endl;
578 dbg<<
"Norm(dS/dg2 - numeric dS/dg2) = "<<Norm(dSdg2-dSdg2_num)<<std::endl;
579 dbg<<
"dE*Norm(d2S_num) = "<<dE*Norm(d2S_num)<<std::endl;
580 test(Norm(dSdg2-dSdg2_num) < 30.*dE*Norm(d2S_num),
"dSdg2");
596 DMatrix dTdx_num = (T2-T1)/(2.*dE);
597 DMatrix d2T_num = (T2+T1-2.*T0)/(dE*dE);
599 dbg<<
"Gx = "<<Gxbig.TMV_subMatrix(0,6,0,6)<<std::endl;
600 dbg<<
"T = "<<Tbig.TMV_subMatrix(0,6,0,6)<<std::endl;
601 dbg<<
"dTdx = "<<dTdx.TMV_subMatrix(0,6,0,6)<<std::endl;
602 dbg<<
"dTdx_num = "<<dTdx_num.TMV_subMatrix(0,6,0,6)<<std::endl;
603 dbg<<
"Norm(dT/dx - numeric dT/dx) = "<<Norm(dTdx-dTdx_num)<<std::endl;
604 dbg<<
"dE*Norm(d2T_num) = "<<dE*Norm(d2T_num)<<std::endl;
605 test(Norm(dTdx-dTdx_num) < 30.*dE*Norm(d2T_num),
"dTdx");
615 DMatrix dTdy_num = (T2-T1)/(2.*dE);
616 d2T_num = (T2+T1-2.*T0)/(dE*dE);
618 dbg<<
"Gy = "<<Gybig.TMV_subMatrix(0,6,0,6)<<std::endl;
619 dbg<<
"T = "<<Tbig.TMV_subMatrix(0,6,0,6)<<std::endl;
620 dbg<<
"dTdy = "<<dTdy.TMV_subMatrix(0,6,0,6)<<std::endl;
621 dbg<<
"dTdy_num = "<<dTdy_num.TMV_subMatrix(0,6,0,6)<<std::endl;
622 dbg<<
"Norm(dT/dy - numeric dT/dy) = "<<Norm(dTdy-dTdy_num)<<std::endl;
623 dbg<<
"dE*Norm(d2T_num) = "<<dE*Norm(d2T_num)<<std::endl;
624 test(Norm(dTdy-dTdy_num) < 30.*dE*Norm(d2T_num),
"dTdy");
635 dbg<<
"x0 = "<<x0<<std::endl;
636 dbg<<
"b0 = "<<bx0<<std::endl;
637 dbg<<
"f0 = "<<f0<<std::endl;
639 for(
int k=0;k<5;k++) {
640 dbg<<
"k = "<<k<<std::endl;
646 dbg<<
"x1 = "<<x1<<std::endl;
647 dbg<<
"b1 = "<<b1<<std::endl;
648 dbg<<
"f1 = "<<f1<<std::endl;
654 dbg<<
"x2 = "<<x2<<std::endl;
655 dbg<<
"b2 = "<<b2<<std::endl;
656 dbg<<
"f2 = "<<f2<<std::endl;
658 DVector dbdE_num = (b2-b1)/(2.*dE);
659 dbg<<
"dbdE_num = "<<dbdE_num<<std::endl;
660 DVector d2bdE = (b2+b1-2.*bx0)/(dE*dE);
661 dbg<<
"d2bdE = "<<d2bdE<<std::endl;
667 dbg<<
"Db0 = "<<Db0<<std::endl;
669 case 0: SDb0 =
S *
Db0;
670 dbg<<
"SDb0 = "<<SDb0<<std::endl;
672 dbg<<
"GDb0 = "<<GDb0<<std::endl;
674 dbg<<
"dbdx = "<<
dbdE<<std::endl;
675 dbg<<
"dbdx = dTdx S D b0 = "<<
678 case 1: SDb0 =
S *
Db0;
679 dbg<<
"SDb0 = "<<SDb0<<std::endl;
681 dbg<<
"GDb0 = "<<GDb0<<std::endl;
683 dbg<<
"dbdy = "<<
dbdE<<std::endl;
684 dbg<<
"dbdy = dTdy S D b0 = "<<
688 dbg<<
"GDb0 = "<<GDb0<<std::endl;
690 dbg<<
"GthDb0 = "<<GthDb0<<std::endl;
692 dbg<<
"GDb0 = "<<GDb0<<std::endl;
694 dbg<<
"SDb0 = "<<SDb0<<std::endl;
696 dbg<<
"dbdg1 = "<<
dbdE<<std::endl;
697 dbg<<
"dbdg1 = T dSdg1 D b0 = "<<
701 dbg<<
"GDb0 = "<<GDb0<<std::endl;
703 dbg<<
"GthDb0 = "<<GthDb0<<std::endl;
705 dbg<<
"GDb0 = "<<GDb0<<std::endl;
707 dbg<<
"SDb0 = "<<SDb0<<std::endl;
709 dbg<<
"dbdg2 = "<<
dbdE<<std::endl;
710 dbg<<
"dbdg2 = T dSdg2 D b0 = "<<
714 dbg<<
"GDb0 = "<<GDb0<<std::endl;
716 dbg<<
"Db0 = "<<Db0<<std::endl;
718 dbg<<
"SDb0 = "<<SDb0<<std::endl;
720 dbg<<
"dbdmu = "<<
dbdE<<std::endl;
721 dbg<<
"dbdmu = T S dDdmu b0 = "<<
726 dbg<<
"dbdE = "<<
dbdE<<std::endl;
727 dbg<<
"Norm(dbdE-dbdE_num) = "<<Norm(
dbdE-dbdE_num)<<std::endl;
728 dbg<<
"dE*Norm(d2bdE) = "<<dE*Norm(d2bdE)<<std::endl;
729 test(Norm(
dbdE-dbdE_num) < 30.*dE*Norm(d2bdE),
"dbdE");
733 dbg<<
"dfdE from dbdE = "<<dfdE<<std::endl;
734 dbg<<
"dfdE from dbdE_num = "<<
735 (dbdE_num.TMV_subVector(1,6) - dbdE_num(0) * f0)/
bx(0)<<std::endl;;
736 dbg<<
"dfdE in J = "<<J.col(k)<<std::endl;
737 J_num.col(k) = (f2-f1)/(2.*dE);
738 dbg<<
"dfdE_num = "<<J_num.col(k)<<std::endl;
739 DVector d2f_num = (f2+f1-2.*f0)/(dE*dE);
740 dbg<<
"d2f_num = "<<d2f_num<<std::endl;
741 dbg<<
"Norm(dfdE_num-dfdE) ="<<Norm(J_num.col(k)-dfdE)<<std::endl;
742 dbg<<
"dE*Norm(d2f_num) ="<<dE*Norm(d2f_num)<<std::endl;
743 test(Norm(J_num.col(k)-dfdE) < 30.*dE*Norm(d2f_num),
"dfdE");
746 dbg<<
"J_num = "<<J_num<<std::endl;
747 dbg<<
"analytic: "<<J<<std::endl;
750 if (
fixmu) J_num.col(4).setZero();
751 dbg<<
"J_num => "<<J_num<<std::endl;
752 dbg<<
"Norm(diff) = "<<Norm(J-J_num)<<std::endl;
753 dbg<<
"dE*Norm(J) = "<<dE*Norm(J_num)<<std::endl;
754 test(Norm(J_num-J) < dE*Norm(J_num),
"dfdE");
void augmentMuTransformCols(double mu, int order, DMatrix &D)
void setupGg1(DMatrix &Gg1, int order1, int order2)
void setupGth(DMatrix &Gth, int order1, int order2)
void calculateGTransform(std::complex< double > g, int order1, int order2, DMatrix &S)
void augmentZTransformCols(std::complex< double > z, int order, DMatrix &T)
void calculateZTransform(std::complex< double > z, int order1, int order2, DMatrix &T)
#define TMV_colRange(m, j1, j2)
void setupGx(DMatrix &Gx, int order1, int order2)
void setupGy(DMatrix &Gy, int order1, int order2)
void augmentGTransformCols(std::complex< double > g, int order, DMatrix &S)
void setupGg2(DMatrix &Gg2, int order1, int order2)
#define TMV_rowRange(m, i1, i2)
std::ostream *const dbgout
void doF3(const DVector &x, DVector &f) const
void calculateMuTransform(double mu, int order1, int order2, DMatrix &D)
void setupGmu(DMatrix &Gmu, int order1, int order2)