dorbdb2 function

Implementation

void dorbdb2(
  final int M,
  final int P,
  final int Q,
  final Matrix<double> X11_,
  final int LDX11,
  final Matrix<double> X21_,
  final int LDX21,
  final Array<double> THETA_,
  final Array<double> PHI_,
  final Array<double> TAUP1_,
  final Array<double> TAUP2_,
  final Array<double> TAUQ1_,
  final Array<double> WORK_,
  final int LWORK,
  final Box<int> INFO,
) {
  final X11 = X11_.having(ld: LDX11);
  final X21 = X21_.having(ld: LDX21);
  final THETA = THETA_.having();
  final PHI = PHI_.having();
  final TAUP1 = TAUP1_.having();
  final TAUP2 = TAUP2_.having();
  final TAUQ1 = TAUQ1_.having();
  final WORK = WORK_.having();
  const NEGONE = -1.0, ONE = 1.0;
  double C = 0, S = 0;
  int I, ILARF = 0, IORBDB5 = 0, LLARF, LORBDB5 = 0, LWORKMIN, LWORKOPT;
  bool LQUERY;
  final CHILDINFO = Box(0);

  // Test input arguments

  INFO.value = 0;
  LQUERY = LWORK == -1;

  if (M < 0) {
    INFO.value = -1;
  } else if (P < 0 || P > M - P) {
    INFO.value = -2;
  } else if (Q < 0 || Q < P || M - Q < P) {
    INFO.value = -3;
  } else if (LDX11 < max(1, P)) {
    INFO.value = -5;
  } else if (LDX21 < max(1, M - P)) {
    INFO.value = -7;
  }

  // Compute workspace

  if (INFO.value == 0) {
    ILARF = 2;
    LLARF = max(P - 1, max(M - P, Q - 1));
    IORBDB5 = 2;
    LORBDB5 = Q - 1;
    LWORKOPT = max(ILARF + LLARF - 1, IORBDB5 + LORBDB5 - 1);
    LWORKMIN = LWORKOPT;
    WORK[1] = LWORKOPT.toDouble();
    if (LWORK < LWORKMIN && !LQUERY) {
      INFO.value = -14;
    }
  }
  if (INFO.value != 0) {
    xerbla('DORBDB2', -INFO.value);
    return;
  } else if (LQUERY) {
    return;
  }

  // Reduce rows 1, ..., P of X11 and X21

  for (I = 1; I <= P; I++) {
    if (I > 1) {
      drot(Q - I + 1, X11(I, I).asArray(), LDX11, X21(I - 1, I).asArray(),
          LDX21, C, S);
    }
    dlarfgp(
        Q - I + 1, X11.box(I, I), X11(I, I + 1).asArray(), LDX11, TAUQ1.box(I));
    C = X11[I][I];
    X11[I][I] = ONE;
    dlarf('R', P - I, Q - I + 1, X11(I, I).asArray(), LDX11, TAUQ1[I],
        X11(I + 1, I), LDX11, WORK(ILARF));
    dlarf('R', M - P - I + 1, Q - I + 1, X11(I, I).asArray(), LDX11, TAUQ1[I],
        X21(I, I), LDX21, WORK(ILARF));
    S = sqrt(pow(dnrm2(P - I, X11(I + 1, I).asArray(), 1), 2) +
        pow(dnrm2(M - P - I + 1, X21(I, I).asArray(), 1), 2));
    THETA[I] = atan2(S, C);

    dorbdb5(
        P - I,
        M - P - I + 1,
        Q - I,
        X11(I + 1, I).asArray(),
        1,
        X21(I, I).asArray(),
        1,
        X11(I + 1, I + 1),
        LDX11,
        X21(I, I + 1),
        LDX21,
        WORK(IORBDB5),
        LORBDB5,
        CHILDINFO);
    dscal(P - I, NEGONE, X11(I + 1, I).asArray(), 1);
    dlarfgp(
        M - P - I + 1, X21.box(I, I), X21(I + 1, I).asArray(), 1, TAUP2.box(I));
    if (I < P) {
      dlarfgp(
          P - I, X11.box(I + 1, I), X11(I + 2, I).asArray(), 1, TAUP1.box(I));
      PHI[I] = atan2(X11[I + 1][I], X21[I][I]).toDouble();
      C = cos(PHI[I]);
      S = sin(PHI[I]);
      X11[I + 1][I] = ONE;
      dlarf('L', P - I, Q - I, X11(I + 1, I).asArray(), 1, TAUP1[I],
          X11(I + 1, I + 1), LDX11, WORK(ILARF));
    }
    X21[I][I] = ONE;
    dlarf('L', M - P - I + 1, Q - I, X21(I, I).asArray(), 1, TAUP2[I],
        X21(I, I + 1), LDX21, WORK(ILARF));
  }

  // Reduce the bottom-right portion of X21 to the identity matrix

  for (I = P + 1; I <= Q; I++) {
    dlarfgp(
        M - P - I + 1, X21.box(I, I), X21(I + 1, I).asArray(), 1, TAUP2.box(I));
    X21[I][I] = ONE;
    dlarf('L', M - P - I + 1, Q - I, X21(I, I).asArray(), 1, TAUP2[I],
        X21(I, I + 1), LDX21, WORK(ILARF));
  }
}