/* 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>

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


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

int main(void)
{

  Integer exit_status_ex1=0;
  Integer exit_status_ex2=0;

  Vprintf("nag_complex_svd (f02xec) Example Program Results\n");
  Vscanf(" %*[^\n]");  /* Skip heading in data file */

  exit_status_ex1 = ex1();
  exit_status_ex2 = ex2();

  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(void)
{
  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);


  Vprintf("Example 1\n\n");
  Vscanf(" %*[^\n]");   /* Skip heading in data file */
  if (scanf("%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)) )
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }
          tda = n;
          tdb = ncolb;
          tdph = n;
        }
      else
        {
          Vprintf("Invalid m or n.\n");
          exit_status = 1;
          return exit_status;
        }
      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)
        for (j=0; j<ncolb; ++j)
          Vscanf("%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)
        {
          Vprintf("Error from nag_complex_svd (f02xec).\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }

      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)
        for (j=0; j<ncolb; ++j)
          Vprintf("%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(void)
{
  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);


  Vprintf("\nExample 2\n\n");
  Vscanf(" %*[^\n]");   /* Skip heading in data file */
  if (scanf("%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)) )
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }
          tda = n;
          tdq = m;
        }
      else
        {
          Vprintf("Invalid m or n.\n");
          exit_status = 1;
          return exit_status;
        }
      for (i=0; i<m; ++i)
        for (j=0; j<n; ++j)
          if (scanf("%lf%lf", COMPLEX(&A(i,j)))!= 2)
            {
              Vprintf("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)
        {
          Vprintf("Error from nag_complex_svd (f02xec).\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }

      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": "   ");

    }
 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;
}