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

public class FTriLow extends FTriFull {

  public FTriLow() {
    super();
  }
  public FTriLow(int i, int j) {
    Mat = new float[i][j];
    rows = i;
    cols = j;
  }
  public FTriLow(float[][] F) {
    rows = F.length;
    cols = (F[0]).length;
    Mat = F;
  }
  public FTriLow(FTriLow M) {
    rows = M.rows;
    cols = M.cols;
    Mat = new float[rows][cols];
    for(int i=0;i<rows;i++)
      for(int j=0;j<cols;j++) 
	Mat[i][j] = M.Mat[i][j];
  }

  protected float oneNorm() {
    float tnorm = 0;
    for(int j=0;j<cols;j++) {
      int q = cols-j;
      tnorm = Math.max(tnorm,this.asum(q,1,j,j));
    }
    return tnorm;
  }

  public LNumber condition(Vector Ze) 
       throws SingularMatrixException, WrongDataTypeException
  {
 
    float[] Z = Ze.getFloatArray();

    float S;
    float tnorm = this.oneNorm();
    float ynorm;

    try {
      this.zeroDiag();
    } finally {
      
      float ek = 1;
      for(int j=0;j<cols;j++) 
	Z[j] = 0;
      
      for(int k=cols-1;k>=0;k--) {
	if(Z[k] != 0) 
	  ek = FUtil.signOf(ek,-Z[k]); 
	
	if(Math.abs(ek-Z[k]) > Math.abs(this.Mat[k][k])) {
	  S = Math.abs(this.Mat[k][k])/Math.abs(ek-Z[k]);
	  FUtil.scal(cols,S,Z,1);
	  ek *= S;
	}
	float wk = ek - Z[k];
	float wkm = -ek - Z[k];
	S = Math.abs(wk);
	float SM = Math.abs(wkm);
      
	if(this.Mat[k][k] != 0) {
	  wk = wk/this.Mat[k][k];
	  wkm = wkm/this.Mat[k][k];
	}
	else {
	  wk = 1;
	  wkm = 1;
	}
	
	if(k > 0) {
	  for(int j=0;j<k;j++) {
	    SM += Math.abs(Z[j] + wkm*this.Mat[k][j]);
	    Z[j] += wk*this.Mat[k][j];
	    S += Math.abs(Z[j]);
	  }
	  if(S < SM) {
	    float T = wkm - wk;
	    wk = wkm;
	    for(int j=0;j<k;j++) 
	      Z[j] += T*this.Mat[k][j];
	  }
      }
	Z[k] = wk;
      }
      S = 1/FUtil.asum(cols,Z,1);
      FUtil.scal(cols,S,Z,1);

      ynorm = this.solveTZ(Z);
    }
    LFloat R;
    if(tnorm != 0) 
      R = new LFloat(ynorm/tnorm);
    else
      R = new LFloat(tnorm);
    return R;
  }
  private float solveTZ(float[] Z) {
    float ynorm = 1;
    
    for(int k=0;k<cols;k++) {
      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] = Z[k]/this.Mat[k][k];
      else 
	Z[k] = 1;
      
      int q = k+1;
      if(q < cols) 
	this.axpy(cols-q,-Z[k],1,q,k,Z,1,q);
    }
    float S = 1/FUtil.asum(cols,Z,1);
    FUtil.scal(cols,S,Z,1);
    ynorm *= S;

    return ynorm;
  }

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

    if(job != 0) {
      int nm = cols-1;
      B[nm] = B[nm]/this.Mat[nm][nm];
      if(cols > 0) {
	for(int jj=1;jj<cols;jj++) {
	  int j = cols - jj - 1;
	  B[j] -= this.dot(jj,1,j+1,j,B,1,j+1);
	  B[j] = B[j]/this.Mat[j][j];
	}
      }
    }
    else {
      B[0] = B[0]/this.Mat[0][0];
      for(int j=1;j<cols;j++) {
	float T = -B[j-1];
	this.axpy(cols-j, T, 1,j,j-1,B,1,j);
	B[j] = B[j]/this.Mat[j][j];    
      }
    }
  }
  
  public void inverse() {
    for(int k=cols-1;k>=0;k--) {
      if(this.Mat[k][k] == 0) 
	return;
      else {
	this.Mat[k][k] = 1/this.Mat[k][k];
	float T = -this.Mat[k][k];
	this.scal(cols-k-1,T,1,k+1,k);
	for(int j=0;j<k;j++) {
	  T = this.Mat[k][j];
	  this.Mat[k][j] = 0;
	  this.axpy(cols-k,T,1,k,k,this,1,k,j);
	}
      }
    }
  }
}