/*
 ##  class DTriUp
 ##  
 ##  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.DMatrix;
import edu.rice.linpack.util.*;
import edu.rice.linpack.Vector.*;
import edu.rice.linpack.LNumber.*;

public class DTriUp extends DTriFull {

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

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

  /*  All of the methods are tested and work well except perhaps
      condition.  
	*/

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

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

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

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

    ynorm = this.solveUZ(Z,1);
    }
    LDouble R;
    
    if(tnorm != 0) 
      R = new LDouble(ynorm/tnorm);
    else
      R = new LDouble(tnorm);
    return R;
  }

  public void solve(Vector Be, int job) 
       throws WrongDataTypeException
  {

    double[] B = Be.getDoubleArray();

    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--) {
	  double 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];
	double 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);
	}
      }
    }
  }
}