/* 
 ##  class NPoBand
 ##
 ##  This class contains the routines for positive definite matrices
 ##  in banded form.
 ##
 */

package edu.rice.linpack.Matrix.NMatrix;
import edu.rice.linpack.LNumber.*;
import edu.rice.linpack.util.*;

public class NPoBand extends NBanded {
  
  private int M;

  public NPoBand() {
    super();
  }
  public  NPoBand(int c, int M2) {
    super(0,M2,c);
    mu = M2;
    ml = 0;
    M = M2;
  }
  public NPoBand(LNumber[][] f, int ML, int MU) {
    mu = MU;
    M = mu;
    ml = 0;
    rows = M+1;
    cols = f[0].length;
    Mat = band(f,M);
  }
  public NPoBand(LNumber[][] f, int M2) {
    mu = M2;
    M = M2;
    ml = 0;
    rows = M+1;
    cols = f[0].length;
    Mat = f;
  }
  public NPoBand(NPoFull F, int M2) {
    mu = M2;
    M = M2;
    ml = 0;
    cols = F.numofCols();
    rows = M+1;
    Mat = band(F.Mat,M);
  }
  public NPoBand(NPoBand F) {
    super(F);
    M = F.M;
  }

  private LNumber[][] band(LNumber[][] f,int M) {
    int n = f[0].length;
    LNumber[][] ret = new LNumber[M+1][n];
    for(int j=0;j<n;j++) {
      int i1 = Math.max(0,j-M);
      for(int i=i1;i<=j;i++) {
	int k = i-j+M;
	ret[k][j] = f[i][j];
      }
    }
    return ret;
  }
  public NPoFull poUnband() {
    NPoFull Ret = new NPoFull(cols,cols);
    for(int j=0;j<cols;j++) {
      int i1 = Math.max(0,j-mu);
      for(int i=i1;i<=j;i++) {
	int k = i-j+mu;
	Ret.Mat[i][j] = this.Mat[k][j];
      }
    }
    return Ret;
  }

  public LNumber oneNorm() {
    
    LNumber[] Z = new LNumber[cols];

    for(int j=0;j<cols;j++) {
      int l = Math.min(j,M);
      int mu = Math.max(M-j,0);
      Z[j] = this.asum(l+1, 1,mu,j);
      int k = j - l;
      for(int i=mu;i<M;i++,k++) 
	Z[k].addTo(Mat[i][j].abs());
    }
    LNumber anorm = Z[0];
    for(int j=1;j<cols;j++) 
      anorm.maxTo(Z[j]);

    return anorm.Clone();
  }

  public void factor() 
       throws SingularMatrixException
  {
    LNumber S = Mat[M][0].Clone();

    for(int j=0;j<cols;j++) {
      S.setZero();
      int ik = M;
      int jk = Math.max(j-M,0);
      int mu = Math.max(M-j,0);
      
      for(int k=mu;k<M;k++,ik--,jk++) {
	LNumber T = Mat[k][j].sub(this.dot(k-mu,1,ik,jk,this,1,mu,j));
	T.divTo(Mat[M][jk]);
	Mat[k][j] = T;
	S.addTo(T.square());
      }
      S = Mat[M][j].sub(S);    
      if(!S.greaterThan(0)) 
	throw new SingularMatrixException(j+1);
      else 
	Mat[M][j] = S.sqrt();
    }
  }
  
  public LNumber condition() 
       throws SingularMatrixException
  {
    LNumber[] Z = new LNumber[cols];
    return this.condition(Z);
  }
  public LNumber condition(LNumber[] Z) 
       throws SingularMatrixException
  {

    LNumber anorm = this.oneNorm();
   
    this.factor();
    
    LNumber ek = Mat[M][0].Clone();
    ek.setOne();
    for(int j=0;j<cols;j++) {
      Z[j] = ek.Clone();
      Z[j].setZero();
    }
    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();
      wk.divTo(Mat[M][k]);
      wkm.divTo(Mat[M][k]);
      int kp = k+1;
      int J2 = Math.min(kp+M,cols);
      int i = M-1;
      
      if(kp < J2) {
	for(int j=kp;j<J2;j++,i--) {
	  SM.addTo((Z[j].add(wkm.mult(Mat[i][j]))).abs());
	  Z[j].addTo(wk.mult(Mat[i][j]));
	  S.addTo(Z[j].abs());
	}
	if(SM.greaterThan(S)) {
	  LNumber T = wkm.sub(wk);
	  wk = wkm;
	  i = M - 1;
	  for(int j=kp;j<J2;j++,i--)
	    Z[j].addTo(T.mult(Mat[i][j]));
	}
      }
      Z[k] = wk;
    }
    LNumber S = (NUtil.asum(cols,Z,1)).invTo();
    NUtil.scal(cols,S,Z,1);
    
    this.solveUZ(M,Z);
    
    LNumber ynorm = this.solveTrans(M,Z);
      
    ynorm = this.solveUZ(M,Z,ynorm);
  
    if(anorm.equals(0)) 
      return anorm;

    ynorm.divTo(anorm);
    return ynorm;
  }
  void solveUZ(int M, LNumber[] Z) {
    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);
      }
      if(!Mat[M][k].equals(0)) 
	Z[k] = Z[k].div(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);
  }
  private LNumber solveTrans (int M, LNumber[] Z) {
    
    LNumber ynorm = Z[0].Clone();
    ynorm.setOne();

    for(int k=0;k<cols;k++) {
      int LM = Math.min(k,M);
      int la = M - LM;
      int lb = k - LM;
      Z[k].subTo(this.dot(LM,1,la,k,Z,1,lb));
      if((Z[k].abs()).greaterThan(Mat[M][k])) {
	LNumber S = Mat[M][k].div(Z[k].abs());
	NUtil.scal(cols,S,Z,1);
	ynorm.multTo(S);
      }
      Z[k].divTo(Mat[M][k]);
    }
    LNumber S = (NUtil.asum(cols,Z,1)).invTo();
    NUtil.scal(cols,S,Z,1);
    ynorm.multTo(S);
    
    return ynorm;
  }

  public void solve(LNumber[] B, int J) {
    this.solve(B);
  }
  public void solve(LNumber[] B) {
    this.solveTransRY(M,B);
    this.solveRX(M,B);
  }

  public LNumber[] determ() {

    LNumber[] Det = new LNumber[2];
    
    Det[0] = Mat[M][0].Clone();
    Det[1] = Det[0].Clone();
    Det[0].setOne();
    Det[1].setZero();
    
    for(int i=0;i<cols;i++) {
      Det[0].multTo(Mat[M][i].square());
      
      if(Det[0].equals(0)) 
	return Det;
      else 
	NUtil.detNorm(Det);
    }
    return Det;
  }
}