/* nag_complex_svd(f02xec) Example Program
 *
 * Copyright 1990 Numerical Algorithms Group.
 *
 * Mark 1, 1990.
 */

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

#define COMPLEX(A) A.re, A.im
#define COMPLEX_CONJ(A) A.re, -A.im

#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("f02xec 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 tdph = EX1_NMAX;

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

  Vprintf("Example 1\n\n");
  Vscanf(" %*[^\n]");   /* Skip heading in data file */
  if (scanf("%ld%ld", &m, &n) != EOF)
    {
      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
        {
          for (i=0; i<m; ++i)
            for (j=0; j<n; ++j)
              Vscanf("%lf%lf", COMPLEX(&a[i][j]));
          for (i = 0; i < m; ++i)
            Vscanf("%lf%lf", COMPLEX(&b[i]));

          /* Find the SVD of A. */

          wantq = TRUE;
          wantp = TRUE;
          f02xec(m, n, &a[0][0], tda, (Integer)1, b, (Integer)1, wantq,
                 dummy, (Integer)1, sv, wantp, &ph[0][0], tdph, &iter,
                 e, &failinfo, &fail);
          if (fail.code != NE_NOERROR) exit(EXIT_FAILURE);

          Vprintf("Singular value decomposition of A\n\nSingular values\n");
          for (i=0; i<n; ++i)
            Vprintf("%9.4f%s", sv[i], (i%5==4 || i==n-1) ? "\n": " ");
          Vprintf("\nLeft-hand singular vectors, by column\n");
          for (i=0; i<m; ++i)
            for (j=0; j<n; ++j)
              Vprintf("%7.4f %7.4f%s", COMPLEX(a[i][j]),
                      (j%3==2 || j==n-1) ? "\n": "   ");
          Vprintf("\nRight-hand singular vectors, by column\n");
          for (i=0; i<n; ++i)
            for (j=0; j<n; ++j)
              Vprintf("%7.4f %7.4f%s", COMPLEX_CONJ(ph[j][i]),
                      (j%3==2 || j==n-1) ? "\n": "   ");
          Vprintf("\nVector conjg(Q')*B\n");
          for (i=0; i<m; ++i)
            Vprintf("%7.4f %7.4f%s", COMPLEX(b[i]),
                    (i%3==2 || i==m-1) ? "\n":"   ");
        }
    }
}

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

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

  Vprintf("\nExample 2\n\n");
  Vscanf(" %*[^\n]");   /* Skip heading in data file */
  if (scanf("%ld%ld", &m, &n) != EOF)
    {
      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
        {
          for (i=0; i<m; ++i)
            for (j=0; j<n; ++j)
              if (scanf("%lf%lf", COMPLEX(&a[i][j]))!= 2)
                {
                  Vfprintf(stderr,"Data input error: program terminated.\n");
                  exit(EXIT_FAILURE);
                }

          /* Find the SVD of A. */

          wantq = TRUE;
          wantp = TRUE;
          ncolb = 0;

          f02xec(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\nSingular values\n");
          for (i=0; i<m; ++i)
            Vprintf("%9.4f%s", sv[i], (i%5==4 || i==m-1) ? "\n": " ");
          Vprintf("\nLeft-hand singular vectors, by column\n");
          for (i=0; i<m; ++i)
            for (j=0; j<m; ++j)
              Vprintf("%7.4f %7.4f%s", COMPLEX(q[i][j]),
                      (j%3==2 || j==n-1) ? "\n": "   ");
          Vprintf("\nRight-hand singular vectors, by column\n");
          for (i=0; i<n; ++i)
            for (j=0; j<m; ++j)
              Vprintf("%7.4f %7.4f%s", COMPLEX_CONJ(a[j][i]),
                      (j%3==2 || j==n-1) ? "\n": "   ");
        }
    }
}