using System;
using System.Text;
using System.Runtime.InteropServices;/* Provides mappings between C# and native code */

namespace NumericalAlgorithmsGroup
{
  class LUSolve
    {   
      [StructLayout(LayoutKind.Sequential)]/* Force sequential mapping */
	public struct complex
	{
	  public double r;
	  public double i;
	};
      /*     F07APF Example Program Text */
      /*     Mark 21 Release. NAG Copyright 2004. */
      
      [DllImport("FLDLL214Z_nag.dll")]/* Import from DLL, the C# compiler provides a rudimentry check of the signature */
        public static extern void ZGESVX(string fact, int fact_len, string trans, int trans_len, ref int n, ref int nrhs, 
                                         complex [,]a, ref int lda, complex [,]af, ref int ldaf, int []ipiv, StringBuilder
                                         equed, int len_equed, double []r, double [] c, complex [,] b, ref int ldb, 
                                         double [,,] x, ref int ldx, ref double rcond, double [] ferr, double [] berr,
                                         complex [] work, double [] rwork, ref int info);
      [DllImport("FLDLL214Z_nag.dll")]
        public static extern void X04DAF(string matrix, int matrix_len, string diag, int diag_len, ref int m, ref int n, 
                                         double [,,] a, ref int lda, string title, int title_len, ref int ifail);
      public static void Main()
        {
          /* Scalars */
          double rcond=0.0;
          int i, ifail, info=99, j, k, n, nrhs;
          
          /* Character args */
          /*string equed="X";*/
          StringBuilder equed = new StringBuilder("X"); /* StringBuilder type used as equed is both input and output.
							   String type is immutable.
							*/

							   
          string complexadata ="-1.34,2.55,0.28,3.17,-6.39,-2.20,0.72,-0.92,-1.70,-14.10,33.10,-1.50,-1.50,13.40,12.90,13.80,-3.29,-2.39,-1.91,4.42,-0.14,-1.35,1.72,1.35,2.41,0.39,-0.56,1.47,-0.83,-0.69,-1.96,0.67";
          string complexbdata="26.26,51.78,31.32,-6.70,64.30,-86.80,158.60,-14.20,-5.75,25.31,-2.15,30.19,1.16,2.57,-2.56,7.55";
          Console.WriteLine("F07APF Example Program Results\n");
          n = 4;
          nrhs = 2;
          complex [,] a = new complex [n,n];
          complex [,] at =new complex [n,n];
          complex [,] af = new complex [n, n];
          complex [,] b = new complex [n,nrhs];
          complex [,] bt = new complex [nrhs,n];
          double [] berr = new double [nrhs];
          double [] c = new double [n];
          double [] ferr = new double [nrhs];
          double [] r = new double [n];
          complex [] work = new complex [4*n];
          double [,,] x = new double [nrhs,n,2]; /* Note Fortran is expecting complex array
	                                            This provides an equivalent double array type
						 */

          int [] ipiv = new int [n];
          double [] rwork = new double [2*n];

          /* Read data for a*/
          string [] complexdataarray = complexadata.Split(',');
          int index=0;
          for (i=0; i<n; i++)
            {
              for (j=0; j<n; j++) 
                {
                  a[i,j].r = double.Parse(complexdataarray[index++]);
                  a[i,j].i = double.Parse(complexdataarray[index++]);
                }
            }

          /* Print data */
	  /*          for (i=0; i<n; i++)
            {
              for (j=0; j<n; j++) 
                {
                  Console.Write("({0, 5:f2}, {1, 5:f2} ) ", a[i,j].r, a[i,j].i);
                }
	      Console.WriteLine();
            }
	  */
          /* Transpose A */
          for (i = 0; i < n; ++i) {
            for (j = 0; j < n; ++j) {
              at[j,i] = a[i,j]; 
            }
          }


          /* Read data for b*/
          complexdataarray = complexbdata.Split(',');
	  index=0;
          for (i=0; i<n; i++)
            {
              for (j=0; j<nrhs; j++) 
                {
                  b[i,j].r = double.Parse(complexdataarray[index++]);
                  b[i,j].i = double.Parse(complexdataarray[index++]);
                }
            }

          /* Print data */
	  /*
          for (i=0; i<n; i++)
            {
              for (j=0; j<nrhs; j++) 
                {
                  Console.Write("({0, 5:f2}, {1, 5:f2} ) ", b[i,j].r, b[i,j].i);
                }
	      Console.WriteLine();
            }

	  */
          /* Transpose b */
          for (i = 0; i < n; ++i) {
            for (j = 0; j < nrhs; ++j) {
              bt[j,i] = b[i,j]; /* bt_ref macros carries out the transform */
            }
          }

          /*        Solve the equations AX = B for X */

          ZGESVX("Equilibration", 13, "No transpose", 12, ref n, ref nrhs, at, ref n,
                 af, ref n, ipiv,  equed, 1, r, c, bt, ref n, x, ref n, ref rcond, ferr,
                 berr, work, rwork, ref info);
          if (info == 0 || info == n + 1) {

            /*  Print solution, error bounds, condition number, the form */
            /*  of equilibration and the pivot growth factor */

            /* x is returned from Fortran and so is 'transposed' already */
	    /* ifail = 0 */
	    /*  X04DAF("General", 7, " ", 1, ref n, ref nrhs, x, ref n, "Solution(s)", 11,
                ref ifail);
	    */
	    /* Print without using X04DAF. This shows how to map a transposed 2D double array into complex */
	    Console.WriteLine("   Solution 1         Solution 2");
	    for(i=0; i<n; i++)
	      {
		for(j=0;j<nrhs;j++)
		  {
		    Console.Write("(");
		    for (k=0; k<2;k++)
		      Console.Write(k==0?"{0, 5:f2}, ":"{0, 5:f2} ", x[j,i, k]); /* Printing a transposed double array as complex */
		    Console.Write(")   ");
		  }
		Console.WriteLine();
	      }
	  Console.WriteLine("\nBackward errors (machine-dependent");
	    
	  for (j = 1; j <= nrhs; ++j) {
	    Console.WriteLine("{0, 11:e1}{1}", berr[j - 1], j%7 == 0 || j == nrhs ?"\n":"");
	  }
	  Console.WriteLine();
            
	  Console.WriteLine("Estimated forward error bounds (machine-dependent)");
	  for (j = 1; j <= nrhs; ++j) {
	    Console.WriteLine("{0, 11:e1}{1}", ferr[j - 1], j%7 == 0 || j == nrhs ?"\n":"");
	  }
            
	  /*	  Console.WriteLine("equed = {0}", equed.ToString());*/
	  if (equed.ToString() == "N") {
	    Console.WriteLine("A has not been equilibrated");
	  } else if (equed.ToString() == "R") {
	    Console.WriteLine("A has been row scaled as diag(R)*A");
	  } else if (equed.ToString() == "C") {
	    Console.WriteLine("A has been column scaled as A*diag(C)");
	  } else if (equed.ToString() == "B") {
	    Console.WriteLine("A has been row and column scaled as diag(R)*A*diag(C)");
	  }
	  
	  Console.WriteLine("\nReciprocal condition number estimate of scaled matrix");
	  Console.WriteLine("{0, 11:e1}", rcond);
	  
	  Console.WriteLine("\nEstimate of reciprocal pivot growth factor");
	  Console.WriteLine("{0, 11:e1}",  rwork[0]);
	  if (info == n + 1) {
	    
	    Console.WriteLine("The matrix A is singular to working precision");
	  }
	} else {
	  Console.WriteLine("The ( {0}) element of the factor U is zero", info);
	}
          return;
    }
}
}