//      f08bf Example Program Text
//      C# version, NAG Copyright 2008
using System;
using NagLibrary;
using System.Globalization;
namespace NagDotNetExamples
{
  public class F08BFE
  {
    const double one=1.00e0;
    const double zero=0.00e0;
    static string datafile = "ExampleData/f08bfe.d";
    static void Main(String[] args)
    {
      if (args.Length == 1)
      {
        datafile = args[0];
      }
      StartExample();
    }
    public static void StartExample()
    {
      try
      {
        DataReader sr = new DataReader(datafile);
        double tol=0.0;
        int i,   info,   j,   k,   m,   n,   nrhs; sr.Reset();
        sr.Reset();
        m = int.Parse(sr.Next());
        n = int.Parse(sr.Next());
        nrhs = int.Parse(sr.Next());
        double[,] a = new double[m, n];
        double[,] b = new double[m, nrhs];
        double[] rnorm = new double[n];
        double[] tau = new double[Math.Min(m,n)];
        double[] work = new double[n];
        int[] jpvt = new int[n];
        int ifail;
        Console.WriteLine("f08bf Example Program Results");
        Console.WriteLine("");
        //      Skip heading in data file
        double [] ArraySlice_b = new double [n];
        if (m > 0 && n > 0 && m >= n && nrhs > 0)
        {
          // 
          //         Read A and B from data file
          // 
          sr.Reset();
          for (i = 1 ; i <= m ; i++)
          {
            for (j = 1 ; j <= n ; j++)
            {
              a[i - 1 , j - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture);
            }
          }
          sr.Reset();
          for (i = 1 ; i <= m ; i++)
          {
            for (j = 1 ; j <= nrhs ; j++)
            {
              b[i - 1 , j - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture);
            }
          }
          // 
          //         Initialize jpvt to be zero so that all columns are free
          // 
          F06.f06db(n, 0, jpvt, 1, out ifail);

         if (ifail < 0)
          {
              return;
          }

          // 
          //         Compute the QR factorization of a
          // 
          F08.f08bf(m, n, a, jpvt, tau, out info);
          if (info < 0)
          {
              return;
          }

          // 
          //         Compute c = (c1) = (q**t)*b, storing the result in B
          //                     (c2)
          // 
          F08.f08ag("Left", "Transpose", m, nrhs, n, a, tau, b, out info);
          // 
          //         Choose tol to reflect the relative accuracy of the input data
          // 
          tol = 0.010e0;
          // 
          //         Determine and print the rank, k, of r relative to tol
          // 
          for (k = 1 ; k <= n ; k++)
          {
            if (Math.Abs(a[k - 1 , k - 1]) <= tol * Math.Abs(a[0 , 0]))
            {
              goto L40;
            }
          }
          L40: ;
          k = k - 1;
          // 
          Console.WriteLine(" {0}","Tolerance used to estimate the rank of A");
          Console.WriteLine("     {0,11:e2}",tol);
          Console.WriteLine(" {0}","Estimated rank of A");
          Console.WriteLine("  {0,8}",k);
          Console.WriteLine("");
          // 
          //         Compute least-squares solutions by backsubstitution in
          //         r(1:k,1:k)*y = c1, storing the result in b
          // 
          F06.f06yj("Left", "Upper", "No transpose", "Non-Unit", k, nrhs, one, a, b, out ifail);
          // 
          //         Compute estimates of the square roots of the residual sums of
          //         squares (2-norm of each of the columns of c2)
          // 
          for (j = 1 ; j <= nrhs ; j++)
          {
            for (int ind=0; ind < m-k; ind++)
            ArraySlice_b[ind]=b[k+ind,j-1];
            rnorm[j - 1] = F06.f06ej(m - k, ArraySlice_b, 1, out ifail);
          }
          // 
          //         Set the remaining elements of the solutions to zero (to give
          //         the basic solutions)
          // 
          double [, ] ArraySlice_b_2d = new double [n-k,  nrhs];
          F06.f06qh("General", n - k, nrhs, zero, zero, ArraySlice_b_2d, out ifail);
          for (i=0; i<n-k;i++)
          for (j=0; j <nrhs; j++)
          b[k+i, j] = ArraySlice_b_2d[i,j];
          // 
          //         Permute the least-squares solutions stored in B to give x = p*y
          // 
          for (j = 1 ; j <= nrhs ; j++)
          {
            for (i = 1 ; i <= n ; i++)
            {
              work[jpvt[i - 1] - 1] = b[i - 1 , j - 1];
            }
            F06.f06ef(n, work, 1, ArraySlice_b, 1, out ifail);
            for (int ind = 0; ind < n; ind++)
            b[ind, j-1] = ArraySlice_b[ind];
          }
          // 
          //         Print least-squares solutions
          // 
          X04.x04ca("General", "", n, nrhs, b, "Least-squares solution(s)", out ifail);
          // 
          //         Print the square roots of the residual sums of squares
          // 
          Console.WriteLine(" {0}","Square root(s) of the residual sum(s) of squares");
          for (j = 1 ; j <= nrhs ; j++)
          {
            Console.Write("      {0, 11:e2}", rnorm[j - 1]);
          }
          Console.WriteLine();
        }
        // 
      }
      catch (Exception e)
      {
        Console.WriteLine(e.Message);
        Console.WriteLine("Exception Raised");
      }
    }
  }
}