/*
 ##  class FTriUp
 ##  
 ##  This class holds the routines for a upper triangular matrix stored
 ##  in full form.  There is no factor routine; triangular matrices are
 ##  already in factored form.
 ##
 */

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

public class FTriUp extends FTriFull {

  public FTriUp() {
    super();
  }
  public FTriUp(int i, int j) {
    Mat = new float[i][j];
    rows = i;
    cols = j;
  }
  public FTriUp(float[][] F) {
    rows = F.length;
    cols = (F[0]).length;
    Mat = F;
  }
  public FTriUp(FTriUp 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++) 
      tnorm = Math.max(tnorm,this.asum(j+1,1,0,j));
    return tnorm;
  }

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

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

    try {
      this.zeroDiag();
    } finally {
    
    for(int j=0;j<cols;j++) 
      Z[j] = 0;
        
    this.solveTransUW(Z);

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

    ynorm = this.solveUZ(Z,1);
    }
    LFloat R = new LFloat();
    
    if(tnorm != 0) 
      R.setFloat(ynorm/tnorm);

    return R;
  }

  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 j=nm-1;j>=0;j--) {
	  float T = -B[j+1];
	  this.axpy(j+1, T, 1,0,j+1,B,1,0);
	  B[j] = B[j]/this.Mat[j][j];    
	}
      }
    }
    else {
      B[0] = B[0]/this.Mat[0][0];
      if(cols > 1) {
	for(int j=1;j<cols;j++) {
	  B[j] -= this.dot(j,1,0,j,B,1,0);
	  B[j] = B[j]/this.Mat[j][j];
	}
      }
    }
  }
  
  public void inverse() {
    for(int k=0;k<cols;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(k,T,1,0,k);
	int kp = k+1;
	for(int j=kp;j<cols;j++) {
	  T = this.Mat[k][j];
	  this.Mat[k][j] = 0;
	  this.axpy(k+1,T,1,0,k,this,1,0,j);
	}
      }
    }
  }
}