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

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

public class DPoPack extends DPacked {

  public DPoPack() {
    super();
    order = 0;
  }  
  public DPoPack(int o) {
    order = o;
    double r = .5*(o*(o+1));
    rows = (int) r;
    cols = 1;
    Mat = new double[rows][1];
  }
  public DPoPack(double[] f, int o) {
    rows = f.length;
    cols = 1;
    Mat = new double[rows][cols];
    for(int i=0;i<rows;i++) 
      Mat[i][0] = f[i];
    order = o;
  }
  public DPoPack(double[][] f, int o) {
    rows = f.length;
    cols = 1;
    order = o;
    Mat = new double[rows][1];
    for(int i=0;i<rows;i++) 
      Mat[i][0] = f[i][0];
  }
  public DPoPack(double[][] f) {
    cols = 1;
    order = f.length;
    double r = .5*(order*(order + 1));
    rows = (int) r;
    Mat = new double[rows][1];
    Mat = pack(f);
  }
  public DPoPack(DPoPack F) {
    rows = F.rows;
    cols = 1;
    Mat = new double[rows][1];
    for(int i=0;i<rows;i++)
      for(int j=0;j<cols;j++) 
	Mat[i][j] = F.Mat[i][j];
    order = F.order;
  }
  public DPoPack(DPoFull F) {
    cols = 1;
    order = F.numofRows();
    double r = .5*(order*(order + 1));
    rows = (int)r;
    Mat = new double[rows][1];
    Mat = pack(F.Mat);
  }

  public DPoFull unpack() {
    int k = 0;
    int n = this.order;
    DPoFull R = new DPoFull(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;
  }

  /* All have been completely tested
     
   */

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

    for(int j=0;j<order;j++) {
      S = 0;
      if(j > 0) {
	int kj = (int)(.5*(j*(j+1)));
	int kk = -1;
	for(int k=0;k<j;k++) {
	  double T = this.Mat[kj][0] - this.dot(k,1,kk+1,0,this,1,jj+1,0);
	  kk += k+1;
	  T = T/this.Mat[kk][0];
	  this.Mat[kj][0] = T;
	  S += T*T;
	  kj++;
	}
      }
      jj += j+1;
      S = this.Mat[jj][0] - S;
      if(S <= 0) 
	throw new SingularMatrixException(j+1);
      else 
	this.Mat[jj][0] = Math.sqrt(S);
    }
  } 

  public LNumber condition() 
       throws SingularMatrixException, WrongDataTypeException
  {
    Vector Z = new DVector(order);
    return this.condition(Z);
  }
  public LNumber condition(Vector Ze) 
       throws SingularMatrixException, WrongDataTypeException
  {
    double[] Z = Ze.getDoubleArray();

    double anorm = this.oneNorm();

    //  Factor  //

    this.factor();

    for(int j=0;j<order;j++) 
      Z[j] = 0;
    
    this.solveTransUW(Z);
      
    double S = 1/DUtil.asum(order,Z,1);
    DUtil.scal(order,S,Z,1);
     
    this.solveUZ(Z);

    double ynorm = this.solveTransRV(Z);
    
    ynorm = this.solveUZ(Z,ynorm);
     
    LDouble R;

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

    return R;
  }
  private void solveTransUW(double[] Z) {
    double ek = 1;
    int kk = -1;
    for(int k=0;k<order;k++) {
      int kp = k+1;
      kk += kp;
      if(Z[k] != 0) 
	ek = DUtil.signOf(ek,-Z[k]);
      if(Math.abs(ek-Z[k]) > this.Mat[kk][0]) {
	double S = this.Mat[kk][0]/Math.abs(ek-Z[k]);
	DUtil.scal(order,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[kk][0];
      wkm = wkm/this.Mat[kk][0];
      int kj = kk + kp;
      if(kp < order) {
	for(int j=kp;j<order;j++) {
	  SM += Math.abs(Z[j] + wkm*this.Mat[kj][0]);
	  Z[j] += wk*this.Mat[kj][0];
	  S += Math.abs(Z[j]);
	  kj += j+1;
	}
	if(S < SM) {
	  double T = wkm - wk;
	  wk = wkm;
	  kj = kk + kp;
	  for(int j=kp;j<order;j++) {
	    Z[j] += T*this.Mat[kj][0];
	    kj += j+1;
	  }
	}
      }
      Z[k] = wk;
    }
  }
  private double solveUZ(double[] Z, double ynorm) {
    
    int kk = order*(order+1)/2 - 1;
    for(int k=order-1;k>=0;k--) {
      if(Math.abs(Z[k]) > this.Mat[kk][0]) {
	double S = this.Mat[kk][0]/Math.abs(Z[k]);
	DUtil.scal(order,S,Z,1);
	ynorm *= S;
      }
      Z[k] = Z[k]/this.Mat[kk][0];
      kk = kk - k - 1;
      this.axpy(k,-Z[k],1,kk+1,0,Z,1,0);
    }
    double S = 1/DUtil.asum(order,Z,1);
    DUtil.scal(order,S,Z,1);
    ynorm *= S;

    return ynorm;
  }
  private void solveUZ(double[] Z) {
    
    int kk = order*(order+1)/2 - 1;
    for(int k=order-1;k>=0;k--) {
      if(Math.abs(Z[k]) > this.Mat[kk][0]) {
	double S = this.Mat[kk][0]/Math.abs(Z[k]);
	DUtil.scal(order,S,Z,1);
      }
      Z[k] = Z[k]/this.Mat[kk][0];
      kk = kk - k - 1;
      this.axpy(k,-Z[k],1,kk+1,0,Z,1,0);
    }
    double S = 1/DUtil.asum(order,Z,1);
    DUtil.scal(order,S,Z,1);
  }
  private double solveTransRV(double[] Z) {
    
    int kk = -1;
    double ynorm = 1;
    for(int k=0;k<order;k++) {
      Z[k] -= this.dot(k,1,kk+1,0,Z,1,0);
      kk += k+1;
      if(Math.abs(Z[k]) > this.Mat[kk][0]) {
	double S = this.Mat[kk][0]/Math.abs(Z[k]);
	DUtil.scal(order,S,Z,1);
	ynorm *= S;
      }
      Z[k] = Z[k]/this.Mat[kk][0];
    }
    double S = 1/DUtil.asum(order,Z,1);
    DUtil.scal(order,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();

    int kk = -1;
    for(int k = 0;k<order;k++) {
      double T = this.dot(k,1,kk+1,0,B,1,0);
      kk += k+1;
      B[k] = (B[k] - T)/this.Mat[kk][0];
    }
    for(int kb=0;kb<order;kb++) {
      int k = order - kb - 1;
      B[k] = B[k]/this.Mat[kk][0];
      kk -= k+1;
      this.axpy(k,-B[k],1,kk+1,0,B,1,0);
    }
  }
  
  public Vector determ() {
    
    double[] Det = new double[2];
    
    Det[0] = 1;
    Det[1] = 0;
    
    int ii = 0;
    for(int k=0;k<order;k++) {
      Det[0] *= Math.pow(this.Mat[ii][0],2); 
      
      if(Det[0] == 0) {
	DVector V = new DVector(Det);
	return V;
      }
      else 
	DUtil.detNorm(Det);
      ii += k+2;
    }
    DVector V2 = new DVector(Det);
    return V2;
  }

  public void inverse() {
    int kk = 0;
    int k1 = 0;
    for(int k=0;k<order;k++) {
      this.Mat[kk][0] = 1/this.Mat[kk][0];
      double T = -this.Mat[kk][0];
      this.scal(k,T,1,k1,0);
      int kp = k+1;
      int j1 = kk+1;
      int kj = kk+k+1;
      if(order > kp) {
	for(int j=kp;j<order;j++) {
	  T = this.Mat[kj][0];
	  Mat[kj][0] = 0;
	  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;
      if(j > 0) {
	for(int k=0;k<j;k++) {
	  double T = this.Mat[kj][0];
	  this.axpy(k+1,T,1,j1,0,this,1,k1,0);
	  k1 += k+1;
	  kj++;
	}
      }
      double T = this.Mat[jj][0];
      this.scal(j+1,T,1,j1,0);
      j1 = jj + 1;
      jj += j+2;
    }
  }
}