// e04yb Example Program Text // C# version, NAG Copyright 2008 using System; using NagLibrary; using System.Globalization; using System.IO; namespace NagDotNetExamples { public class E04YBE { static int mdec = 15; static int ndec = 3; public static double[,] t = new double[mdec, ndec]; public static double[] y = new double[mdec]; static string datafile = "ExampleData/e04ybe.d"; // as a command line argument. It defaults to the named file specified below otherwise. // static void Main(String[] args) { if (args.Length == 1) { datafile = args[0]; } StartExample(); } public static void StartExample() { E04.E04YA_LSQFUN lsqfunE04YA = new E04.E04YA_LSQFUN(lsqfun); E04.E04YB_LSQFUN lsqfunE04YB = new E04.E04YB_LSQFUN(lsqfun); E04.E04YB_LSQHES lsqhesE04YB = new E04.E04YB_LSQHES(lsqhes); try { DataReader sr = new DataReader(datafile); int i, ifail, j, k, m, n; Console.WriteLine("e04yb Example Program Results"); // Skip heading in data file sr.Reset(); m = mdec; n = ndec; int lb = n * (n + 1) / 2; double[] b = new double[lb]; double[,] fjac = new double[m, n]; double[] fvec = new double[m]; double[] x = new double[n]; int[] iw = new int[1]; // Observations of tj (j = 1, 2, 3) are held in t[i-1, j-1] //(i = 1, 2, . . . , 15) for (i = 1; i <= m; i++) { sr.Reset(); y[i - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture); for (j = 1; j <= n; j++) { t[i - 1, j - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture); } } // Set up an arbitrary point at which to check the derivatives x[0] = 0.190e0; x[1] = -1.340e0; x[2] = 0.880e0; // Check the 1st derivatives E04.e04ya(m, n, lsqfunE04YA, x, fvec, fjac, out ifail); if (ifail != 0) { Console.WriteLine(" ** e04ya returned with ifail = {0, 3}", ifail); goto L60; } // Console.WriteLine(""); Console.WriteLine(" {0}", "The test point is"); for (j = 1; j <= n; j++) { Console.Write(" " + " {0, 10:f5}", x[j - 1]); } Console.WriteLine(" "); // // Check the evaluation of B E04.e04yb(m, n, lsqfunE04YB, lsqhesE04YB, x, fvec, fjac, b, out ifail); // Console.WriteLine(""); if (ifail == 1) { Console.WriteLine(" {0}", "A parameter is outside its expected range"); } else if (ifail >= 0) { if (ifail == 0) { Console.WriteLine(" {0}", "The matrix B is consistent with 1st derivatives"); } else if (ifail == 2) { Console.WriteLine(" {0}", "Probable error in calculation of the matrix B"); } Console.WriteLine(""); Console.WriteLine(" {0}", "At the test point, LSQFUN gives"); Console.WriteLine(""); Console.WriteLine("{0}\n", " Residuals 1st derivatives"); Console.WriteLine(""); for (i = 1; i <= m; i++) { Console.Write(" " + "{0,15:e3}", fvec[i - 1]); for (j = 1; j <= n; j++) { Console.Write(" " + "{0,15:e3}", fjac[i - 1, j - 1]); } Console.WriteLine(" "); } Console.WriteLine(""); Console.WriteLine(" {0}", "and LSQHES gives the lower triangle of the matrix B"); Console.WriteLine(""); k = 1; for (i = 1; i <= n; i++) { for (j = k; j <= k + i - 1; j++) { Console.Write(" " + " {0, 15:e3}", b[j - 1]); } Console.WriteLine(" "); k = k + i; } } else if (ifail == -399) { Console.WriteLine(" ** e04yb returned with ifail = {0, 3}", ifail); } else { Console.WriteLine(" {0}{1,3}{2}", "IFLAG was set to ", ifail, "in LSQFUN or LSQHES"); } L60: ; // } catch (Exception e) { Console.WriteLine(e.Message); Console.WriteLine( "Exception Raised"); } } // public static void lsqfun(ref int iflag, int m, int n, double[] xc, double[] fvec, double[,] fjac) { // Routine to evaluate the residuals and their 1st derivatives double denom, dummy; int i = 0; for (i = 1; i <= m; i++) { denom = xc[1] * t[i - 1, 1] + xc[2] * t[i - 1, 2]; fvec[i - 1] = xc[0] + t[i - 1, 0] / denom - y[i - 1]; fjac[i - 1, 0] = 1.00e0; dummy = -1.00e0 / (denom * denom); fjac[i - 1, 1] = t[i - 1, 0] * t[i - 1, 1] * dummy; fjac[i - 1, 2] = t[i - 1, 0] * t[i - 1, 2] * dummy; } } // public static void lsqhes(ref int iflag, int m, int n, double[] fvec, double[] xc, double[] b) { // Routine to compute the lower triangle of the matrix B // (stored by rows in the array B) double dummy, sum22, sum32, sum33; int i = 0; b[0] = 0.00e0; b[1] = 0.00e0; sum22 = 0.00e0; sum32 = 0.00e0; sum33 = 0.00e0; for (i = 1; i <= m; i++) { dummy = 2.00e0 * t[i - 1, 0] / Math.Pow((xc[1] * t[i - 1, 1] + xc[2] * t[i - 1, 2]), 3); sum22 = sum22 + fvec[i - 1] * dummy * (t[i - 1, 1]) * (t[i - 1, 1]); sum32 = sum32 + fvec[i - 1] * dummy * t[i - 1, 1] * t[i - 1, 2]; sum33 = sum33 + fvec[i - 1] * dummy * (t[i - 1, 2]) * (t[i - 1, 2]); } b[2] = sum22; b[3] = 0.00e0; b[4] = sum32; b[5] = sum33; } } }