dggsvd3 function

Implementation

void dggsvd3(
  final String JOBU,
  final String JOBV,
  final String JOBQ,
  final int M,
  final int N,
  final int P,
  final Box<int> K,
  final Box<int> L,
  final Matrix<double> A_,
  final int LDA,
  final Matrix<double> B_,
  final int LDB,
  final Array<double> ALPHA_,
  final Array<double> BETA_,
  final Matrix<double> U_,
  final int LDU,
  final Matrix<double> V_,
  final int LDV,
  final Matrix<double> Q_,
  final int LDQ,
  final Array<double> WORK_,
  final int LWORK,
  final Array<int> IWORK_,
  final Box<int> INFO,
) {
  final A = A_.having(ld: LDA);
  final B = B_.having(ld: LDB);
  final ALPHA = ALPHA_.having();
  final BETA = BETA_.having();
  final U = U_.having(ld: LDU);
  final V = V_.having(ld: LDV);
  final Q = Q_.having(ld: LDQ);
  final WORK = WORK_.having();
  final IWORK = IWORK_.having();
  bool WANTQ, WANTU, WANTV, LQUERY;
  int I, IBND, ISUB, J, LWKOPT;
  double ANORM, BNORM, SMAX, TEMP, TOLA = 0, TOLB = 0, ULP, UNFL;
  final NCYCLE = Box(0);

  // Decode and test the input parameters

  WANTU = lsame(JOBU, 'U');
  WANTV = lsame(JOBV, 'V');
  WANTQ = lsame(JOBQ, 'Q');
  LQUERY = (LWORK == -1);
  LWKOPT = 1;

  // Test the input arguments

  INFO.value = 0;
  if (!(WANTU || lsame(JOBU, 'N'))) {
    INFO.value = -1;
  } else if (!(WANTV || lsame(JOBV, 'N'))) {
    INFO.value = -2;
  } else if (!(WANTQ || lsame(JOBQ, 'N'))) {
    INFO.value = -3;
  } else if (M < 0) {
    INFO.value = -4;
  } else if (N < 0) {
    INFO.value = -5;
  } else if (P < 0) {
    INFO.value = -6;
  } else if (LDA < max(1, M)) {
    INFO.value = -10;
  } else if (LDB < max(1, P)) {
    INFO.value = -12;
  } else if (LDU < 1 || (WANTU && LDU < M)) {
    INFO.value = -16;
  } else if (LDV < 1 || (WANTV && LDV < P)) {
    INFO.value = -18;
  } else if (LDQ < 1 || (WANTQ && LDQ < N)) {
    INFO.value = -20;
  } else if (LWORK < 1 && !LQUERY) {
    INFO.value = -24;
  }

  // Compute workspace

  if (INFO.value == 0) {
    dggsvp3(JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, TOLB, K, L, U, LDU,
        V, LDV, Q, LDQ, IWORK, WORK, WORK, -1, INFO);
    LWKOPT = N + WORK[1].toInt();
    LWKOPT = max(2 * N, LWKOPT);
    LWKOPT = max(1, LWKOPT);
    WORK[1] = LWKOPT.toDouble();
  }

  if (INFO.value != 0) {
    xerbla('DGGSVD3', -INFO.value);
    return;
  }
  if (LQUERY) {
    return;
  }

  // Compute the Frobenius norm of matrices A and B

  ANORM = dlange('1', M, N, A, LDA, WORK);
  BNORM = dlange('1', P, N, B, LDB, WORK);

  // Get machine precision and set up threshold for determining
  // the effective numerical rank of the matrices A and B.

  ULP = dlamch('Precision');
  UNFL = dlamch('Safe Minimum');
  TOLA = max(M, N) * max(ANORM, UNFL) * ULP;
  TOLB = max(P, N) * max(BNORM, UNFL) * ULP;

  // Preprocessing

  dggsvp3(JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, TOLB, K, L, U, LDU,
      V, LDV, Q, LDQ, IWORK, WORK, WORK(N + 1), LWORK - N, INFO);

  // Compute the GSVD of two upper "triangular" matrices

  dtgsja(JOBU, JOBV, JOBQ, M, P, N, K.value, L.value, A, LDA, B, LDB, TOLA,
      TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, NCYCLE, INFO);

  // Sort the singular values and store the pivot indices in IWORK
  // Copy ALPHA to WORK, then sort ALPHA in WORK

  dcopy(N, ALPHA, 1, WORK, 1);
  IBND = min(L.value, M - K.value);
  for (I = 1; I <= IBND; I++) {
    // Scan for largest ALPHA[K+I]

    ISUB = I;
    SMAX = WORK[K.value + I];
    for (J = I + 1; J <= IBND; J++) {
      TEMP = WORK[K.value + J];
      if (TEMP > SMAX) {
        ISUB = J;
        SMAX = TEMP;
      }
    }
    if (ISUB != I) {
      WORK[K.value + ISUB] = WORK[K.value + I];
      WORK[K.value + I] = SMAX;
      IWORK[K.value + I] = K.value + ISUB;
    } else {
      IWORK[K.value + I] = K.value + I;
    }
  }

  WORK[1] = LWKOPT.toDouble();
}