/// <summary>Performs the QR factorization.</summary>
		protected override void InternalCompute()
		{
			int m = matrix.Rows;
			int n = matrix.Columns;

#if MANAGED
			int minmn = m < n ? m : n;
			r_ = new DoubleMatrix(matrix); // create a copy
			DoubleVector[] u = new DoubleVector[minmn];
			for (int i = 0; i < minmn; i++)
			{
				u[i] = Householder.GenerateColumn(r_, i, m - 1, i);
				Householder.UA(u[i], r_, i, m - 1, i + 1, n - 1);
			}
			q_ = DoubleMatrix.CreateIdentity(m);
			for (int i = minmn - 1; i >= 0; i--)
			{
				Householder.UA(u[i], q_, i, m - 1, i, m - 1);
			}
#else
      qr = new double[matrix.data.Length];
      Array.Copy(matrix.data, qr, matrix.data.Length);
      jpvt = new int[n];
      jpvt[0] = 1;
      Lapack.Geqp3.Compute(m, n, qr, m, jpvt, out tau);
      r_ = new DoubleMatrix(m, n);
      // Populate R
      for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
          if (i <= j) {
            r_.data[j * m + i] = qr[(jpvt[j]-1) * m + i];
          }
          else {
            r_.data[j * m + i] = 0.0;
          }
        }
      }
      q_ = new DoubleMatrix(m, m);
      for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
          if (j < n)
            q_.data[j * m + i] = qr[j * m + i];
          else
            q_.data[j * m + i] = 0.0;
        }
      }

      if( m < n ){
        Lapack.Orgqr.Compute(m, m, m, q_.data, m, tau);
      } else{
        Lapack.Orgqr.Compute(m, m, n, q_.data, m, tau);
      }
#endif
			for (int i = 0; i < m; i++)
			{
				if (q_[i, i] == 0)
					isFullRank = false;
			}
		}

 public void CtorInitialValues()
 {
   DoubleVector test = new DoubleVector(2,1);
   
   Assert.AreEqual(test.Length, 2);
   Assert.AreEqual(test[0],1);
   Assert.AreEqual(test[1],1);
 }

 public void CtorDimensions()
 {
   DoubleVector test = new DoubleVector(2);
   
   Assert.AreEqual(test.Length, 2);
   Assert.AreEqual(test[0],0);
   Assert.AreEqual(test[1],0);
 }

 public override double Value (DoubleVector x) 
 {
   double retvalue=0;
   for (int i=1; i<x.Length; i++) 
   {
     retvalue = retvalue + 100*System.Math.Pow((x[i] - System.Math.Pow(x[i-1],2)),2) + System.Math.Pow((1-x[i-1]),2);
   }
   return retvalue;
 }

		public void CurrentException2()
		{
			DoubleVector test = new DoubleVector(new double[2] { 1, 2 });
			IEnumerator enumerator = test.GetEnumerator();
			enumerator.MoveNext();
			enumerator.MoveNext();
			enumerator.MoveNext();
			object value = enumerator.Current;
		}

 public void CtorArray()
 {
   double[] testvector = new double[2]{0,1};
   
   DoubleVector test = new DoubleVector(testvector);
   Assert.AreEqual(test.Length,testvector.Length);
   Assert.AreEqual(test[0],testvector[0]);
   Assert.AreEqual(test[1],testvector[1]);
 }

 public override DoubleVector Gradient(DoubleVector x) 
 {
   DoubleVector retvalue = new DoubleVector(x.Length,0.0);
   retvalue[0] = -400*x[0]*(x[1]-System.Math.Pow(x[0],2))-2*(1-x[0]);
   retvalue[x.Length-1] = 200*(x[x.Length-1]-System.Math.Pow(x[x.Length-2],2));
   if (x.Length>2) 
   {
     for (int i=1; i<x.Length-1; i++)
       retvalue[i] = 200*(x[i]-System.Math.Pow(x[i-1],2))-400*x[i]*(x[i+1]-System.Math.Pow(x[i],2))-2*(1-x[i]);
   }
   return retvalue;
 }

		public double Update(DoubleVector solution, DoubleVector direction, double beta)
		{
			DoubleVector newSolution;
			double newbeta = beta;
			for (int i = 0; i < 200; i++)
			{
				newSolution = solution + newbeta * direction;
				if (Check(newSolution))
				{
					return newbeta;
				}
				newbeta *= 0.5;
			}
			throw new OptimizationException("Beta couldn't be found to satisfy constraint");
		}

 public override DoubleMatrix Hessian(DoubleVector x) 
 {
   DoubleMatrix ret = new DoubleMatrix(x.Length,x.Length,0.0);
   
   for (int i=0; i<x.Length-1; i++)
   {
     ret[i,i+1] = -400*x[i];
     ret[i+1,i] = -400*x[i];
   }
   ret[0,0] = System.Math.Pow(1200*x[0],2)-400*x[1]+2;
   ret[x.Length-1,x.Length-1] = 200;
   for (int i=1; i<x.Length-1; i++)
     ret[i,i] = 202 + System.Math.Pow(1200*x[i],2) - 400*x[i+1];
   return ret;
 } 

		///<summary> Minimize the given cost function </summary>
		public override DoubleVector Search(DoubleVector x, DoubleVector d, double stp)
		{
			DoubleVector ret = new DoubleVector(x);
			double j = 0;
			double delta_d = d.GetDotProduct(d);
			double alpha;
			do
			{
				alpha = -GradientEvaluation(ret).GetDotProduct(d) /
					d.GetDotProduct(HessianEvaluation(ret) * d);
				ret = ret + alpha * d;
				j++;
			} while ((j < maxIteration) && (alpha * alpha * delta_d > tolerance * tolerance));
			return ret;
		}

		public void Current()
		{
			DoubleVector test = new DoubleVector(new double[2] { 1, 2 });
			IEnumerator enumerator = test.GetEnumerator();
			bool movenextresult;

			movenextresult = enumerator.MoveNext();
			Assert.IsTrue(movenextresult);
			Assert.AreEqual(enumerator.Current, test[0]);

			movenextresult = enumerator.MoveNext();
			Assert.IsTrue(movenextresult);
			Assert.AreEqual(enumerator.Current, test[1]);

			movenextresult = enumerator.MoveNext();
			Assert.IsFalse(movenextresult);
		}

 ///<summary>Method to override to calculate the grad_f, the first derivative of 
 /// the cost function with respect to x</summary>
 public virtual DoubleVector Gradient(DoubleVector x)
 {
   double eps = 1e-8;
   double fp, fm;
   DoubleVector grad = new DoubleVector(x.Length,0.0);
   
   DoubleVector xx = new DoubleVector(x);
   for (int i=0; i<x.Length; i++) 
   {
     xx[i] += eps;
     fp = this.Value(xx);
     xx[i] -= 2.0*eps;
     fm = this.Value(xx);
     grad[i] = 0.5*(fp - fm)/eps;
     xx[i] = x[i];
   }
   return grad;
 }

		public void TestRosenbrock()
		{
			Rosenbrock cf = new Rosenbrock();
			EndCriteria ec = new EndCriteria();
			ConjugateGradient optim = new ConjugateGradient(cf, ec);
			//  new SecantLineSearch(cf,ec));

			DoubleVector x0 = new DoubleVector(new double[5] { 1.3, 0.7, 0.8, 1.9, 1.2 });

			optim.Minimize(x0);

			Assert.AreEqual(optim.SolutionValue, 0.0, 0.1);
			Assert.AreEqual(optim.SolutionVector[0], 1.0, 0.1);
			Assert.AreEqual(optim.SolutionVector[1], 1.0, 0.1);
			Assert.AreEqual(optim.SolutionVector[2], 1.0, 0.1);
			Assert.AreEqual(optim.SolutionVector[3], 1.0, 0.2);
			Assert.AreEqual(optim.SolutionVector[4], 1.0, 0.4);
		}

 ///<summary> Minimize the given cost function </summary>
 public override DoubleVector Search(DoubleVector x, DoubleVector d, double step) 
 {
   DoubleVector ret = new DoubleVector(x);
   double j=0;
   double eta;
 
   double delta_d = d.GetDotProduct(d);
   double alpha = -sigma_0;
   double eta_prev = d.GetDotProduct(GradientEvaluation(ret+sigma_0*d));
   do 
   {
     eta = d.GetDotProduct(GradientEvaluation(ret));
     alpha = alpha*(eta/(eta_prev-eta));
     ret = ret + alpha*d;
     eta_prev = eta;
     j++;
   } while ((j<maxIteration) && (alpha*alpha*delta_d > tolerance*tolerance));
   return ret;
 }

 public void TestReflection() 
 {
   Poly cf = new Poly();
   NelderMead optim = new NelderMead(cf);
   
   DoubleVector[] simplex = new DoubleVector[3];
   simplex[0] = new DoubleVector(new double[2]{1,1});
   simplex[1] = new DoubleVector(new double[2]{1,-1});
   simplex[2] = new DoubleVector(new double[2]{2,0});
   
   optim.Rho = 1.5;
   optim.InitializeMethod(simplex);
   optim.IterateMethod();
   
   DoubleVector xr = (1+optim.Rho)*(new DoubleVector(new double[2]{1,0})) - optim.Rho*simplex[2];
   
   Assert.IsTrue(optim.LastStep == NelderMead.Step.Reflection);
   Assert.AreEqual(optim.Simplex[0][0],xr[0]);
   Assert.AreEqual(optim.Simplex[0][1],xr[1]);
 }

		public override LinearAlgebra.DoubleVector Search(LinearAlgebra.DoubleVector x, LinearAlgebra.DoubleVector direction, double step)
		{
			DoubleVector retx = new DoubleVector(x);
			double oldVal = FunctionEvaluation(retx);
			double newVal = oldVal;

			// First find the initial direction
			double valPos = FunctionEvaluation(retx + direction * step);
			double valNeg = FunctionEvaluation(retx - direction * step);
			if (valPos >= oldVal && valNeg < oldVal) // we reverse the direction only if the other direction really gives the smaller result
			{
				retx -= direction * step;
				oldVal = valNeg;
				step = -step;
			}
			else if (valPos < oldVal)
			{
				retx += direction * step;
				oldVal = valPos;
			}

			// now iterate
			for (; ; )
			{
				retx += direction * step;
				newVal = FunctionEvaluation(retx);

				if (newVal > oldVal)
				{
					step /= -2;
				}
				else if (!(newVal != oldVal))
				{
					break;
				}
				oldVal = newVal;
			}

			return retx;
		}

 public void TestInitializeMethod()
 {
   Rosenbrock cf = new Rosenbrock();
   NelderMead optim = new NelderMead(cf);
   DoubleVector x0 = new DoubleVector(new double[4]{0,1,2,3});
   
   optim.SimplexDelta = 0.1;
   optim.SimplexZeroDelta = 0.0001;
   
   optim.InitializeMethod(x0);
   
   Assert.AreEqual(optim.Simplex.Length,5);
   for (int i=0; i<optim.Simplex.Length; i++) 
   {
     Assert.AreEqual(optim.Simplex[i][0],x0[0],optim.SimplexZeroDelta);
     Assert.AreEqual(optim.Simplex[i][1],x0[1],optim.SimplexDelta*x0[1]+0.001);
     Assert.AreEqual(optim.Simplex[i][2],x0[2],optim.SimplexDelta*x0[2]+0.001);
     Assert.AreEqual(optim.Simplex[i][3],x0[3],optim.SimplexDelta*x0[3]+0.001);
   }
   for (int i=1; i<optim.Simplex.Length; i++) 
   {
     Assert.IsTrue(cf.Value(optim.Simplex[i-1])<cf.Value(optim.Simplex[i]));
   }
 }

 public void SetDiagonal()
 {
   DoubleMatrix a = new DoubleMatrix(2,2);
   DoubleVector b = new DoubleVector(2);
   b[0] = 1;
   b[1] = 2;
   a.SetDiagonal(b);
   Assert.AreEqual(b[0], a[0,0]);
   Assert.AreEqual(b[1], a[1,1]);
 }

 ///<summary> Constructor </summary>
 public DoubleVectorEnumerator (DoubleVector vector) 
 {
   v=vector;
   index=-1;
   length=v.Length;
 }

 public void SetColumnOutOfRange()
 {
   DoubleMatrix a = new DoubleMatrix(2,2);
   DoubleVector b = new DoubleVector(2);
   a.SetColumn(2,b);
 }