package edu.rice.linpack.Matrix.FMatrix;
import edu.rice.linpack.util.*;
import edu.rice.linpack.Vector.*;
import edu.rice.linpack.LNumber.*;

public class FSiFull extends FFull {

  public FSiFull() {
    super();
  }
  public FSiFull(int n, int j) {
    super(n,j);
  }
  public FSiFull(float[][] f) {
    super(f);
  }
  public FSiFull(FSiFull F) {
    super(F);
  }

  public FSiPack pack() {
    FSiPack S = new FSiPack(this.Mat);
    return S;
  }

  public void transpose() {
  }

  public float getElem(int i, int j) {
    if(i > j)
      return super.getElem(j,i);
    return super.getElem(i,j);
  }
  public float[] getRow(int i) {
    return getRow(i,0);
  }
  public float[] getColumn(int i) {
    return getRow(i,0);
  }
  public float[] getColumn(int i, int j) {
    return getRow(i,j);
  }
  public float[] getRow(int i, int j) {
    float[] F = new float[rows - j];
    for(int k=0;k<i && k<(rows-j);k++) {
      F[k] = Mat[j+k][i];
    }
    for(int k=i;k<(rows-j);k++) {
      F[k] = Mat[j+k][i];
    }
    return F;
  }
  public void setColumn(int c, float[] F) {
    for(int r=0;r<c;r++) 
      Mat[r][c] = F[r];
    for(int r=c;r<rows;r++) 
      Mat[r][c] = F[r];
  }
  public void setElem(int i, int c, float L) {
    if(i > c)
      super.setElem(c,i,L);
    else
      super.setElem(i,c,L);
  }

  protected float oneNorm() {
    float[] Z = new float[cols];
    
    for(int j=0;j<cols;j++) {
      Z[j] = this.asum(j+1,1,0,j);
      for(int i=0;i<j;i++) {
	Z[i] += Math.abs(this.Mat[i][j]);
      }
    }
    float anorm = 0;
    for(int j=0;j<cols;j++) 
      anorm = Math.max(anorm,Z[j]);
    return anorm;
  }


  /*
      inertia returns values different from matlab in sparse mats.
      */ 

  public void factor() 
       throws SingularMatrixException
  {    
    int info = 0;
    int kstep = 1;
    boolean swap;
    
    float alpha = (float)(1+Math.sqrt(17))/(float)8;

    for(int k=cols-1;k>=0;k-=kstep) {
      if (k == 0) {
	pivot[0] = 1;
	if(this.Mat[0][0] == 0) 
	  throw new SingularMatrixException(1);
      }
      else {
	int km = k - 1;
	
	float absakk = Math.abs(this.Mat[k][k]);
	int imax = this.i_amax(k,1,0,k);
	float colmax = Math.abs(this.Mat[imax][k]);
	
	if(absakk >= alpha*colmax) {
	  kstep = 1;
	  swap = false;
	}
	else {
	  float rowmax = 0;
	  for(int j=imax+1;j<=k;j++) {
	    rowmax = Math.max(rowmax,Math.abs(this.Mat[imax][j]));
	  }
	  if(imax != 0) {
	    int jmax = this.i_amax(imax,1,0,imax);
	    rowmax = Math.max(rowmax,Math.abs(this.Mat[jmax][imax]));
	  }
	  if(Math.abs(this.Mat[imax][imax]) >= alpha*rowmax) {
	    kstep = 1;
	    swap = true;
	  }
	  else {
	    if(absakk >= alpha*colmax*(colmax/rowmax)) {
	      kstep = 1;
	      swap = false;
	    }
	    else {
	      kstep = 2;
	      swap = (imax != km);
	    }
	  }
	}
	if(Math.max(absakk,colmax) == 0) {
	  pivot[k] = k+1;
	  info = k+1;
	}
	else {
	  if(kstep == 1) {
	    if(swap) {
	      this.swap(imax+1,1,0,imax,this,1,0,k);
	      for(int j=k;j>=imax;j--) 
		this.swapElem(j,k,imax,j);
	    }
	    for(int j=km;j>=0;j--) {
	      float mulk = -this.Mat[j][k]/this.Mat[k][k];
	      this.axpy(j+1,mulk,1,0,k,this,1,0,j);
	      this.Mat[j][k] = mulk;
	    }
	    if(swap) 
	      pivot[k] = imax+1;
	    else 
	      pivot[k] = k+1;
	  }
	  else {
	    if(swap) {
	      this.swap(imax+1,1,0,imax,this,1,0,km);
	      for(int j=km;j >= imax;j--) 
		this.swapElem(j,km,imax,j);
	      this.swapElem(km,k,imax,k);
	    }
	    if(k > 1) {
	      float ak = this.Mat[k][k]/this.Mat[km][k];
	      float akm = this.Mat[km][km]/this.Mat[km][k];
	      float dnom = 1 - ak*akm;
	      for(int j=km-1;j>=0;j--) 
		this.mulk(k,j,ak,akm,dnom);
	    }
	    if(swap) 
	      pivot[k] = -imax - 1;
	    else 
	      pivot[k] =  -k;
	    pivot[km] = pivot[k];
	  }
	}
      }
    }
    if(info != 0) 
      throw new SingularMatrixException(info);
  }
  private void mulk(int k, int j, float ak, float akm, float dnom) {
    int km = k-1;
    float Bk = this.Mat[j][k]/this.Mat[km][k];
    float Bkm = this.Mat[j][km]/this.Mat[km][k];
    float Mulk = (akm*Bk - Bkm)/dnom;
    float Mulkm = (ak*Bkm - Bk)/dnom;
    this.axpy(j+1,Mulk,1,0,k,this,1,0,j);
    this.axpy(j+1,Mulkm,1,0,km,this,1,0,j);
    this.Mat[j][k] = Mulk;
    this.Mat[j][km] = Mulkm;
  }
  
  public LNumber condition() 
       throws SingularMatrixException, WrongDataTypeException
  {
    Vector Z = new FVector(cols);
    return this.condition(Z);
  }
  public LNumber condition(Vector Ze) 
       throws SingularMatrixException, WrongDataTypeException
  {
    float[] Z = Ze.getFloatArray();
    float S;

    float anorm = this.oneNorm();
    float ynorm;

    //  Factor  //

    try {
      this.factor();
    } finally {

      //  Compute the inverse of the condition  //
      
      for(int j=0;j<cols;j++) 
	Z[j] = 0;
      
      this.solveUDW(Z);
      
      //  Solve trans(U)*Y = W  //
      
      this.solveTransUY(Z);

      S = 1/FUtil.asum(cols,Z,1);
      FUtil.scal(cols,S,Z,1);

      ynorm = this.solveUDV(Z);
    
      this.solveTransUY(Z);

      S = 1/FUtil.asum(cols,Z,1);
      FUtil.scal(cols,S,Z,1);
      ynorm *= S;
    }
    LFloat R;

    if(anorm != 0)
      R = new LFloat(ynorm/anorm);
    else
      R = new LFloat(anorm);

    return R;
  }
  private void solveTransUY(float[] Z) {

    int ks = 1;
    for(int k=0;k<cols;k+=ks) {
      ks = 1;
      if(pivot[k] < 0) 
	ks = 2;
      if(k != 0) {
	Z[k] += this.dot(k,1,0,k,Z,1,0);
	if(ks == 2) 
	  Z[k+1] += this.dot(k,1,0,k+1,Z,1,0);
	int kp = Math.abs(pivot[k]) - 1;
	if(kp != k) 
	  FUtil.swapElems(Z,k,kp);
      }
    }
  }
  private float solveUDV(float[] Z) {
    
    int ks = 1;
    float ynorm = 1;

    for(int k=cols-1;k>=0;k-=ks) {
      ks = 1;
      if(pivot[k] < 0) 
	ks = 2;
      if(k != ks-1) {
	int kp = Math.abs(pivot[k]) - 1;
	int kps = k - ks + 1;
	if(kp != kps) 
	  FUtil.swapElems(Z,kp,kps);
	this.axpy(k-ks+1,Z[k],1,0,k,Z,1,0);
	if(ks == 2) 
	  this.axpy(k-ks+1,Z[k-1],1,0,k-1,Z,1,0);
      }
      if(ks == 1) {
	if(Math.abs(Z[k]) > Math.abs(this.Mat[k][k])) {
	  float S = Math.abs(this.Mat[k][k])/Math.abs(Z[k]);
	  FUtil.scal(cols,S,Z,1);
	  ynorm *= S;
	}
	if(this.Mat[k][k] == 0) 
	  Z[k] = 1;
	else 
	  Z[k] = Z[k]/this.Mat[k][k];
      }
      else 
	this.Zfixer(k, Z);
    }
    float S = 1/FUtil.asum(cols,Z,1);
    FUtil.scal(cols,S,Z,1);
    ynorm *= S;
    
    return ynorm;
  }
  private void solveUDW(float[] Z) {
    
    int ks = 1;
    float ek = 1;
    for(int k=cols-1;k>=0;k-=ks) {
      ks = 1;
      if(pivot[k] < 0) 
	ks = 2;
      int kp = Math.abs(pivot[k]) - 1;
      int kps = k + 1 - ks;
      if(kp != kps) 
	FUtil.swapElems(Z,kp,kps);
      if(Z[k] != 0) 
	ek = FUtil.signOf(ek,Z[k]);
      Z[k] += ek;
      this.axpy(k-ks+1,Z[k],1,0,k,Z,1,0);
      if(ks == 1) {
	if(Math.abs(Z[k]) > Math.abs(this.Mat[k][k])) {
	  float S = Math.abs(this.Mat[k][k])/Math.abs(Z[k]);
	  FUtil.scal(cols,S,Z,1);
	  ek *= S;
	}
	if(this.Mat[k][k] == 0) 
	  Z[k] = 1;
	else 
	  Z[k] = Z[k]/this.Mat[k][k];
      }
      else {
	if(Z[k-1] != 0) 
	  ek = FUtil.signOf(ek,Z[k-1]);
	Z[k-1] += ek;
	axpy(k-ks+1, Z[k-1], 1,0,k-1,Z,1,0);
	this.Zfixer(k,Z);
      }
    }
    float S = 1/FUtil.asum(cols,Z,1);
    FUtil.scal(cols,S,Z,1);
  }
  private void Zfixer(int k, float[] Z) {
    int km = k-1;
    float ak = this.Mat[k][k]/this.Mat[km][k];
    float akm = this.Mat[km][km]/this.Mat[km][k];
    float bk = Z[k]/this.Mat[km][k];
    float bkm = Z[km]/this.Mat[km][k];
    float denom = ak*akm - 1;	
    Z[k] = (akm*bk - bkm)/denom;
    Z[km] = (ak*bkm - bk)/denom;
  }

  public void solve(Vector B, int Job) 
       throws WrongDataTypeException
  {
    this.solve(B);
  }
  public void solve(Vector Be) 
       throws WrongDataTypeException
  {
    
    float[] B = Be.getFloatArray();

    //  Loop backward applying the transformations and D inv to B  //
    this.loopback(B);

    //  Loop forward applying the transformations  // 
    this.loopforward(B);
  }
  private void loopback(float[] B) {
    int kstep = 1;
    for(int k=cols-1;k>=0;k-=kstep) {
      if(pivot[k] >= 0) {
	if(k != 0) {
	  int kp = pivot[k] - 1;
	  if(kp != k) 
	    FUtil.swapElems(B,k,kp);
	  this.axpy(k,B[k],1,0,k,B,1,0);
	}
	B[k] = B[k]/this.Mat[k][k];
	kstep = 1;
      }
      else {
	int km = k-1;
	if(k != 1) {
	  int kp = Math.abs(pivot[k]) - 1;
	  if(kp != km) 
	    FUtil.swapElems(B,km,kp);
	  this.axpy(km,B[k],1,0,k,B,1,0);
	  this.axpy(km,B[km],1,0,km,B,1,0);
	}
	this.Zfixer(k,B);
	kstep = 2;
      }
    }
  }
  private void loopforward(float[] B) {
    int kstep = 1;
    for(int k=0;k<cols;k+=kstep) {
      if(pivot[k] >= 0) 
	kstep = 1;
      else
	kstep = 2;
      if(k != 0) {
	B[k] += this.dot(k,1,0,k,B,1,0);
	if(kstep == 2)
	  B[k+1] += this.dot(k,1,0,k+1,B,1,0);
	int kp = Math.abs(pivot[k]) - 1;
	if(kp != k) 
	  FUtil.swapElems(B,k,kp);
      }
    }
  }

  public Vector determ() {
    
    float[] Det = new float[2];
    
    Det[0] = 1;
    Det[1] = 0;
    
    TD G = new TD();
    for(int k=0;k<cols;k++) {

      this.TDer(G,k);
      Det[0] *= G.D;
      if(Det[0] != 0) 
	FUtil.detNorm(Det);
    }
    FVector V = new FVector(Det);
    return V;
  }

  public int[] inertia() {

    int[] iner = new int[3];
    
    iner[0] = 0;
    iner[1] = 0;
    iner[2] = 0;
    
    TD A = new TD();
    
    for(int k=0;k<cols;k++) {
      this.TDer(A,k);

      if(A.D > 0) 
	iner[0] += 1;
      else if(A.D < 0) 
	iner[1] += 1;
      else 
	iner[2] += 1;
    } 
    return iner;
  }
  private void TDer(TD A, int k) {
    A.D = this.Mat[k][k];
      
    if(pivot[k] < 0) {
      if(A.T == 0) {
	A.T = Math.abs(this.Mat[k][k+1]);
	A.D = (A.D/A.T)*this.Mat[k+1][k+1] - A.T;
      }
      else {
	A.D = A.T;
	A.T = 0;
      }
    }
  }

  public void inverse() {
    
    float[] Work = new float[cols];
    int kstep;

    for(int k=0;k<cols;k+=kstep) {
      int kp = k + 1;
      if(pivot[k] >= 0) {
	this.Mat[k][k] = 1/this.Mat[k][k];
	
	if(k > 0) {
	  this.copy(k,1,0,k,Work,1);
	  for(int j=0;j<k;j++) {
	    this.Mat[j][k] = this.dot(j+1,1,0,j,Work,1,0);
	    this.axpy(j,Work[j],1,0,j,this,1,0,k);
	  }
	  this.Mat[k][k] += this.dot(k,1,0,k,Work,1,0);
	}
	kstep = 1;
      }
      else {
	this.invAdj(k,kp);
	if(k > 0) {
	  this.copy(k,1,0,kp,Work,1);
	  
	  for(int j=0;j<k;j++) {
	    this.Mat[j][kp] = this.dot(j+1,1,0,j,Work,1,0);
	    this.axpy(j,Work[j],1,0,j,this,1,0,kp);
	  }
	  this.Mat[kp][kp] += this.dot(k,1,0,kp,Work,1,0);
	  this.Mat[k][kp] += this.dot(k,1,0,k,this,1,0,kp);
	  this.copy(k,1,0,k,Work,1);
	  
	  for(int j=0;j<k;j++) {
	    this.Mat[j][k] = this.dot(j+1,1,0,j,Work,1,0);
	    this.axpy(j,Work[j],1,0,j,this,1,0,k);
	  }
	  this.Mat[k][k] += this.dot(k,1,0,k,Work,1,0);
	}
	kstep = 2;
      }
      int ks = Math.abs(pivot[k]) - 1;
      if(ks != kp) {
	this.swap(ks+1,1,0,ks,this,1,0,k);
	for(int j=k;j>=ks;j--) 
	  this.swapElem(j,k,ks,j);
	if(kstep != 1) 
	  this.swapElem(k,kp,ks,kp);
      }
    }
  }
  private void invAdj(int k, int kp) {
    float T = Math.abs(this.Mat[k][kp]);
    float ak = this.Mat[k][k]/T;
    float akp = this.Mat[kp][kp]/T;
    float akkp = this.Mat[k][kp]/T;
    float D = T*(ak*akp - 1);
    this.Mat[k][k] = akp/D;
    this.Mat[kp][kp] = ak/D;
    this.Mat[k][kp] = -akkp/D;
  }
}