/*
 ##  class DBanded
 ##  
 ##  This class holds the frame for matrices stored in banded form.
 ##  It also contains the Gaussian routines for banded matrices.
 ##
 */


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

public class DBanded extends DMatrix {

  protected int mu, ml;

  public DBanded() {
    super();
    cons();
  }
  public DBanded(int ML, int MU, int columns) {
    cons(2*ML+MU+1,columns);
    mu = MU;
    ml = ML;
  }
  // This is for unbanded double arrays //
  public DBanded(double[][] f, int ML, int MU) {
    mu = MU;
    ml = ML;
    Mat = band(f,ml,mu);
    rows = Mat.length;
    cols = Mat[0].length;
    pivot = new int[cols];
  }
  // This is for banded double arrays //
  public DBanded(double[][] f, int ML) {
    cons(f);
    ml = ML;
    mu = f.length - 1 - 2*ML;
  }
  public DBanded(DFull F, int Ml, int Mu) {
    mu = Mu;
    ml = Ml;
    Mat = band(F.Mat,ml,mu);
    rows = Mat.length;
    cols = Mat[0].length;
    pivot = new int[cols];
  }
  public DBanded(DBanded F) {
    cons(F);
    mu = F.mu;
    ml = F.ml;
  }

  private double[][] band(double[][] f, int Ml, int Mu) {
    int n = f[0].length;
    double[][] Dband = new double[2*ml+mu+1][n];
    int M = Ml + Mu;
    for(int j=0;j<n;j++) {
      int i1 = Math.max(0,j-Mu);
      int i2 = Math.min(n,j+Ml+1);
      for(int i=i1;i<i2;i++) {
	int k = i - j + M;
	Dband[k][j] = f[i][j];
      }
    }
    return Dband;
  }
  public DFull unband() {
    DFull Fband = new DFull(cols,cols);
    int M = ml + mu;
    for(int j=0;j<cols;j++) {
      int i1 = Math.max(0,j-mu);
      int i2 = Math.min(cols,j+ml+1);
      for(int i=i1;i<i2;i++) {
	int k = i - j + M;
	Fband.Mat[i][j] = this.Mat[k][j];
      }
    }
    return Fband;
  }  
  
  public double oneNorm() {
    double anorm = 0;
    int L = ml;
    int is = L + mu;
    for(int j=0;j<cols;j++) {
      anorm = Math.max(anorm,this.asum(L+1,1,is,j));
      if(is > ml) 
	is--;
      if(j < mu)
	L++;
      if(j >= cols-ml-1)
	L--;  
    }
    return anorm;
  }

  /*  Pivots are java based  */

  public void factor() 
       throws SingularMatrixException
  {
    
    int info = 0;
    int M = ml + mu + 1;
    double T;
    
    int Ji = Math.min(cols,M) - 1;
    int io;
    for(int Jz=mu+1;Jz < Ji;Jz++) {
      io = M - 1 - Jz;
      for(int i=io;i < ml;i++) 
	this.Mat[i][Jz] = 0;
    }
    
    int Jz = Ji;
    int Ju = 0;
    
    /* Gaussian elimination with partial pivoting */

    int colsm = cols - 1;
    if(cols > 1) {
      int kp;
      for(int k=0;k<colsm;k++,Jz++) {
	kp = k+1;
	if(Jz < cols) 
	  for(int i=0;i<ml;i++) 
	    this.Mat[i][Jz] = 0;
	
	int lm = Math.min(ml, colsm-k) - 1;
	int L = this.i_amax(lm+2,1,M-1,k) + M - 1;
	pivot[k] = L + kp - M;

	if(this.Mat[L][k] != 0) {
	  if(L != M-1) 
	    this.swapElem(L,k,M-1,k);
	  
	  T = -1.0/this.Mat[M-1][k];
	  this.scal(lm+1,T,1,M,k);
	  
	  /*  ROW ELIMINATION WITH COLUMN INDEXING */

	  Ju = Math.min(Math.max(Ju,mu+pivot[k]),colsm);
	  int MM = M-2;
	  for(int j=kp;j<=Ju;j++) {
	    L--;
	    T = this.Mat[L][j];
	    if(L != MM) 
	      this.swapElem(L,j,MM,j);
	    this.axpy(lm+1, T,1,M,k,this,1,MM+1,j);
	    MM--;
	  }
	}
	else 
	  info = k+1;
      }
    }
    pivot[cols-1] = cols-1;
    if(this.Mat[M-1][cols-1] == 0) 
      info = cols;
    if(info != 0)
      throw new SingularMatrixException(info);
  }

  public LNumber condition() 
       throws SingularMatrixException, WrongDataTypeException
  {
    Vector Z = new DVector(cols);
    return this.condition(Z);
  }
  public LNumber condition(Vector Ze) 
       throws SingularMatrixException, WrongDataTypeException
  {
    double[] Z = Ze.getDoubleArray();
    
    // Compute 1-Norm of this //

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

    // Factor //

    try {
      this.factor();
    } finally {
    
      // Rcond = 1/(norm(this) * (Estimate of norm(inv(this)))) //
      
      for(int j=0;j<cols;j++)
	Z[j] = 0;
      
      int M = ml + mu;
      
      this.solveTransUW(M,Z);

      double S = 1/DUtil.asum(cols,Z,1);
      DUtil.scal(cols,S,Z,1);
      
      this.solveLV(M,Z);
      
      ynorm = this.solveTransLY(M,Z);

      // Solve U*Z = W //
      
      ynorm = this.solveUZ(M,Z,ynorm);
    }   
    LDouble R;

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

    return R;
  }
  private void solveTransUW(int M, double[] Z) {
    double ek = 1;
    int Ju = 0;
    for(int k=0;k<cols;k++) {
      if(Z[k] != 0) 
	ek = DUtil.signOf(ek,-Z[k]);
      if(Math.abs(ek-Z[k]) > Math.abs(this.Mat[M][k])) {
	double S = Math.abs(this.Mat[M][k])/Math.abs(ek-Z[k]);
	DUtil.scal(cols,S,Z,1);
	ek *= S;
      }
      double wk = ek - Z[k];
      double wkm = -ek - Z[k];
      double S = Math.abs(wk);
      double SM = Math.abs(wkm);
      if(this.Mat[M][k] != 0) {
	wk = wk/this.Mat[M][k];
	wkm = wkm/this.Mat[M][k];
      }
      else {
	wk = 1;
	wkm = 1;
      }
      int kp = k + 1;
      Ju = Math.min(Math.max(Ju,mu+pivot[k]+1),cols);
      int MM = M-1;
      if(kp < Ju) {
	for(int j=kp;j<Ju; j++,MM--) {
	  SM += Math.abs(Z[j]+wkm*this.Mat[MM][j]);
	  Z[j] += wk*this.Mat[MM][j];
	  S += Math.abs(Z[j]);
	}
	if(S < SM) {
	  double T = wkm - wk;
	  wk = wkm;
	  MM = M-1;
	  for(int j=kp;j<Ju;j++,MM--) 
	    Z[j] += T*this.Mat[MM][j];
	}
      }
      Z[k] = wk;
    }
  }
  private void solveLV(int M, double[] Z) {
    
    for(int k=cols-1;k>=0;k--) {
      int lm = Math.min(ml,cols-k-1);
      if(k < cols-1)
	Z[k] += this.dot(lm, 1,M+1,k, Z, 1, k+1);
      if(Math.abs(Z[k]) > 1) {
	double S = 1.0/Math.abs(Z[k]);
	DUtil.scal(cols,S,Z,1);
      }
      int L = pivot[k];
      DUtil.swapElems(Z,L,k);
    }
    
    double S = 1.0/DUtil.asum(cols,Z,1);
    DUtil.scal(cols,S,Z,1);
  }
  private double solveTransLY(int M, double[] Z) {
    double ynorm = 1;

    for(int k=0;k<cols;k++) {
      int L = pivot[k];
      DUtil.swapElems(Z,L,k);
      int lm = Math.min(ml,cols-k-1);
      if(k < cols)
	this.axpy(lm,Z[k],1,M+1,k, Z, 1, k+1);
      if(Math.abs(Z[k]) > 1) {
	double S = 1.0/Math.abs(Z[k]);
	DUtil.scal(cols,S,Z,1);
	ynorm *= S;
      }
    }
    double S = 1.0/DUtil.asum(cols,Z,1);
    DUtil.scal(cols,S,Z,1);
    ynorm *= S;
    
    return ynorm;
  }
  double solveUZ(int M, double[] Z, double ynorm){
    
    for(int k=cols-1;k>=0;k--) {
      if(Math.abs(Z[k]) > Math.abs(this.Mat[M][k])) {
	double S = Math.abs(this.Mat[M][k])/Math.abs(Z[k]);
	DUtil.scal(cols,S,Z,1);
	ynorm *= S;
      }
      
      if(this.Mat[M][k] != 0)
	Z[k] /= this.Mat[M][k];
      else
	Z[k] = 1;
      int lm = Math.min(k,M);
      int la = M - lm;
      int lz = k - lm;
      double T = -Z[k];
      this.axpy(lm,T,1,la,k,Z,1,lz);
    }

    double S = 1.0/DUtil.asum(cols,Z,1);
    DUtil.scal(cols,S,Z,1);
    ynorm *= S;
  
    return ynorm;
  }
  
  public void solve(Vector B) 
       throws WrongDataTypeException
  {
    this.solve(B,0);
  }
  public void solve(Vector Be, int Job) 
       throws WrongDataTypeException
  {
    
    double[] B = Be.getDoubleArray();

    int M = mu+ml+1;
    int nm = cols-1;
    double T;

    if(Job == 0) {
      if(ml != 0) {
	if(nm > 0) {
	  for(int k=0;k<nm;k++) {
	    int lm = Math.min(ml,nm-k);
	    int l = pivot[k];
	    T = B[l];
	    if(l != k) 
	      DUtil.swapElems(B,l,k);
	    this.axpy(lm,T,1,M,k,B,1,k+1);
	  }
	}
      }

      //  Solve U*X=Y  //

      this.solveRX(M-1,B);

    }
    else {

      //  Job is non-zero, solve trans(A)*X=B  //

      this.solveTransRY(M-1,B);

      if(ml != 0) {
	for(int k=cols-2;k>=0;k--) {
	  int lm = Math.min(ml,nm-k);
	  B[k] += this.dot(lm,1,M,k,B,1,k+1);
	  int l = pivot[k];
	  if(l != k) 
	    DUtil.swapElems(B,l,k);
	}
      }
    }
    return; 
  }
  void solveTransRY(int M, double[] B) {
    for(int k=0;k<cols;k++) {
      int lm = Math.min(k,M);
      int la = M-lm;
      int lb = k-lm;
      double T = this.dot(lm,1,la,k,B,1,lb);
      B[k] = (B[k] - T)/this.Mat[M][k];
    }
  }
  void solveRX(int M, double[] B){
    for(int k=cols-1;k>=0;k--) {
      int LM = Math.min(k,M);
      int la = M - LM;
      int lb = k - LM;
      B[k] /= this.Mat[M][k];
      double T = -B[k];
      this.axpy(LM,T,1,la,k,B,1,lb);
    }
  }
	    
  public Vector determ() {
    
    double[] Det = new double[2];
    int M = mu + ml + 1;
    
    Det[0] = 1;
    Det[1] = 0;

    for(int i=0;i<cols;i++) {
      if(pivot[i] != i) 
	Det[0] = -Det[0];
      Det[0] *= this.Mat[M-1][i];
      
      if(Det[0] == 0) {
	DVector V = new DVector(Det);
	return V;
      }      
      else 
	DUtil.detNorm(Det);
    }
    DVector V2 = new DVector(Det);
    return V2;
  }
}