/* nag_lapackeig_dbdsdc (f08mdc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>

int main(void) {
  /* Scalars */
  double alpha, beta, eps, norm;
  Integer abi, i, j, k1, k2, leniq, lenq, mlvl, n, pdab, pdb, pdu, pdvt;
  Integer exit_status = 0, smlsiz = 25;

  /* Arrays */
  double *ab = 0, *b = 0, *d = 0, *e = 0, *q = 0, *u = 0, *vt = 0;
  Integer *iq = 0;
  char nag_enum_arg[40];

  /* Nag Types */
  NagError fail;
  Nag_OrderType order;
  Nag_UploType uplo;
  Nag_ComputeSingularVecsType compq;

#ifdef NAG_COLUMN_MAJOR
#define B(I, J) b[(J - 1) * pdb + I - 1]
#define U(I, J) u[(J - 1) * pdu + I - 1]
  order = Nag_ColMajor;
#else
#define B(I, J) b[(I - 1) * pdb + J - 1]
#define U(I, J) u[(I - 1) * pdu + J - 1]
  order = Nag_RowMajor;
#endif

  INIT_FAIL(fail);

  printf("nag_lapackeig_dbdsdc (f08mdc) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n]");
  scanf("%" NAG_IFMT "%*[^\n]", &n);
  if (n < 0) {
    printf("Invalid n\n");
    exit_status = 1;
    goto END;
    ;
  }
  scanf(" %39s%*[^\n]", nag_enum_arg);
  /* nag_enum_name_to_value (x04nac).
   * Converts NAG enum member name to value
   */
  uplo = (Nag_UploType)nag_enum_name_to_value(nag_enum_arg);
  scanf(" %39s%*[^\n]", nag_enum_arg);
  /* Starting index for main diagonal in banded storage format = abi. */
  if ((order == Nag_ColMajor && uplo == Nag_Lower) ||
      (order == Nag_RowMajor && uplo == Nag_Upper)) {
    abi = 0;
  } else {
    abi = 1;
  }

  compq = (Nag_ComputeSingularVecsType)nag_enum_name_to_value(nag_enum_arg);
  /* size of u, vt, q and iq depends on value of compq input */
  if (compq == Nag_SingularVecs) {
    pdu = n;
    pdvt = n;
  } else {
    pdu = 1;
    pdvt = 1;
  }
  if (compq == Nag_PackedSingularVecs) {
    mlvl = (Integer)(log(n / (smlsiz + 1.0)) / log(2.0)) + 1;
    if (mlvl < 1)
      mlvl = 1;
    lenq = MAX(n * n + 5 * n, n * (3 + 2 * smlsiz + 8 * mlvl));
    leniq = n * 3 * mlvl;
  } else {
    lenq = 1;
    leniq = 1;
  }

  pdb = n;
  pdab = 2;
  /* Allocate memory */
  if (!(b = NAG_ALLOC(n * n, double)) || !(ab = NAG_ALLOC(pdab * n, double)) ||
      !(d = NAG_ALLOC(n, double)) || !(e = NAG_ALLOC(n - 1, double)) ||
      !(q = NAG_ALLOC(lenq, double)) || !(u = NAG_ALLOC(pdu * pdu, double)) ||
      !(vt = NAG_ALLOC(pdvt * pdvt, double)) ||
      !(iq = NAG_ALLOC(leniq, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read the bidiagonal matrix B from data file,
   * first the diagonal elements, and then the off-diagonal elements.
   */
  for (i = 0; i < n; ++i)
    scanf("%lf", &d[i]);
  scanf("%*[^\n]");
  for (i = 0; i < n - 1; ++i)
    scanf("%lf", &e[i]);
  scanf("%*[^\n]");

  /* Store diagonal arrays in banded format in ab for printing */
  for (i = 0; i < n; i++)
    ab[2 * i + abi] = d[i];
  for (i = 0; i < n - 1; i++)
    ab[2 * i + abi + 1] = e[i];

  /* k1 = lower bandwidth, k2 = upper bandwidth */
  k1 = (uplo == Nag_Upper ? 0 : 1);
  k2 = 1 - k1;

  /* Print Bidiagonal Matrix B stored in ab.
   * nag_file_print_matrix_real_band (x04cec).
   * Print real packed banded matrix (easy-to-use)
   */
  fflush(stdout);
  nag_file_print_matrix_real_band(order, n, n, k1, k2, ab, 2, "Matrix B", 0,
                                  &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_band (x04cec).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Calculate the singular values and left and right singular vectors of B
   * using nag_lapackeig_dbdsdc (f08mdc).
   */
  nag_lapackeig_dbdsdc(order, uplo, compq, n, d, e, u, pdu, vt, pdvt, q, iq,
                       &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_lapackeig_dbdsdc (f08mdc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\nSingular values\n");
  for (i = 0; i < n; ++i)
    printf(" %8.4f%s", d[i], i % 8 == 7 ? "\n" : "");
  printf("\n\n");

  if (compq == Nag_SingularVecs) {
    /* Reconstruct bidiagonal matrix from decomposition:
     * first, U <- U*S, then Compute B = U*S*V^T.
     * nag_blast_dgemm (f16yac).
     */
    for (i = 1; i <= n; i++)
      for (j = 1; j <= n; j++)
        U(i, j) = U(i, j) * d[j - 1];
    alpha = 1.0;
    beta = 0.0;
    nag_blast_dgemm(order, Nag_NoTrans, Nag_NoTrans, n, n, n, alpha, u, pdu, vt,
                    pdvt, beta, b, pdb, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_blast_dgemm (f16yac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    /* Subtract original bidiagonal matrix:
     * this should give a matrix close to zero.
     */
    for (i = 1; i <= n; i++)
      B(i, i) -= ab[2 * i - 2 + abi];
    if (uplo == Nag_Upper)
      for (i = 1; i <= n - 1; i++)
        B(i, i + 1) -= ab[2 * i - 1 + abi];
    else
      for (i = 1; i <= n - 1; i++)
        B(i + 1, i) -= ab[2 * i - 1 + abi];

    /* nag_blast_dge_norm (f16rac): Find norm of matrix B and print warning if
     * it is too large.
     */
    nag_blast_dge_norm(order, Nag_OneNorm, n, n, b, pdb, &norm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_blast_dge_norm (f16rac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    /* Get the machine precision, using nag_machine_precision (x02ajc) */
    eps = nag_machine_precision;
    if (norm > pow(eps, 0.8)) {
      printf("Norm of B-(U*S*V^T) is much greater than 0.\nSchur "
             "factorization has failed.\n norm = %13.4e\n",
             norm);
    }
  }
END:
  NAG_FREE(ab);
  NAG_FREE(b);
  NAG_FREE(d);
  NAG_FREE(e);
  NAG_FREE(q);
  NAG_FREE(u);
  NAG_FREE(vt);
  NAG_FREE(iq);

  return exit_status;
}