/*  
 ##  class NPoPack
 ##
 ##  This class contains the routines for positive definite matrices
 ##  in packed form.
 ##
 */

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

public class NPoPack extends NPacked {

  public NPoPack() {
    super();
    order = 0;
  }  
  public NPoPack(int o) {
    order = o;
    float r = (float).5*(o*(o+1));
    rows = (int) r;
    cols = 1;
    Mat = new LNumber[rows][1];
  }
  public NPoPack(LNumber[] f, int o) {
    rows = f.length;
    cols = 1;
    Mat = new LNumber[rows][cols];
    for(int i=0;i<rows;i++) 
      Mat[i][0] = f[i];
    order = o;
  }
  public NPoPack(LNumber[][] f, int o) {
    rows = f.length;
    cols = 1;
    order = o;
    Mat = new LNumber[rows][1];
    Mat = f;
  }
  public NPoPack(LNumber[][] f) {
    cols = 1;
    order = f.length;
    float r = (float).5*(order*(order + 1));
    rows = (int) r;
    Mat = new LNumber[rows][1];
    Mat = pack(f);
  }
  public NPoPack(NPoPack F) {
    rows = F.numofRows();
    cols = 1;
    Mat = new LNumber[rows][1];
    for(int i=0;i<rows;i++)
      for(int j=0;j<cols;j++) 
	Mat[i][j] = (F.getElem(i,j)).Clone();
    order = F.getOrder();
  }
  public NPoPack(NPoFull F) {
    cols = 1;
    order = F.numofRows();
    float r = (float).5*(order*(order + 1));
    rows = (int)r;
    Mat = new LNumber[rows][1];
    Mat = pack(F.Mat);
  }

  public NPoFull unpack() {
    int k = 0;
    int n = this.order;
    NPoFull R = new NPoFull(n,n);
    for(int j=0;j<n;j++) {
      for(int i=0;i<=j;i++) {
	R.Mat[i][j] = this.Mat[k][0];
	if(j != i)
	  R.Mat[j][i] = this.Mat[k][0];
	k++;
      }
    }	
    return R;
  }

  public void factor() 
       throws SingularMatrixException
  {
    int jj = -1;

    LNumber S = Mat[0][0].Clone();

    for(int j=0;j<order;j++) {
      S.setZero();
      if(j > 0) {
	int kj = (int)(.5*(j*(j+1)));
	int kk = -1;
	for(int k=0;k<j;k++) {
	  LNumber T = Mat[kj][0].sub(this.dot(k,1,kk+1,0,this,1,jj+1,0));
	  kk += k+1;
	  T.divTo(Mat[kk][0]);
	  Mat[kj][0] = T;
	  S.addTo(T.square());
	  kj++;
	}
      }
      jj += j+1;
      S = Mat[jj][0].sub(S);
      if(!S.greaterThan(0)) 
	throw new SingularMatrixException(j+1);
      else 
	Mat[jj][0] = S.sqrt();
    }
  } 

  public LNumber condition() 
       throws SingularMatrixException
  {
    LNumber[] Z = new LNumber[order];
    return this.condition(Z);
  }
  public LNumber condition(LNumber[] Z) 
       throws SingularMatrixException
  {
    LNumber anorm = this.oneNorm();

    //  Factor  //

    this.factor();

    for(int j=0;j<order;j++) {
      Z[j] = Mat[0][0].Clone();
      Z[j].setZero();
    }    
    this.solveTransUW(Z);
    
    LNumber S = (NUtil.asum(order,Z,1)).invTo();
    NUtil.scal(order,S,Z,1);
     
    this.solveUZ(Z);

    LNumber ynorm = this.solveTransRV(Z);

    ynorm = this.solveUZ(Z,ynorm);
      
    if(anorm.equals(0))
      return anorm;
    
    ynorm.divTo(anorm);
    return ynorm;
  }
  private void solveTransUW(LNumber[] Z) {
    LNumber ek = Z[0].Clone();
    ek.setOne();
    int kk = -1;
    for(int k=0;k<order;k++) {
      int kp = k+1;
      kk += kp;
      if(!Z[k].equals(0)) 
	ek = NUtil.signOfA(ek,Z[k].negate());
      if(((ek.sub(Z[k])).abs()).greaterThan(Mat[kk][0])) {
	LNumber S = Mat[kk][0].div((ek.sub(Z[k])).abs());
	NUtil.scal(order,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[kk][0]);
      wkm.divTo(Mat[kk][0]);
      int kj = kk + kp;
      if(kp < order) {
	for(int j=kp;j<order;j++) {
	  SM.addTo((Z[j].add(wkm.mult(Mat[kj][0]))).abs());
	  Z[j].addTo(wk.mult(Mat[kj][0]));
	  S.addTo(Z[j].abs());
	  kj += j+1;
	}
	if(SM.greaterThan(S)) {
	  LNumber T = wkm.sub(wk);
	  wk = wkm.Clone();
	  kj = kk + kp;
	  for(int j=kp;j<order;j++) {
	    Z[j].addTo(T.mult(Mat[kj][0]));
	    kj += j+1;
	  }
	}
      }
      Z[k] = wk.Clone();
    }
  }
  private LNumber solveUZ(LNumber[] Z, LNumber ynorm) {
    
    int kk = order*(order+1)/2 - 1;
    for(int k=order-1;k>=0;k--) {
      if((Z[k].abs()).greaterThan(Mat[kk][0])) {
	LNumber S = Mat[kk][0].div(Z[k].abs());
	NUtil.scal(order,S,Z,1);
	ynorm.multTo(S);
      }
      Z[k].divTo(Mat[kk][0]);
      kk = kk - k - 1;
      this.axpy(k,Z[k].negate(),1,kk+1,0,Z,1,0);
    }
    LNumber S = (NUtil.asum(order,Z,1)).invTo();
    NUtil.scal(order,S,Z,1);
    ynorm.multTo(S);

    return ynorm;
  }
  private void solveUZ(LNumber[] Z) {
    int kk = order*(order+1)/2 - 1;
    for(int k=order-1;k>=0;k--) {
      if((Z[k].abs()).greaterThan(Mat[kk][0])) {
	LNumber S = Mat[kk][0].div(Z[k].abs());
	NUtil.scal(order,S,Z,1);
      }
      Z[k].divTo(Mat[kk][0]);
      kk = kk - k - 1;
      this.axpy(k,Z[k].negate(),1,kk+1,0,Z,1,0);
    }
    LNumber S = (NUtil.asum(order,Z,1)).invTo();
    NUtil.scal(order,S,Z,1);
  }
  private LNumber solveTransRV(LNumber[] Z) {
    
    int kk = -1;
    LNumber ynorm = Z[0].Clone();
    ynorm.setOne();
    for(int k=0;k<order;k++) {
      Z[k].subTo(this.dot(k,1,kk+1,0,Z,1,0));
      kk += k+1;
      if((Z[k].abs()).greaterThan(Mat[kk][0])) {
	LNumber S = Mat[kk][0].div(Z[k].abs());
	NUtil.scal(order,S,Z,1);
	ynorm.multTo(S);
      }
      Z[k].divTo(Mat[kk][0]);
    }
    LNumber S = (NUtil.asum(order,Z,1)).invTo();
    NUtil.scal(order,S,Z,1);
    ynorm.multTo(S);
    
    return ynorm;
  }
  
  
  public void solve(LNumber[] B, int J) {
    this.solve(B);
  }
  public void solve(LNumber[] B) {
    int kk = -1;
    for(int k = 0;k<order;k++) {
      LNumber T = this.dot(k,1,kk+1,0,B,1,0);
      kk += k+1;
      B[k] = (B[k].sub(T)).div(Mat[kk][0]);
    }
    for(int kb=0;kb<order;kb++) {
      int k = order - kb - 1;
      B[k].divTo(Mat[kk][0]);
      kk -= k+1;
      this.axpy(k,B[k].negate(),1,kk+1,0,B,1,0);
    }
  }
  
  public LNumber[] determ() {
    
    LNumber[] Det = new LNumber[2];
    
    Det[0] = Mat[0][0].Clone();
    Det[1] = Det[0].Clone();
    Det[0].setOne();
    Det[1].setZero();
    
    int ii = 0;
    for(int k=0;k<order;k++) {
      Det[0].multTo(Mat[ii][0].square()); 
      if(Det[0].equals(0)) 
	return Det;
      else 
	NUtil.detNorm(Det);
      ii += k+2;
    }
    return Det;
  }

  public void inverse() {
    int kk = 0;
    int k1 = 0;
    for(int k=0;k<order;k++) {
      Mat[kk][0].invTo();
      LNumber T = Mat[kk][0].negate();
      this.scal(k,T,1,k1,0);
      int kp = k+1;
      int j1 = kk+1;
      int kj = kk+k+1;
      for(int j=kp;j<order;j++) {
	T = Mat[kj][0].Clone();
	Mat[kj][0].setZero();
	this.axpy(k+1,T,1,k1,0,this,1,j1,0);
	j1 += j+1;
	kj += j+1;
      }
      k1 = kk+1;
      kk += k+2;
    }
    int jj = 0;
    int j1 = 0;
    for(int j=0;j<order;j++) {
      int jm = j-1;
      k1 = 0;
      int kj = j1;
      for(int k=0;k<j;k++) {
	LNumber T = Mat[kj][0].Clone();
	this.axpy(k+1,T,1,j1,0,this,1,k1,0);
	k1 += k+1;
	kj++;
      }
      LNumber T = Mat[jj][0].Clone();
      this.scal(j+1,T,1,j1,0);
      j1 = jj + 1;
      jj += j+2;
    }
  }
}