//      g02hb Example Program Text
//      C# version, NAG Copyright 2008
using System;
using NagLibrary;
using System.Globalization;
using System.IO;
namespace NagDotNetExamples
{
  public class G02HBE
  {
    static string datafile = "ExampleData/g02hbe.d";
    static void Main(String[] args)
    {
      if (args.Length == 1)
      {
        datafile = args[0];
      }
      StartExample();
    }
    public static void StartExample()
    {
      try
      {
        G02.G02HB_UCV ucvG02HB = new G02.G02HB_UCV(ucv);
        DataReader sr = new DataReader(datafile);
        double bd,   bl,   tol; int i,  j,  k,  l1,  l2,  m,  maxit,  mm,  n,  nit,  nitmon;
        int ifail;
        Console.WriteLine("g02hb Example Program Results");
        //      Skip heading in data file
        sr.Reset();
        //      Read in the dimensions of X
        sr.Reset();
        n = int.Parse(sr.Next());
        m = int.Parse(sr.Next());
        double[] a = new double[m*(m+1)/2];
        double[,] x = new double[n, m];
        double[] z = new double[n];
        if (n > 0 && m > 0 && m <= n)
        {
          //         Read in the x matrix
          for (i = 1; i <= n; i++)
          {
            sr.Reset();
            for (j = 1; j <= m; j++)
            {
              x[i - 1, j - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture);
            }
          }
          //         Read in the initial value of A
          mm = (m + 1) * m / 2;
          sr.Reset();
          for (j = 1; j <= mm; j++)
          {
            a[j - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture);
          }
          //         Set the values remaining parameters
          bl = 0.90e0;
          bd = 0.90e0;
          maxit = 50;
          tol = 0.50e-4;
          // Change nitmon to a positive value if monitoring information
          // is required
          nitmon = 0;
          // 
          G02.g02hb(ucvG02HB, n, m, x, a, z, bl, bd, tol, maxit, nitmon, out nit, out ifail);
          // 
          if (ifail == 0)
          {
            Console.Write("  {0}{1,4}{2}", "g02hb required ", nit, " iterations to converge");
            Console.WriteLine(" ");
            Console.Write(" {0}", "Matrix A");
            Console.WriteLine(" ");
            l2 = 0;
            for (j = 1; j <= m; j++)
            {
              l1 = l2 + 1;
              l2 = l2 + j;
              for (k = l1; k <= l2; k++)
              {
                Console.Write(" {0, 9:f4}", a[k - 1]);
              }
              Console.WriteLine(" ");
            }
            Console.WriteLine(" ");
            Console.WriteLine(" {0}", "Vector Z");
            for (i = 1; i <= n; i++)
            {
              Console.WriteLine("  {0,9:f4}", z[i - 1]);
            }
            //            Calculate Krasker-Welsch weights
            Console.WriteLine(" ");
            Console.WriteLine(" {0}", "Vector of weights");
            for (i = 1; i <= n; i++)
            {
              z[i - 1] = 1.00e0 / z[i - 1];
              Console.WriteLine("  {0,9:f4}", z[i - 1]);
            }
          }
          else
          {
            Console.WriteLine(" ");
            Console.WriteLine("** g02hb failed with ifail = {0,5}", ifail);
          }
          // 
        }
        // 
      }
      catch (Exception e)
      {
        Console.WriteLine(e.Message);
        Console.Write("Exception Raised");
      }
    }
    // 
    public static double ucv(double t)
    {
      const double ucvc = 2.50;
      double ucvValue = 0.0;
      // ucv function for Krasker-Welsch weights
      double pc,   pd,   q,   q2;
      ucvValue = 1.00e0;
      if (t != 0.00e0)
      {
        q = ucvc / t;
        q2 = q * q;
        pc = S.s15ab(q);
        if (q2 < -Math.Log(X02.x02ak()))
        {
          pd = Math.Exp(-q2 / 2.00e0) / Math.Sqrt(X01.x01aa() * 2.00e0);
        }
        else
        {
          pd = 0.00e0;
        }
        ucvValue = (2.00e0 * pc - 1.00e0) * (1.00e0 - q2) + q2 - 2.00e0 * q * pd;
      }
      return ucvValue;
    }
  }
}