package dmrgvisu;

import org.netlib.lapack.DSYEV;
import org.netlib.util.*;

public class DMRGsystem
{
	block b1;
	block b2;
	block b3;
	block b4;
	public DMRGsystem(block bb1, block bb2, block bb3, block bb4)
	{
		b1 = bb1; b2 = bb2; b3 = bb3; b4 = bb4;
	}
	
	public double GetGroundState(WaveFunction p, DMRG DMRGref)
	{
		int dim = 4;
		double [][]H= new double [dim][dim];
		for (int i=0;i<dim;i++)
			for (int j=0; j<dim; j++)
				H[i][j]=0.0;
		H[0][0] = b1.H11;
		if (DMRGref.sweep%2 == 1)
		{
			H[1][1] = b2.H11+DMRGref.V(DMRGref.i-2);
			H[2][2] = b3.H11+DMRGref.V(DMRGref.i-1);
		}
		else
		{
			H[1][1] = b2.H11+DMRGref.V(DMRGref.i+1);
			H[2][2] = b3.H11+DMRGref.V(DMRGref.i+2);
		}
		
		H[3][3] = b4.H11;
		
		H[0][1] = -b1.L_inner;
		H[1][0] = -b1.L_inner;
		
		H[1][2] = -1.0/Math.pow(DMRGref.h, 2.0);
		H[2][1] = -1.0/Math.pow(DMRGref.h, 2.0);
		
		H[2][3] = -b4.L_inner;
		H[3][2] = -b4.L_inner;
		
		double []eva =new double[dim];
	    int lwork = dim*dim;
	    double []work = new double[lwork];
	    intW info = new intW(0);
		DSYEV.DSYEV("V", "U", dim, H, eva, work, lwork, info);
		for (int i = 0; i < dim; i++)
			p.V[i]=H[i][0];
		if (p.V[0] < 0.0 || p.V[1] < 0.0 || p.V[2] < 0.0 || p.V[3] < 0.0)
		{
			p.V[0] *= -1;
			p.V[1] *= -1;
			p.V[2] *= -1;
			p.V[3] *= -1;
		}
		return eva[0];
	}
}
