/* nag_real_form_q(f01qec) Example Program
 *
 * Copyright 1990 Numerical Algorithms Group.
 *
 * Mark 1, 1990.
 */

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

#define MMAX 20
#define NMAX 10

int main(void)
{
  Integer tda = NMAX;
  Integer tdq = MMAX;

  double zeta[NMAX], a[MMAX][NMAX], q[MMAX][MMAX];
  Integer i, j, m, n, ncolq;

  Vprintf("f01qec Example Program Results\n");
  Vscanf(" %*[^\n]"); /* skip headings in data file */
  Vscanf(" %*[^\n]");

  Vscanf("%ld%ld", &m, &n);
  if (m > MMAX || n > 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)   /* Read matrix data */
        for (j = 0; j < n; ++j)
          Vscanf("%lf", &a[i][j]);

      /* Find the QR factorization of A */
      f01qcc(m, n, &a[0][0], tda, zeta, NAGERR_DEFAULT);

      /* Copy the array a into q and form the m by m matrix Q */
      for (j = 0; j < n; ++j)
        for (i = 0; i < m; ++i)
          q[i][j] = a[i][j];
      ncolq = m;
      f01qec(Nag_ElementsSeparate, m, n, ncolq, &q[0][0], tdq,
             zeta, NAGERR_DEFAULT);

      Vprintf("Matrix  Q\n");
      for (i = 0; i < m; ++i)
        {
          for (j = 0; j < ncolq; ++j)
            Vprintf(" %8.4f", q[i][j]);
          Vprintf("\n");
        }
    }
  return EXIT_SUCCESS;
}