/* nag_dgejsv (f08khc) Example Program.
 *
 * Copyright 2011 Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */

#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagf08.h>
#include <nagx02.h>
#include <nagx04.h>

int main(void)
{
  /* Scalars */
  double            eps, serrbd;
  Integer           exit_status = 0;
  Integer           pda, pdu, pdv;
  Integer           i, j, m, n, n_uvecs, n_vvecs;
  /* Arrays */
  double            *a = 0, *rcondu = 0, *rcondv = 0, *s = 0, *u = 0, *v = 0;
  double            work[7];
  Integer           iwork[3];
  char              nag_enum_arg[40];
  
  /* Nag Types */
  Nag_OrderType     order;
  Nag_Preprocess    joba;
  Nag_LeftVecsType  jobu;
  Nag_RightVecsType jobv;
  Nag_ZeroCols      jobr;
  Nag_TransType     jobt;
  Nag_Perturb       jobp;
  NagError          fail;

#ifdef NAG_COLUMN_MAJOR
#define A(I, J)   a[(J-1)*pda + I-1]
  order = Nag_ColMajor;
#else
#define A(I, J)   a[(I-1)*pda + J-1]
  order = Nag_RowMajor;
#endif

  INIT_FAIL(fail);

  printf("nag_dgejsv (f08khc) Example Program Results\n\n");

  jobu = Nag_LeftSpan;
  jobv = Nag_RightVecs;
  jobr = Nag_ZeroColsRestrict;
  jobt = Nag_NoTrans;
  jobp = Nag_PerturbOff;

  /* Skip heading in data file*/
  scanf("%*[^\n]");
  scanf("%ld%ld%*[^\n]", &m, &n);

  if (n < 0 || m < n)
    {
      printf("Invalid n or nrhs\n");
      exit_status = 1;
      goto END;;
    }

  /* Read Nag type arguments by name and convert to value */
  scanf(" %s%*[^\n]", nag_enum_arg);
  /* nag_enum_name_to_value(x04nac).
   * Converts NAG enum member name to value
   */
  joba = (Nag_Preprocess) nag_enum_name_to_value(nag_enum_arg);
  scanf(" %s%*[^\n]", nag_enum_arg);
  jobu = (Nag_LeftVecsType) nag_enum_name_to_value(nag_enum_arg);
  scanf(" %s%*[^\n]", nag_enum_arg);
  jobv = (Nag_RightVecsType) nag_enum_name_to_value(nag_enum_arg);
  scanf(" %s%*[^\n]", nag_enum_arg);
  jobr = (Nag_ZeroCols) nag_enum_name_to_value(nag_enum_arg);
  scanf(" %s%*[^\n]", nag_enum_arg);
  jobt = (Nag_TransType) nag_enum_name_to_value(nag_enum_arg);
  scanf(" %s%*[^\n]", nag_enum_arg);
  jobp = (Nag_Perturb) nag_enum_name_to_value(nag_enum_arg);

  /* Size of u and v depends on some of the above Nag type arguments. */
  n_uvecs = 1;
  if (jobu==Nag_LeftVecs) {
    n_uvecs = m;
  } else if (jobu==Nag_LeftSpan) {
    n_uvecs = n;
  } else if (jobu==Nag_NotLeftWork && jobv==Nag_RightVecs && 
	     jobt==Nag_Trans && m==n) {
    n_uvecs = m;
  }
  if (jobv==Nag_NotRightVecs) {
    n_vvecs = 1;
  } else {
    n_vvecs = n;
  }
#ifdef NAG_COLUMN_MAJOR
  pda = m;
  pdu = m;
  pdv = n;
#else
  pda = n;
  pdu = n_uvecs;
  pdv = n_vvecs;
#endif

  if (!(a      = NAG_ALLOC(m*n, double)) ||
      !(rcondu = NAG_ALLOC(m, double)) ||
      !(rcondv = NAG_ALLOC(m, double)) ||
      !(s      = NAG_ALLOC(n, double)) ||
      !(u      = NAG_ALLOC(m*n_uvecs, double)) ||
      !(v      = NAG_ALLOC(n_vvecs*n_vvecs, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Read the m by n matrix A from data file*/
  for (i = 1; i <= m; i++)
    for (j = 1; j <= n; j++) scanf("%lf", &A(i, j));
  scanf("%*[^\n]");

  /* nag_dgejsv (f08khc)
   * Compute the singular values and left and right singular vectors
   * of A (A = U*S*V^T, m>=n).
   */
  nag_dgejsv(order, joba, jobu, jobv, jobr, jobt, jobp, m, n, a, pda, s, u, pdu,
	     v, pdv, work, iwork, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_dgejsv (f08khc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Get the machine precision, eps and compute the approximate
   * error bound for the computed singular values. Note that for
   * the 2-norm, s[0] = norm(A).
   */
  eps = nag_machine_precision;
  serrbd = eps * s[0];

  /* Print (possibly scaled) singular values. */
  if (fabs(work[0] - work[1]) < 2.0 * eps)
    {
      /* No scaling required*/
      printf("Singular values\n");
      for (j = 0; j < n; j++) printf("%8.4f", s[j]);
    }
  else
    {
      printf("Scaled singular values\n");
      for (j = 0; j < n; j++) printf("%8.4f", s[j]);
      printf("\nFor true singular values, multiply by a/b,\n");
      printf("where a = %f and b = %f", work[0], work[1]);
    }
  printf("\n\n");

  /* Print left and right (spanning) singular vectors, if requested.  using
   * nag_gen_real_mat_print (x04cac)
   * Print real general matrix (easy-to-use)
   */
  if (jobu==Nag_LeftVecs || jobu==Nag_LeftSpan) {
    fflush(stdout);
    nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, m, n, u,
			   pdu, "Left singular vectors", 0, &fail);
    if (fail.code != NE_NOERROR)
      {
	printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message);
	exit_status = 1;
	goto END;
      }
  }
  if (jobv==Nag_RightVecs || jobv==Nag_RightVecsJRots) {
    printf("\n");
    fflush(stdout);
    nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, n, v,
			   pdv, "Right singular vectors", 0, &fail);
    if (fail.code != NE_NOERROR)
      {
	printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message);
	exit_status = 1;
	goto END;
      }
  }

  /* nag_ddisna (f08flc)
   * Estimate reciprocal condition numbers for the singular vectors.
   */
  nag_ddisna(Nag_LeftSingVecs, m, n, s, rcondu, &fail);
  if (fail.code == NE_NOERROR)
    nag_ddisna(Nag_RightSingVecs, m, n, s, rcondv, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_ddisna (f08flc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  if (joba==Nag_ColpivRrankCond || joba==Nag_FullpivRrankCond) {
    printf("\n\nEstimate of the condition number of column equilibrated A\n");
    printf("%11.1e", work[2]);
  }

  /* Print the approximate error bounds for the singular values and vectors. */
  printf("\n\nError estimate for the singular values\n%11.1e", serrbd);

  printf("\n\nError estimates for left singular vectors\n");
  for (i = 0; i < n; i++) printf("%11.1e", serrbd/rcondu[i]);

  printf("\n\nError estimates for right singular vectors\n");
  for (i = 0; i < n; i++) printf("%11.1e", serrbd/rcondv[i]);
  printf("\n");

 END:
  if (a)      NAG_FREE(a);
  if (rcondu) NAG_FREE(rcondu);
  if (rcondv) NAG_FREE(rcondv);
  if (s)      NAG_FREE(s);
  if (u)      NAG_FREE(u);
  if (v)      NAG_FREE(v);

  return exit_status;
}