/* nag_rand_orthog_matrix (g05pxc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 9, 2009.
 */
/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg05.h>

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

int main(void)
{
  /* Integer scalar and array declarations */
  Integer         exit_status = 0;
  Integer         i, j, lstate, a_size;
  Integer         *state = 0;
  Integer         pda;

  /* NAG structures */
  NagError        fail;

  /* Double scalar and array declarations */
  double          *a = 0;

  /* Set the size of the matrix to be generated */
  Integer         n = 4;
  Integer         m = 4;

  /* Set A to the identity matrix on input */
  Nag_InitializeA init = Nag_InitializeI;

  /* Multiple on the right hand side */
  Nag_SideType    side = Nag_RightSide;

  /* Choose the base generator */
  Nag_BaseRNG     genid = Nag_Basic;
  Integer         subid = 0;

  /* Set the seed */
  Integer         seed[] = { 1762543 };
  Integer         lseed = 1;

  /* Initialise the error structure */
  INIT_FAIL(fail);

  printf(
          "nag_rand_orthog_matrix (g05pxc) Example Program Results\n\n");

  /* Get the length of the state array */
  lstate = -1;
  nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  pda = n;
  a_size = pda * m;

  /* Allocate arrays */
  if (!(a = NAG_ALLOC(a_size, double)) ||
      !(state = NAG_ALLOC(lstate, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Initialise the generator to a repeatable sequence */
  nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Generate the random orthogonal matrix*/
  nag_rand_orthog_matrix(side, init, m, n, state, a, pda, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_rand_orthog_matrix (g05pxc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Display the results */
  for (i = 0; i < m; i++)
    {
      printf("  ");
      for (j = 0; j < n; j++)
        printf("%8.3f%s", A(i, j), (j+1)%4?" ":"\n");
      if (n%4) printf("\n");
    }

 END:
  NAG_FREE(a);
  NAG_FREE(state);

  return exit_status;
}