/* nag_real_cholesky_skyline_solve(f04mcc) Example Program
 *
 * Copyright 1996 Numerical Algorithms Group.
 *
 * Mark 4, 1996.
 */

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

#define NMAX 6
#define NRHSMAX 2
#define TDB NRHSMAX
#define TDX NRHSMAX
#define LALMAX 14

int main(void)
{
  Integer i, nrhs, k, k1, k2, lal, n;
  double a[LALMAX], al[LALMAX], b[NMAX][TDB], d[NMAX], x[NMAX][TDX];
  Integer row[NMAX];
  Nag_SolveSystem select;
  static NagError fail;

  Vprintf("f04mcc Example Program Results\n");
  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  Vscanf("%ld",&n);
  if (n<1 || n>NMAX)
    {
      Vprintf("\n n is out of range: n = %ld\n", n);
      return EXIT_FAILURE;
    }
  for (i=0; i<n; ++i)
    Vscanf("%ld",&row[i]);
  k2 = 0;
  for (i=0; i<n; ++i)
    {
      k1 = k2;
      k2 = k2 + row[i];
      for (k=k1; k<k2; ++k)
        Vscanf("%lf",&a[k]);
    }
  lal = k2;
  if (lal > LALMAX)
    {
      Vprintf("\n lal is out of range: lal = %ld\n", lal);
      return EXIT_FAILURE;
    }
  Vscanf("%ld",&nrhs);
  if (nrhs<1 || nrhs>NRHSMAX)
    {
      Vprintf("\n nrhs is out of range: nrhs = %ld\n", nrhs);
      return EXIT_FAILURE;
    }
  for (i=0; i<n; ++i)
    for (k=0; k<nrhs; ++k)
      Vscanf("%lf",&b[i][k]);
  fail.print = TRUE;
  f01mcc(n, a, lal, row, al, d, &fail);
  if (fail.code != NE_NOERROR)
    return EXIT_FAILURE;
  select = Nag_LDLTX;
  f04mcc(select, n, nrhs, al, lal, d, row, &b[0][0], (Integer)TDB,
         &x[0][0], (Integer)TDX, &fail);
  if (fail.code != NE_NOERROR)
    return EXIT_FAILURE;
  Vprintf("\n Solution\n");
  for (i=0; i<n; ++i)
    {
      for (k=0; k<nrhs; ++k)
        Vprintf("%9.3f",x[i][k]);
      Vprintf("\n");
    }
  return EXIT_SUCCESS;
}