/* nag_zgeqrf (f08asc) Example Program.

 *

 * Copyright 2001 Numerical Algorithms Group.

 *

 * Mark 7, 2001.

 */



#include <stdio.h>

#include <nag.h>

#include <nag_stdlib.h>

#include <nagf07.h>

#include <nagf08.h>

#include <nagx04.h>



int main(void)

{

  /* Scalars */

  Integer  i, j, m, n, nrhs, pda, pdb, tau_len;

  Integer  exit_status=0;

  NagError fail;

  Nag_OrderType order;

  /* Arrays */

  Complex *a=0, *b=0, *tau=0;



#ifdef NAG_COLUMN_MAJOR

#define A(I,J) a[(J-1)*pda + I - 1]

#define B(I,J) b[(J-1)*pdb + I - 1]

  order = Nag_ColMajor;

#else

#define A(I,J) a[(I-1)*pda + J - 1]

#define B(I,J) b[(I-1)*pdb + J - 1]

  order = Nag_RowMajor;

#endif



  INIT_FAIL(fail);

  Vprintf("nag_zgeqrf (f08asc) Example Program Results\n");

  /* Skip heading in data file */

  Vscanf("%*[^\n] ");

  Vscanf("%ld%ld%ld%*[^\n] ", &m, &n, &nrhs);

#ifdef NAG_COLUMN_MAJOR

  pda = m;

  pdb = m;

#else

  pda = n;

  pdb = nrhs;

#endif

  tau_len = MIN(m,n);



  /* Allocate memory */

  if ( !(a = NAG_ALLOC(m * n, Complex)) ||

       !(b = NAG_ALLOC(m * nrhs, Complex)) ||

       !(tau = NAG_ALLOC(tau_len, Complex)) )

    {

      Vprintf("Allocation failure\n");

      exit_status = -1;

      goto END;

    }

  /* Read A and B from data file */

  for (i = 1; i <= m; ++i)

    {

      for (j = 1; j <= n; ++j)

        Vscanf(" ( %lf , %lf )", &A(i,j).re, &A(i,j).im);

    }

  Vscanf("%*[^\n] ");

  for (i = 1; i <= m; ++i)

    {

      for (j = 1; j <= nrhs; ++j)

        Vscanf(" ( %lf , %lf )", &B(i,j).re, &B(i,j).im);

    }

  Vscanf("%*[^\n] ");



  /* Compute the QR factorization of A */

  /* nag_zgeqrf (f08asc).

   * QR factorization of complex general rectangular matrix

   */

  nag_zgeqrf(order, m, n, a, pda, tau, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_zgeqrf (f08asc).\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }



  /* Compute C = (Q**H)*B, storing the result in B */

  /* nag_zunmqr (f08auc).

   * Apply unitary transformation determined by nag_zgeqrf

   * (f08asc) or nag_zgeqpf (f08bsc)

   */

  nag_zunmqr(order, Nag_LeftSide, Nag_ConjTrans, m, nrhs, n, a, pda,

             tau, b, pdb, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_zunmqr (f08auc).\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }



  /* Compute least-squares solution by backsubstitution in R*X = C */

  /* nag_ztrtrs (f07tsc).

   * Solution of complex triangular system of linear

   * equations, multiple right-hand sides

   */

  nag_ztrtrs(order, Nag_Upper, Nag_NoTrans, Nag_NonUnitDiag, n, nrhs,

             a, pda, b, pdb, &fail);



  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_ztrtrs (f07tsc).\n%s\n", fail.message);

      exit_status = 1;

      goto END;

    }

  /* Print least-squares solution(s) */

  Vprintf("\n");

  /* nag_gen_complx_mat_print_comp (x04dbc).

   * Print complex general matrix (comprehensive)

   */

  nag_gen_complx_mat_print_comp(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n,

                                nrhs, b, pdb, Nag_BracketForm, "%7.4f",

                                "Least-squares solution(s)",

                                Nag_IntegerLabels, 0, Nag_IntegerLabels, 0, 80,

                                0, 0, &fail);

  if (fail.code != NE_NOERROR)

    {

      Vprintf("Error from nag_gen_complx_mat_print_comp (x04dbc).\n%s\n",

              fail.message);

      exit_status = 1;

      goto END;

    }

 END:

  if (a) NAG_FREE(a);

  if (b) NAG_FREE(b);

  if (tau) NAG_FREE(tau);

  return exit_status;

}