ztprfb function

void ztprfb(
  1. String SIDE,
  2. String TRANS,
  3. String DIRECT,
  4. String STOREV,
  5. int M,
  6. int N,
  7. int K,
  8. int L,
  9. Matrix<Complex> V_,
  10. int LDV,
  11. Matrix<Complex> T_,
  12. int LDT,
  13. Matrix<Complex> A_,
  14. int LDA,
  15. Matrix<Complex> B_,
  16. int LDB,
  17. Matrix<Complex> WORK_,
  18. int LDWORK,
)

Implementation

void ztprfb(
  final String SIDE,
  final String TRANS,
  final String DIRECT,
  final String STOREV,
  final int M,
  final int N,
  final int K,
  final int L,
  final Matrix<Complex> V_,
  final int LDV,
  final Matrix<Complex> T_,
  final int LDT,
  final Matrix<Complex> A_,
  final int LDA,
  final Matrix<Complex> B_,
  final int LDB,
  final Matrix<Complex> WORK_,
  final int LDWORK,
) {
  final V = V_.having(ld: LDV);
  final T = T_.having(ld: LDT);
  final A = A_.having(ld: LDA);
  final B = B_.having(ld: LDB);
  final WORK = WORK_.having(ld: LDWORK);
  int I, J, MP, NP, KP;
  bool LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW;

  // Quick return if possible
  if (M <= 0 || N <= 0 || K <= 0 || L < 0) return;

  if (lsame(STOREV, 'C')) {
    COLUMN = true;
    ROW = false;
  } else if (lsame(STOREV, 'R')) {
    COLUMN = false;
    ROW = true;
  } else {
    COLUMN = false;
    ROW = false;
  }

  if (lsame(SIDE, 'L')) {
    LEFT = true;
    RIGHT = false;
  } else if (lsame(SIDE, 'R')) {
    LEFT = false;
    RIGHT = true;
  } else {
    LEFT = false;
    RIGHT = false;
  }

  if (lsame(DIRECT, 'F')) {
    FORWARD = true;
    BACKWARD = false;
  } else if (lsame(DIRECT, 'B')) {
    FORWARD = false;
    BACKWARD = true;
  } else {
    FORWARD = false;
    BACKWARD = false;
  }

  if (COLUMN && FORWARD && LEFT) {
    // Let  W =  [ I ]    (K-by-K)
    //           [ V ]    (M-by-K)
    //
    // Form  H C  or  H**H C  where  C = [ A ]  (K-by-N)
    //                                   [ B ]  (M-by-N)
    //
    // H = I - W T W**H          or  H**H = I - W T**H W**H
    //
    // A -=   T (A + V**H B)  or  A -=   T**H (A + V**H B)
    // B -= V T (A + V**H B)  or  B -= V T**H (A + V**H B)

    MP = min(M - L + 1, M);
    KP = min(L + 1, K);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= L; I++) {
        WORK[I][J] = B[M - L + I][J];
      }
    }
    ztrmm('L', 'U', 'C', 'N', L, N, Complex.one, V(MP, 1), LDV, WORK, LDWORK);
    zgemm('C', 'N', L, N, M - L, Complex.one, V, LDV, B, LDB, Complex.one, WORK,
        LDWORK);
    zgemm('C', 'N', K - L, N, M, Complex.one, V(1, KP), LDV, B, LDB,
        Complex.zero, WORK(KP, 1), LDWORK);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= K; I++) {
        WORK[I][J] += A[I][J];
      }
    }

    ztrmm('L', 'U', TRANS, 'N', K, N, Complex.one, T, LDT, WORK, LDWORK);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= K; I++) {
        A[I][J] -= WORK[I][J];
      }
    }

    zgemm('N', 'N', M - L, N, K, -Complex.one, V, LDV, WORK, LDWORK,
        Complex.one, B, LDB);
    zgemm('N', 'N', L, N, K - L, -Complex.one, V(MP, KP), LDV, WORK(KP, 1),
        LDWORK, Complex.one, B(MP, 1), LDB);
    ztrmm('L', 'U', 'N', 'N', L, N, Complex.one, V(MP, 1), LDV, WORK, LDWORK);
    for (J = 1; J <= N; J++) {
      for (I = 1; I <= L; I++) {
        B[M - L + I][J] -= WORK[I][J];
      }
    }
  } else if (COLUMN && FORWARD && RIGHT) {
    // Let  W =  [ I ]    (K-by-K)
    //           [ V ]    (N-by-K)
    //
    // Form  C H or  C H**H  where  C = [ A B ] (A is M-by-K, B is M-by-N)
    //
    // H = I - W T W**H          or  H**H = I - W T**H W**H
    //
    // A -= (A + B V) T      or  A -= (A + B V) T**H
    // B -= (A + B V) T V**H  or  B -= (A + B V) T**H V**H

    NP = min(N - L + 1, N);
    KP = min(L + 1, K);

    for (J = 1; J <= L; J++) {
      for (I = 1; I <= M; I++) {
        WORK[I][J] = B[I][N - L + J];
      }
    }
    ztrmm('R', 'U', 'N', 'N', M, L, Complex.one, V(NP, 1), LDV, WORK, LDWORK);
    zgemm('N', 'N', M, L, N - L, Complex.one, B, LDB, V, LDV, Complex.one, WORK,
        LDWORK);
    zgemm('N', 'N', M, K - L, N, Complex.one, B, LDB, V(1, KP), LDV,
        Complex.zero, WORK(1, KP), LDWORK);

    for (J = 1; J <= K; J++) {
      for (I = 1; I <= M; I++) {
        WORK[I][J] += A[I][J];
      }
    }

    ztrmm('R', 'U', TRANS, 'N', M, K, Complex.one, T, LDT, WORK, LDWORK);

    for (J = 1; J <= K; J++) {
      for (I = 1; I <= M; I++) {
        A[I][J] -= WORK[I][J];
      }
    }

    zgemm('N', 'C', M, N - L, K, -Complex.one, WORK, LDWORK, V, LDV,
        Complex.one, B, LDB);
    zgemm('N', 'C', M, L, K - L, -Complex.one, WORK(1, KP), LDWORK, V(NP, KP),
        LDV, Complex.one, B(1, NP), LDB);
    ztrmm('R', 'U', 'C', 'N', M, L, Complex.one, V(NP, 1), LDV, WORK, LDWORK);
    for (J = 1; J <= L; J++) {
      for (I = 1; I <= M; I++) {
        B[I][N - L + J] -= WORK[I][J];
      }
    }
  } else if (COLUMN && BACKWARD && LEFT) {
    // Let  W =  [ V ]    (M-by-K)
    //           [ I ]    (K-by-K)
    //
    // Form  H C  or  H**H C  where  C = [ B ]  (M-by-N)
    //                                   [ A ]  (K-by-N)
    //
    // H = I - W T W**H          or  H**H = I - W T**H W**H
    //
    // A -=   T (A + V**H B)  or  A -=   T**H (A + V**H B)
    // B -= V T (A + V**H B)  or  B -= V T**H (A + V**H B)

    MP = min(L + 1, M);
    KP = min(K - L + 1, K);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= L; I++) {
        WORK[K - L + I][J] = B[I][J];
      }
    }

    ztrmm('L', 'L', 'C', 'N', L, N, Complex.one, V(1, KP), LDV, WORK(KP, 1),
        LDWORK);
    zgemm('C', 'N', L, N, M - L, Complex.one, V(MP, KP), LDV, B(MP, 1), LDB,
        Complex.one, WORK(KP, 1), LDWORK);
    zgemm('C', 'N', K - L, N, M, Complex.one, V, LDV, B, LDB, Complex.zero,
        WORK, LDWORK);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= K; I++) {
        WORK[I][J] += A[I][J];
      }
    }

    ztrmm('L', 'L', TRANS, 'N', K, N, Complex.one, T, LDT, WORK, LDWORK);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= K; I++) {
        A[I][J] -= WORK[I][J];
      }
    }

    zgemm('N', 'N', M - L, N, K, -Complex.one, V(MP, 1), LDV, WORK, LDWORK,
        Complex.one, B(MP, 1), LDB);
    zgemm('N', 'N', L, N, K - L, -Complex.one, V, LDV, WORK, LDWORK,
        Complex.one, B, LDB);
    ztrmm('L', 'L', 'N', 'N', L, N, Complex.one, V(1, KP), LDV, WORK(KP, 1),
        LDWORK);
    for (J = 1; J <= N; J++) {
      for (I = 1; I <= L; I++) {
        B[I][J] -= WORK[K - L + I][J];
      }
    }
  } else if (COLUMN && BACKWARD && RIGHT) {
    // Let  W =  [ V ]    (N-by-K)
    //           [ I ]    (K-by-K)
    //
    // Form  C H  or  C H**H  where  C = [ B A ] (B is M-by-N, A is M-by-K)
    //
    // H = I - W T W**H          or  H**H = I - W T**H W**H
    //
    // A -= (A + B V) T      or  A -= (A + B V) T**H
    // B -= (A + B V) T V**H  or  B -= (A + B V) T**H V**H

    NP = min(L + 1, N);
    KP = min(K - L + 1, K);

    for (J = 1; J <= L; J++) {
      for (I = 1; I <= M; I++) {
        WORK[I][K - L + J] = B[I][J];
      }
    }
    ztrmm('R', 'L', 'N', 'N', M, L, Complex.one, V(1, KP), LDV, WORK(1, KP),
        LDWORK);
    zgemm('N', 'N', M, L, N - L, Complex.one, B(1, NP), LDB, V(NP, KP), LDV,
        Complex.one, WORK(1, KP), LDWORK);
    zgemm('N', 'N', M, K - L, N, Complex.one, B, LDB, V, LDV, Complex.zero,
        WORK, LDWORK);

    for (J = 1; J <= K; J++) {
      for (I = 1; I <= M; I++) {
        WORK[I][J] += A[I][J];
      }
    }

    ztrmm('R', 'L', TRANS, 'N', M, K, Complex.one, T, LDT, WORK, LDWORK);

    for (J = 1; J <= K; J++) {
      for (I = 1; I <= M; I++) {
        A[I][J] -= WORK[I][J];
      }
    }

    zgemm('N', 'C', M, N - L, K, -Complex.one, WORK, LDWORK, V(NP, 1), LDV,
        Complex.one, B(1, NP), LDB);
    zgemm('N', 'C', M, L, K - L, -Complex.one, WORK, LDWORK, V, LDV,
        Complex.one, B, LDB);
    ztrmm('R', 'L', 'C', 'N', M, L, Complex.one, V(1, KP), LDV, WORK(1, KP),
        LDWORK);
    for (J = 1; J <= L; J++) {
      for (I = 1; I <= M; I++) {
        B[I][J] -= WORK[I][K - L + J];
      }
    }
  } else if (ROW && FORWARD && LEFT) {
    // Let  W =  [ I V ] ( I is K-by-K, V is K-by-M )
    //
    // Form  H C  or  H**H C  where  C = [ A ]  (K-by-N)
    //                                   [ B ]  (M-by-N)
    //
    // H = I - W**H T W          or  H**H = I - W**H T**H W
    //
    // A -=     T (A + V B)  or  A -=     T**H (A + V B)
    // B -= V**H T (A + V B)  or  B -= V**H T**H (A + V B)

    MP = min(M - L + 1, M);
    KP = min(L + 1, K);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= L; I++) {
        WORK[I][J] = B[M - L + I][J];
      }
    }
    ztrmm('L', 'L', 'N', 'N', L, N, Complex.one, V(1, MP), LDV, WORK, LDB);
    zgemm('N', 'N', L, N, M - L, Complex.one, V, LDV, B, LDB, Complex.one, WORK,
        LDWORK);
    zgemm('N', 'N', K - L, N, M, Complex.one, V(KP, 1), LDV, B, LDB,
        Complex.zero, WORK(KP, 1), LDWORK);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= K; I++) {
        WORK[I][J] += A[I][J];
      }
    }

    ztrmm('L', 'U', TRANS, 'N', K, N, Complex.one, T, LDT, WORK, LDWORK);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= K; I++) {
        A[I][J] -= WORK[I][J];
      }
    }

    zgemm('C', 'N', M - L, N, K, -Complex.one, V, LDV, WORK, LDWORK,
        Complex.one, B, LDB);
    zgemm('C', 'N', L, N, K - L, -Complex.one, V(KP, MP), LDV, WORK(KP, 1),
        LDWORK, Complex.one, B(MP, 1), LDB);
    ztrmm('L', 'L', 'C', 'N', L, N, Complex.one, V(1, MP), LDV, WORK, LDWORK);
    for (J = 1; J <= N; J++) {
      for (I = 1; I <= L; I++) {
        B[M - L + I][J] -= WORK[I][J];
      }
    }
  } else if (ROW && FORWARD && RIGHT) {
    // Let  W =  [ I V ] ( I is K-by-K, V is K-by-N )
    //
    // Form  C H  or  C H**H  where  C = [ A B ] (A is M-by-K, B is M-by-N)
    //
    // H = I - W**H T W            or  H**H = I - W**H T**H W
    //
    // A -= (A + B V**H) T      or  A -= (A + B V**H) T**H
    // B -= (A + B V**H) T V    or  B -= (A + B V**H) T**H V

    NP = min(N - L + 1, N);
    KP = min(L + 1, K);

    for (J = 1; J <= L; J++) {
      for (I = 1; I <= M; I++) {
        WORK[I][J] = B[I][N - L + J];
      }
    }
    ztrmm('R', 'L', 'C', 'N', M, L, Complex.one, V(1, NP), LDV, WORK, LDWORK);
    zgemm('N', 'C', M, L, N - L, Complex.one, B, LDB, V, LDV, Complex.one, WORK,
        LDWORK);
    zgemm('N', 'C', M, K - L, N, Complex.one, B, LDB, V(KP, 1), LDV,
        Complex.zero, WORK(1, KP), LDWORK);

    for (J = 1; J <= K; J++) {
      for (I = 1; I <= M; I++) {
        WORK[I][J] += A[I][J];
      }
    }

    ztrmm('R', 'U', TRANS, 'N', M, K, Complex.one, T, LDT, WORK, LDWORK);

    for (J = 1; J <= K; J++) {
      for (I = 1; I <= M; I++) {
        A[I][J] -= WORK[I][J];
      }
    }

    zgemm('N', 'N', M, N - L, K, -Complex.one, WORK, LDWORK, V, LDV,
        Complex.one, B, LDB);
    zgemm('N', 'N', M, L, K - L, -Complex.one, WORK(1, KP), LDWORK, V(KP, NP),
        LDV, Complex.one, B(1, NP), LDB);
    ztrmm('R', 'L', 'N', 'N', M, L, Complex.one, V(1, NP), LDV, WORK, LDWORK);
    for (J = 1; J <= L; J++) {
      for (I = 1; I <= M; I++) {
        B[I][N - L + J] -= WORK[I][J];
      }
    }
  } else if (ROW && BACKWARD && LEFT) {
    // Let  W =  [ V I ] ( I is K-by-K, V is K-by-M )
    //
    // Form  H C  or  H**H C  where  C = [ B ]  (M-by-N)
    //                                   [ A ]  (K-by-N)
    //
    // H = I - W**H T W          or  H**H = I - W**H T**H W
    //
    // A -=     T (A + V B)  or  A -=     T**H (A + V B)
    // B -= V**H T (A + V B)  or  B -= V**H T**H (A + V B)

    MP = min(L + 1, M);
    KP = min(K - L + 1, K);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= L; I++) {
        WORK[K - L + I][J] = B[I][J];
      }
    }
    ztrmm('L', 'U', 'N', 'N', L, N, Complex.one, V(KP, 1), LDV, WORK(KP, 1),
        LDWORK);
    zgemm('N', 'N', L, N, M - L, Complex.one, V(KP, MP), LDV, B(MP, 1), LDB,
        Complex.one, WORK(KP, 1), LDWORK);
    zgemm('N', 'N', K - L, N, M, Complex.one, V, LDV, B, LDB, Complex.zero,
        WORK, LDWORK);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= K; I++) {
        WORK[I][J] += A[I][J];
      }
    }

    ztrmm('L', 'L ', TRANS, 'N', K, N, Complex.one, T, LDT, WORK, LDWORK);

    for (J = 1; J <= N; J++) {
      for (I = 1; I <= K; I++) {
        A[I][J] -= WORK[I][J];
      }
    }

    zgemm('C', 'N', M - L, N, K, -Complex.one, V(1, MP), LDV, WORK, LDWORK,
        Complex.one, B(MP, 1), LDB);
    zgemm('C', 'N', L, N, K - L, -Complex.one, V, LDV, WORK, LDWORK,
        Complex.one, B, LDB);
    ztrmm('L', 'U', 'C', 'N', L, N, Complex.one, V(KP, 1), LDV, WORK(KP, 1),
        LDWORK);
    for (J = 1; J <= N; J++) {
      for (I = 1; I <= L; I++) {
        B[I][J] -= WORK[K - L + I][J];
      }
    }
  } else if (ROW && BACKWARD && RIGHT) {
    // Let  W =  [ V I ] ( I is K-by-K, V is K-by-N )
    //
    // Form  C H  or  C H**H  where  C = [ B A ] (A is M-by-K, B is M-by-N)
    //
    // H = I - W**H T W            or  H**H = I - W**H T**H W
    //
    // A -= (A + B V**H) T      or  A -= (A + B V**H) T**H
    // B -= (A + B V**H) T V    or  B -= (A + B V**H) T**H V

    NP = min(L + 1, N);
    KP = min(K - L + 1, K);

    for (J = 1; J <= L; J++) {
      for (I = 1; I <= M; I++) {
        WORK[I][K - L + J] = B[I][J];
      }
    }
    ztrmm('R', 'U', 'C', 'N', M, L, Complex.one, V(KP, 1), LDV, WORK(1, KP),
        LDWORK);
    zgemm('N', 'C', M, L, N - L, Complex.one, B(1, NP), LDB, V(KP, NP), LDV,
        Complex.one, WORK(1, KP), LDWORK);
    zgemm('N', 'C', M, K - L, N, Complex.one, B, LDB, V, LDV, Complex.zero,
        WORK, LDWORK);

    for (J = 1; J <= K; J++) {
      for (I = 1; I <= M; I++) {
        WORK[I][J] += A[I][J];
      }
    }

    ztrmm('R', 'L', TRANS, 'N', M, K, Complex.one, T, LDT, WORK, LDWORK);

    for (J = 1; J <= K; J++) {
      for (I = 1; I <= M; I++) {
        A[I][J] -= WORK[I][J];
      }
    }

    zgemm('N', 'N', M, N - L, K, -Complex.one, WORK, LDWORK, V(1, NP), LDV,
        Complex.one, B(1, NP), LDB);
    zgemm('N', 'N', M, L, K - L, -Complex.one, WORK, LDWORK, V, LDV,
        Complex.one, B, LDB);
    ztrmm('R', 'U', 'N', 'N', M, L, Complex.one, V(KP, 1), LDV, WORK(1, KP),
        LDWORK);
    for (J = 1; J <= L; J++) {
      for (I = 1; I <= M; I++) {
        B[I][J] -= WORK[I][K - L + J];
      }
    }
  }
}