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

public class NBanded extends NMatrix {

  protected int mu, ml;
  
  public NBanded() {
    super();
    cons();
  }
  public NBanded(int ML, int MU, int columns) {
    cons(2*ML+MU+1,columns);
    mu = MU;
    ml = ML;
  }
  // This is for unbanded LNumber arrays //
  public NBanded(LNumber[][] 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 LNumber arrays //
  public NBanded(LNumber[][] f, int ML) {
    cons(f);
    ml = ML;
    mu = f.length - 1 - 2*ML;
  }
  public NBanded(NFull 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 NBanded(NBanded F) {
    mu = F.mu;
    ml = F.ml;
    rows = F.numofRows();
    cols = F.numofCols();
    pivot = new int[cols];
    int M = ml+mu;
    Mat = new LNumber[rows][cols];
    for(int j=0;j<cols;j++) {
      int low = Math.max(0,M-j);
      int hi;
      if(ml > 0)
	hi = Math.min(rows,rows+M-j+1);
      else 
	hi = rows;
      for(int i=low;i<hi;i++) 
	Mat[i][j] = (F.getElem(i,j)).Clone();
    }
  }

  private static LNumber[][] band(LNumber[][] f, int Ml, int Mu) {
    int n = f[0].length;
    LNumber[][] Fband = new LNumber[2*Ml+Mu+1][n];
    int M = Ml + Mu;
    for(int j=0;j<n;j++) {
      int i1 = Math.max(0,M-j);
      int i2 = Math.min(n,n+M-j);
      for(int i=0;i<n;i++) {
	int k = i - j + M;
	if(k>=i1 && k<i2 && k<(M+Ml+1)) {
	  if(k>=Ml) 
	    Fband[k][j] = f[i][j];
	  else {
	    LNumber T = f[0][0].Clone();
	    T.setZero();
	    Fband[k][j] = T;
	  }
	}
      }
    }
    return Fband;
  }
  public NFull unband() {
    NFull Fband = new NFull(cols,cols);
    int M = ml + mu;
    for(int j=0;j<cols;j++) {
      int low = Math.max(0,M-j);
      int hi = Math.min(cols,cols+M-j);
      for(int i=0;i<cols;i++) {
	int k = i - j + M;
	if(k>=low && k<hi && k<rows) 
	  Fband.setElem(i,j,Mat[k][j]);
	else {
	  LNumber T = Mat[M][M].Clone();
	  T.setZero();
	  Fband.setElem(i,j,T);
	}
      }
    }
    return Fband;
  }  
  public void Print() {
    System.out.println(); 
    for(int j=0;j<rows;j++) {
      int i1 = Math.max(0,ml+mu-j);
      int i2 = Math.min(cols,cols+mu+ml-j);
      for(int i=0;i<cols;i++) {
	if(i>=i1 && i<i2) {
	  if(j>=ml || Mat[j][i] != null) {
	    Mat[j][i].Print();
	    System.out.print("   ");
	  }
	  else
	    System.out.print("0.00   ");
	}
	else
	  System.out.print("0.00   ");
      }
      System.out.println();
    }
    System.out.println();
  }

  public LNumber oneNorm() {
    int L = ml;
    int is = L + mu;
    LNumber anorm = this.asum(L+1,1,is,0);
    for(int j=0;j<cols;j++) {
      if(is > ml) 
	is--;
      if(j < mu)
	L++;
      if(j >= cols-ml-1)
	L--; 
      if(j < cols-1)
	anorm = anorm.max(this.asum(L+1,1,is,j+1));
    }
    return anorm;
  }

  /*  Pivots are java based  */

  public void factor() 
       throws SingularMatrixException
  {
    
    int info = 0;
    int M = ml + mu + 1;
    LNumber 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++) 
	Mat[i][Jz].setZero();
    }
    
    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++) 
	    Mat[i][Jz].setZero();
	
	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(!(Mat[L][k].equals(0))) {
	  if(L != M-1) 
	    this.swapElem(L,k,M-1,k);
	  T = (Mat[M-1][k].inv()).negateTo();
	  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 = Mat[L][j].Clone();
	    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(Mat[M-1][cols-1].equals(0)) 
      info = cols;
    if(info != 0) 
      throw new SingularMatrixException(info);
  }

  public LNumber condition() 
       throws SingularMatrixException
  {
    LNumber[] Z = new LNumber[cols];
    return this.condition(Z);
  }
  public LNumber condition(LNumber[] Z) 
       throws SingularMatrixException
  {
    // Compute 1-Norm of this //

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

    // Factor //

    try {
      this.factor();
    } finally {

      // Rcond = 1/(norm(this) * (Estimate of norm(inv(this)))) //

      int M = ml + mu;
      
      for(int j=0;j<cols;j++) {
	Z[j] = Mat[M][0].Clone();
	Z[j].setZero();
      }
      
      this.solveTransUW(M,Z);

      LNumber S = (NUtil.asum(cols,Z,1)).invTo();
      NUtil.scal(cols,S,Z,1);

      this.solveLV(M,Z);
    
      ynorm = this.solveTransLY(M,Z);

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

    if(!anorm.equals(0))
      R = ynorm.div(anorm);
    else
      R = anorm;

    return R;
  }
  private void solveTransUW(int M, LNumber[] Z) {
    LNumber ek = Mat[mu+ml][0].Clone();
    ek.setOne();
    int Ju = 0;
    for(int k=0;k<cols;k++) {
      if(!Z[k].equals(0)) 
	ek = NUtil.signOfA(ek,Z[k].negate());
      if(((ek.sub(Z[k])).abs()).greaterThan(Mat[M][k].abs())) {
	LNumber S = (Mat[M][k].abs()).div((ek.sub(Z[k])).abs());
	NUtil.scal(cols,S,Z,1);
	ek.multTo(S);
      }
      
      LNumber wk = ek.sub(Z[k]);
      LNumber wkm = (ek.negate()).sub(Z[k]);
      LNumber S = wk.abs();
      LNumber SM = wkm.abs();
      if(!Mat[M][k].equals(0)) {
	wk.divTo(Mat[M][k]);
	wkm.divTo(Mat[M][k]);
      }
      else {
	wk.setOne();
	wkm.setOne();
      }
      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.addTo((Z[j].add(wkm.mult(Mat[MM][j]))).abs());
	  Z[j].addTo(wk.mult(Mat[MM][j]));
	  S.addTo(Z[j].abs());
	}
	if(SM.greaterThan(S)) {
	  LNumber T = wkm.sub(wk);
	  wk = wkm.Clone();
	  MM = M-1;
	  for(int j=kp;j<Ju;j++,MM--) 
	    Z[j].addTo(T.mult(Mat[MM][j]));
	}
      }
      Z[k] = wk.Clone();
    }
  }
  private void solveLV(int M, LNumber[] Z) {
    
    for(int k=cols-1;k>=0;k--) {
      int lm = Math.min(ml,cols-k-1);
      if(k < cols-1) 
	Z[k].addTo(dot(lm,1,M+1,k,Z,1,k+1));
      if((Z[k].abs()).greaterThan(1)) {
	LNumber S = (Z[k].abs()).invTo();
	NUtil.scal(cols,S,Z,1);
      }
      int L = pivot[k];
      NUtil.swapElems(Z,L,k);
    }
    
    LNumber S = (NUtil.asum(cols,Z,1)).invTo();
    NUtil.scal(cols,S,Z,1);
  }
  private LNumber solveTransLY(int M, LNumber[] Z) {
    LNumber ynorm = Z[0].abs();
    ynorm.setOne();
    
    for(int k=0;k<cols;k++) {
      int L = pivot[k];
      NUtil.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((Z[k].abs()).greaterThan(1)) {
	LNumber S = (Z[k].abs()).invTo();
	NUtil.scal(cols,S,Z,1);
	ynorm.multTo(S);
      }
    }
    LNumber S = (NUtil.asum(cols,Z,1)).invTo();
    NUtil.scal(cols,S,Z,1);
    ynorm.multTo(S);
    
    return ynorm;
  }
  LNumber solveUZ(int M, LNumber[] Z, LNumber ynorm){
    
    for(int k=cols-1;k>=0;k--) {
      if((Z[k].abs()).greaterThan(Mat[M][k].abs())) {
	LNumber S = (Mat[M][k].abs()).div(Z[k].abs());
	NUtil.scal(cols,S,Z,1);
	ynorm.multTo(S);
      }
      
      if(!Mat[M][k].equals(0)) 
	Z[k].divTo(Mat[M][k]);
      else 
	Z[k].setOne();
      int lm = Math.min(k,M);
      int la = M - lm;
      int lz = k - lm;
      LNumber T = Z[k].negate();
      this.axpy(lm,T,1,la,k,Z,1,lz);
    }

    LNumber S = (NUtil.asum(cols,Z,1)).invTo();
    NUtil.scal(cols,S,Z,1);
    ynorm.multTo(S);
  
    return ynorm;
  }
  
  public void solve(LNumber[] B) {
    this.solve(B,0);
  }
  public void solve(LNumber[] B, int Job) {
    
    int M = mu+ml+1;
    int nm = cols-1;
    LNumber 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].Clone();
	    if(l != k) 
	      NUtil.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] = B[k].add(dot(lm,1,M,k,B,1,k+1));
	  int l = pivot[k];
	  if(l != k) 
	    NUtil.swapElems(B,l,k);
	}
      }
    }
    return; 
  }
  void solveTransRY(int M, LNumber[] B) {
    for(int k=0;k<cols;k++) {
      int lm = Math.min(k,M);
      int la = M-lm;
      int lb = k-lm;
      LNumber T = this.dot(lm,1,la,k,B,1,lb);
      B[k] = (B[k].sub(T)).div(Mat[M][k]);
    }
  }
  void solveRX(int M, LNumber[] 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].divTo(Mat[M][k]);
      this.axpy(LM,B[k].negate(),1,la,k,B,1,lb);
    }
  }
	    
  public LNumber[] determ() {
    
    LNumber[] Det = new LNumber[2];
    int M = mu + ml;
    
    Det[0] = Mat[M][0].Clone();
    Det[1] = Det[0].Clone();
    Det[0].setOne();
    Det[1].setZero();

    for(int i=0;i<cols;i++) {
      if(pivot[i] != i) 
	Det[0].multTo(-1);
      Det[0].multTo(Mat[M][i]);
      
      if(Det[0].equals(0)) 
	return Det;
      else 
	NUtil.detNorm(Det);
    }
    return Det;
  }
}