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

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

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

static int ex1(FILE *fpin, FILE *fpout), ex2(FILE *fpin, FILE *fpout);

int main(int argc, char *argv[])
{
  FILE     *fpin, *fpout;
  Integer  exit_status_ex1 = 0;
  Integer  exit_status_ex2 = 0;
  NagError fail;

  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);

  fprintf(fpout, "nag_complex_svd (f02xec) Example Program Results\n");
  fscanf(fpin, " %*[^\n]");  /* Skip heading in data file */

  exit_status_ex1 = ex1(fpin, fpout);
  exit_status_ex2 = ex2(fpin, fpout);

  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  return (exit_status_ex1 == 0 && exit_status_ex2 == 0) ? 0 : 1;
}


#define A(I, J)  a[(I) *tda + J]
#define B(I, J)  b[(I) *tdb + J]
#define PH(I, J) ph[(I) *tdph + J]

static int ex1(FILE *fpin, FILE *fpout)
{
  Nag_Boolean wantp, wantq;
  Complex     *a = 0, *b = 0, *dummy = 0, *ph = 0;
  Integer     exit_status = 0, failinfo, i, iter, j, m, n, ncolb, tda, tdb,
              tdph;
  NagError    fail;
  double      *e = 0, *sv = 0;

  INIT_FAIL(fail);

  fprintf(fpout, "Example 1\n\n");
  fscanf(fpin, " %*[^\n]");   /* Skip heading in data file */
  if (fscanf(fpin, "%ld%ld", &m, &n) != EOF)
    {
      if (m >= 0 && n >= 0)
        {
          ncolb = 1;
          if (!(e = NAG_ALLOC(MIN(m, n)-1, double)) ||
              !(sv = NAG_ALLOC(MIN(m, n), double)) ||
              !(a = NAG_ALLOC(m*n, Complex)) ||
              !(b = NAG_ALLOC(m*ncolb, Complex)) ||
              !(ph = NAG_ALLOC(n*n, Complex)) ||
              !(dummy = NAG_ALLOC(1, Complex)))
            {
              fprintf(fpout, "Allocation failure\n");
              exit_status = -1;
              goto END;
            }
          tda = n;
          tdb = ncolb;
          tdph = n;
        }
      else
        {
          fprintf(fpout, "Invalid m or n.\n");
          exit_status = 1;
          return exit_status;
        }
      for (i = 0; i < m; ++i)
        for (j = 0; j < n; ++j)
          fscanf(fpin, "%lf%lf", COMPLEX(&A(i, j)));
      for (i = 0; i < m; ++i)
        for (j = 0; j < ncolb; ++j)
          fscanf(fpin, "%lf%lf", COMPLEX(&B(i, j)));

      /* Find the SVD of A. */

      wantq = Nag_TRUE;
      wantp = Nag_TRUE;
      /* nag_complex_svd (f02xec).
       * SVD of complex matrix
       */
      nag_complex_svd(m, n, a, tda, ncolb, b, tdb, wantq,
                      dummy, (Integer) 1, sv, wantp, ph, tdph, &iter,
                      e, &failinfo, &fail);
      if (fail.code != NE_NOERROR)
        {
          fprintf(fpout, "Error from nag_complex_svd (f02xec).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

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

  return exit_status;
}

#define A(I, J) a[(I) *tda + J]
#define Q(I, J) q[(I) *tdq + J]

static int ex2(FILE *fpin, FILE *fpout)
{
  Nag_Boolean wantp, wantq;
  Complex     *a = 0, *dummy = 0, *q = 0;
  Integer     exit_status = 0, failinfo, i, iter, j, m, n, ncolb, tda, tdq;
  NagError    fail;
  double      *e = 0, *sv = 0;

  INIT_FAIL(fail);

  fprintf(fpout, "\nExample 2\n\n");
  fscanf(fpin, " %*[^\n]");   /* Skip heading in data file */
  if (fscanf(fpin, "%ld%ld", &m, &n) != EOF)
    {
      if (m >= 0 && n >= 0)
        {
          if (!(e = NAG_ALLOC(MIN(m, n)-1, double)) ||
              !(sv = NAG_ALLOC(MIN(m, n), double)) ||
              !(a = NAG_ALLOC(m*n, Complex)) ||
              !(q = NAG_ALLOC(m*m, Complex)) ||
              !(dummy = NAG_ALLOC(1, Complex)))
            {
              fprintf(fpout, "Allocation failure\n");
              exit_status = -1;
              goto END;
            }
          tda = n;
          tdq = m;
        }
      else
        {
          fprintf(fpout, "Invalid m or n.\n");
          exit_status = 1;
          return exit_status;
        }
      for (i = 0; i < m; ++i)
        for (j = 0; j < n; ++j)
          if (fscanf(fpin, "%lf%lf", COMPLEX(&A(i, j))) != 2)
            {
              fprintf(fpout, "Data input error: program terminated.\n");
              exit_status = 1;
              goto END;
            }

      /* Find the SVD of A. */

      wantq = Nag_TRUE;
      wantp = Nag_TRUE;
      ncolb = 0;

      /* nag_complex_svd (f02xec), see above. */
      nag_complex_svd(m, n, a, tda, ncolb, dummy, (Integer) 1, wantq,
                      q, tdq, sv, wantp, dummy, (Integer) 1, &iter,
                      e, &failinfo, &fail);
      if (fail.code != NE_NOERROR)
        {
          fprintf(fpout, "Error from nag_complex_svd (f02xec).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      fprintf(fpout, "Singular value decomposition of A\n\nSingular values\n");
      for (i = 0; i < m; ++i)
        fprintf(fpout, "%9.4f%s", sv[i], (i%5 == 4 || i == m-1)?"\n":" ");
      fprintf(fpout, "\nLeft-hand singular vectors, by column\n");
      for (i = 0; i < m; ++i)
        for (j = 0; j < m; ++j)
          fprintf(fpout, "%7.4f %7.4f%s", COMPLEX(Q(i, j)),
                  (j%3 == 2 || j == n-1)?"\n":"   ");
      fprintf(fpout, "\nRight-hand singular vectors, by column\n");
      for (i = 0; i < n; ++i)
        for (j = 0; j < m; ++j)
          fprintf(fpout, "%7.4f %7.4f%s", COMPLEX_CONJ(A(j, i)),
                  (j%3 == 2 || j == n-1)?"\n":"   ");
    }
 END:
  if (e) NAG_FREE(e);
  if (sv) NAG_FREE(sv);
  if (a) NAG_FREE(a);
  if (q) NAG_FREE(q);
  if (dummy) NAG_FREE(dummy);

  return exit_status;
}