//      c05rc Example Program Text
//      C# version, NAG Copyright 2008
using System;
using NagLibrary;
using System.Globalization;
namespace NagDotNetExamples
{
  public class C05RCE
  {
    static double[] ruser = new double[5];
    static void Main(String[] args)
    {
      StartExample();
    }
    public static void StartExample()
    {
      try
      {
        const int n = 9;
	      // C05.C05PC_FCN fcnC05PC = new C05.C05PC_FCN(fcn);
        C05.C05RC_FCN fcnC05RC = new C05.C05RC_FCN(fcn);
        double factor,  fnorm,  xtol = 0.0; int j, maxfev,  mode, nfev = 0, njev = 0, nprint = 0;
        double[] diag = new double[n];
        double[,] fjac = new double[n, n];
        double[] fvec = new double[n];
        double[] qtf = new double[n];
        double[] r = new double[n*(n+1)/2];
        double[] x = new double[n];
        int ifail;
        Console.WriteLine("c05rc Example Program Results");
        Console.WriteLine("");
        //      Store the coefficients describing the system
        ruser[0] = -1.00e0;
        ruser[1] = 3.00e0;
        ruser[2] = -2.00e0;
        ruser[3] = -2.00e0;
        ruser[4] = -1.00e0;
        //
        //      The following starting values provide a rough solution.
        for (j = 1; j <= n; j++)
        {
          x[j - 1] = -1.00e0;
        }
        xtol = Math.Sqrt(X02.x02aj());
        for (j = 1; j <= n; j++)
        {
          diag[j - 1] = 1.00e0;
        }
        maxfev = 1000;
        mode = 2;
        factor = 100.00e0;
        //
	/*C05.c05pc(fcnC05PC, n, x, fvec, fjac, xtol, maxfev, diag, mode, factor, nprint, out nfev,
        out njev, r, qtf, out ifail);*/
        C05.c05rc(fcnC05RC, n, x, fvec, fjac, xtol, maxfev, mode, diag, factor, nprint, out nfev,
        out njev, r, qtf, out ifail);
        //
        if (ifail == 0)
        {
          fnorm = F06.f06ej(n, fvec, 1, out ifail);
          Console.WriteLine("  {0}{1,12:e4}", "Final 2-norm of the residuals =", fnorm);
          Console.WriteLine("");
          Console.WriteLine("  {0}{1,10}", "Number of function evaluations =", nfev);
          Console.WriteLine("");
          Console.WriteLine("  {0}{1,10}", "Number of Jacobian evaluations =", njev);
          Console.WriteLine("");
          Console.WriteLine(" {0}", "Final approximate solution");
          Console.WriteLine("");
          for (j = 1; j <= n; j++)
          {
            Console.Write(" {0, 10:f4}{1}", x[j - 1], j%3==0?"\n":"");
          }
          Console.WriteLine("");
        }
        else if (ifail < 0)
        {
          Console.WriteLine("");
          Console.WriteLine("** c05rc failed with ifail = {0,5}", ifail);
        }
        else
        {
          Console.WriteLine("  {0}{1,5}", "ifail = ", ifail);
          if (ifail > 2)
          {
            Console.WriteLine("");
            Console.WriteLine(" {0}", "Approximate solution:");
            Console.WriteLine("");
            for (j = 1; j <= n; j++)
            {
              Console.Write(" {0, 10:f4}{1}", x[j - 1], j%3==0?"\n":"");
            }
            Console.WriteLine("");
          }
        }
        //
      }
      catch (Exception e)
      {
        Console.WriteLine(e.Message);
        Console.WriteLine("Exception Raised");
      }
      Console.WriteLine("");
    }
    //
    public static void fcn(int n, double[] x, double[] fvec, double[,] fjac,
    ref int iflag)
    {
      const double zero = 0.00e0;
      if (iflag == 0)
      {
        //
        //         Insert print statements here when nprint is positive.
        //
      }
      else if (iflag != 2)
      {
        for (int k = 1; k <= n; k++)
        {
          fvec[k - 1] = (ruser[1] + ruser[2] * x[k - 1]) * x[k - 1] - ruser[4];
          if (k > 1)
          {
            fvec[k - 1] = fvec[k - 1] + ruser[0] * x[k - 2];
          }
          if (k < n)
          {
            fvec[k - 1] = fvec[k - 1] + ruser[3] * x[k];
          }
        }
      }
      else
      {
        for (int k = 1; k <= n; k++)
        {
          for (int j = 1; j <= n; j++)
          {
            fjac[k - 1, j - 1] = zero;
          }
          fjac[k - 1, k - 1] = ruser[1] + 2.00e0 * ruser[2] * x[k - 1];
          if (k > 1)
          {
            fjac[k - 1, k - 2] = ruser[0];
          }
          if (k < n)
          {
            fjac[k - 1, k] = ruser[3];
          }
        }
      }
    }
  }
}