/* nag_real_svd(f02wec) Example Program
 *
 * Copyright 1990 Numerical Algorithms Group.
 *
 * Mark 1, 1990.
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nagf02.h>

#define EX1_MMAX 20
#define EX1_NMAX 10

#define EX2_MMAX 10
#define EX2_NMAX 20

static void ex1(void), ex2(void);

int main(void)
{

  Vprintf("f02wec Example Program Results\n");
  Vscanf(" %*[^\n]");  /* Skip heading in data file */
  ex1();
  ex2();
  return EXIT_SUCCESS;
}

static void ex1(void)
{
  Integer tda = EX1_NMAX;
  Integer tdpt = EX1_NMAX;

  double a[EX1_MMAX][EX1_NMAX], b[EX1_MMAX], e[EX1_NMAX-1];
  double pt[EX1_NMAX][EX1_NMAX], sv[EX1_NMAX], dummy[1];
  Integer i, j, m, n, iter, failinfo;
  Boolean wantp, wantq;
  static NagError fail;

  Vprintf("Example 1\n");
  Vscanf(" %*[^\n]");  /* Skip Example 1 heading */
  Vscanf(" %*[^\n]");

  Vscanf("%ld%ld", &m, &n);
  if (m > EX1_MMAX || n > EX1_NMAX)
    {
      Vprintf("m or n is out of range.\n");
      Vprintf("m = %2ld, n = %2ld\n", m, n);
    }
  else
    {
      Vscanf(" %*[^\n]");
      for (i = 0; i < m; ++i)
        for (j = 0; j < n; ++j)
          Vscanf("%lf", &a[i][j]);
      Vscanf(" %*[^\n]");
      for (i = 0; i < m; ++i)
        Vscanf("%lf", &b[i]);

      /* Find the SVD of A. */
      wantq = TRUE;
      wantp = TRUE;
      fail.print = TRUE;
      f02wec(m, n, &a[0][0], tda, (Integer)1, b, (Integer)1, wantq,
             dummy, (Integer)1, sv, wantp, &pt[0][0], tdpt, &iter,
             e, &failinfo, &fail);
      if (fail.code != NE_NOERROR) exit(EXIT_FAILURE);

      Vprintf("Singular value decomposition of A\n\n");
      Vprintf("Singular values\n");
      for (i = 0; i < n; ++i)
        Vprintf(" %8.4f", sv[i]);
      Vprintf("\n\n");
      Vprintf("Left-hand singular vectors, by column\n");
      for (i = 0; i < m; ++i)
        {
          for (j = 0; j < n; ++j)
            Vprintf(" %8.4f", a[i][j]);
          Vprintf("\n");
        }
      Vprintf("\n");
      Vprintf("Right-hand singular vectors, by column\n");
      for (i = 0; i < n; ++i)
        {
          for (j = 0; j < n; ++j)
            Vprintf(" %8.4f", pt[j][i]);
          Vprintf("\n");
        }
      Vprintf("\n");
      Vprintf("Vector Q'*B\n");
      for (i = 0; i < m; ++i)
        Vprintf(" %8.4f", b[i]);
      Vprintf("\n\n");
    }
}


static void ex2(void)
{
  Integer tda = EX2_NMAX;
  Integer tdq = EX2_MMAX;

  double a[EX2_MMAX][EX2_NMAX], e[EX2_NMAX-1];
  double q[EX2_MMAX][EX2_MMAX], sv[EX2_MMAX], dummy[1];
  Integer i, j, m, n, iter, ncolb, failinfo;
  Boolean wantp, wantq;
  static NagError fail;

  Vprintf("\nExample 2\n");
  Vscanf(" %*[^\n]");  /* Skip Example 2 heading */
  Vscanf(" %*[^\n]");

  Vscanf("%ld%ld", &m, &n);
  if (m > EX2_MMAX || n > EX2_NMAX)
    {
      Vprintf("m or n is out of range.\n");
      Vprintf("m = %2ld, n = %2ld\n", m, n);
    }
  else
    {
      Vscanf(" %*[^\n]");
      for (i = 0; i < m; ++i)
        for (j = 0; j < n; ++j)
          Vscanf("%lf", &a[i][j]);

      /* Find the SVD of A. */
      wantq = TRUE;
      wantp = TRUE;
      ncolb = 0;
      fail.print = TRUE;
      f02wec(m, n, &a[0][0], tda, ncolb, dummy, (Integer)1, wantq,
             &q[0][0], tdq, sv, wantp, dummy, (Integer)1, &iter,
             e, &failinfo, &fail);
      if (fail.code != NE_NOERROR) exit(EXIT_FAILURE);

      Vprintf("Singular value decomposition of A\n\n\n");
      Vprintf("Singular values\n\n");
      for (i = 0; i < m; ++i)
        Vprintf(" %8.4f", sv[i]);
      Vprintf("\n\n");
      Vprintf("Left-hand singular vectors, by column\n\n");
      for (i = 0; i < m; ++i)
        {
          for (j = 0; j < m; ++j)
            Vprintf(" %8.4f", q[i][j]);
          Vprintf("\n");
        }
      Vprintf("Right-hand singular vectors, by column\n\n");
      for (i = 0; i < n; ++i)
        {
          for (j = 0; j < m; ++j)
            Vprintf(" %8.4f", a[j][i]);
          Vprintf("\n");
        }
    }
}