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

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

public class DPoBand extends DBanded {
  
  private int M;

  public DPoBand() {
    super();
  }
  public  DPoBand(int c, int M2) {
    super(0,M2,c);
    mu = M2;
    ml = 0;
    M = M2;
  }
  public DPoBand(double[][] f, int ML, int MU) {
    mu = MU;
    M = mu;
    ml = 0;
    rows = M+1;
    cols = f[0].length;
    Mat = band(f);
  }
  public DPoBand(double[][] f, int M2) {
    mu = M2;
    M = M2;
    ml = 0;
    rows = M+1;
    cols = f[0].length;
    Mat = new double[rows][cols];
    for(int i=0;i<rows;i++) 
      for(int j=0;j<cols;j++) 
	Mat[i][j] = f[i][j];
  }
  public DPoBand(DPoFull F, int M2) {
    mu = M2;
    M = M2;
    ml = 0;
    cols = F.numofCols();
    rows = M+1;
    Mat = band(F.Mat);
  }
  public DPoBand(DPoBand F) {
    super(F);
    M = F.M;
  }

  private double[][] band(double[][] f) {
    int n = f[0].length;
    double[][] ret = new double[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 DPoFull poUnband() {
    DPoFull Ret = new DPoFull(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 double oneNorm() {
    
    double[] Z = new double[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] += Math.abs(this.Mat[i][j]);
    }
    double anorm = 0;
    for(int j=0;j<cols;j++) 
      anorm = Math.max(anorm, Z[j]);

    return anorm;
  }

  public void factor() 
       throws SingularMatrixException
  {
    double S;

    for(int j=0;j<cols;j++) {
      S = 0;
      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++) {
	double T = this.Mat[k][j] - this.dot(k-mu,1,ik,jk,this,1,mu,j);
	T = T/this.Mat[M][jk];
	this.Mat[k][j] = T;
	S += T*T;
      }
      S = this.Mat[M][j] - S;    
      if(S <= 0) 
	throw new SingularMatrixException(j+1);
      else 
	this.Mat[M][j] = Math.sqrt(S);
    }
  }

  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();
    double anorm = this.oneNorm();
   
    this.factor();
    
    double ek = 1;
    for(int j=0;j<cols;j++) 
      Z[j] = 0;
      
    for(int k=0;k<cols;k++) {
      if(Z[k] != 0) 
	ek = DUtil.signOf(ek,-Z[k]);
      if(Math.abs(ek - Z[k]) > this.Mat[M][k]) {
	double S = 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);
      wk = wk/this.Mat[M][k];
      wkm = wkm/this.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 += Math.abs(Z[j] + wkm*this.Mat[i][j]);
	  Z[j] += wk*this.Mat[i][j];
	  S += Math.abs(Z[j]);
	}
	if(S < SM) {
	  double T = wkm - wk;
	  wk = wkm;
	  i = M - 1;
	  for(int j=kp;j<J2;j++,i--) 
	    Z[j] += T*this.Mat[i][j];
	}
      }
      Z[k] = wk;
    }
    double S = 1/DUtil.asum(cols,Z,1);
    DUtil.scal(cols,S,Z,1);
    
    this.solveUZ(M,Z);
    
    double ynorm = this.solveTrans(M,Z);
      
    ynorm = this.solveUZ(M,Z,ynorm);
  
    LDouble R;

    if(anorm != 0) 
      R = new LDouble(ynorm/anorm);
    else
      R = new LDouble(anorm);
    return R;
  }
  void solveUZ(int M, double[] Z) {
    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);
      }
      
      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);
  }
  private double solveTrans (int M, double[] Z) {
    
    double ynorm = 1;

    for(int k=0;k<cols;k++) {
      int LM = Math.min(k,M);
      int la = M - LM;
      int lb = k - LM;
      Z[k] -= this.dot(LM,1,la,k,Z,1,lb);
      if(Math.abs(Z[k]) > this.Mat[M][k]) {
	double S = this.Mat[M][k]/Math.abs(Z[k]);
	DUtil.scal(cols,S,Z,1);
	ynorm *= S;
      }
      Z[k] /= this.Mat[M][k];
    }
    double S = 1/DUtil.asum(cols,Z,1);
    DUtil.scal(cols,S,Z,1);
    ynorm *= S;
    
    return ynorm;
  }

  public void solve(Vector B, int J) 
       throws WrongDataTypeException
  {
    this.solve(B);
  }
  public void solve(Vector Be) 
       throws WrongDataTypeException
  {
    double[] B = Be.getDoubleArray();
    this.solveTransRY(M,B);
    this.solveRX(M,B);
  }

  public Vector determ() {

    double[] Det = new double[2];
    
    Det[0] = 1;
    Det[1] = 0;
    
    for(int i=0;i<cols;i++) {
      Det[0] *= Math.pow(this.Mat[M][i],2);
      
      if(Det[0] == 0) {
	DVector V = new DVector(Det);
	return V;
      }
      else 
	DUtil.detNorm(Det);
    }
    DVector V2 = new DVector(Det);
    return V2;
  }
}