/*  
 ##  class DPoFull
 ##  
 ##  This class holds the routines for positive definite matrices in full
 ##  form.  It also holds the Cholesky routines for decomposing.  
 ##
 */

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

public class DPoFull extends DFull {

  public DPoFull() {
    super();
  }
  public DPoFull(int n, int q) {
    super(n, q);
  }
  public DPoFull(double[][] f) {
    super(f);
  }
  public DPoFull(DPoFull F) {
    super(F);
  }

  public DPoPack pack() {
    DPoPack P = new DPoPack(this.Mat);
    return P;
  }
  
  public void setPivots(int[] p) {
    int q = Math.min(p.length,cols);
    for(int i=0;i<q;i++) 
      pivot[i] = p[i];
  }
  public void setPivotElem(int i, int piv) {
    pivot[i] = piv;
  }

  public void transpose() {
  }

  public double getElem(int i, int j) {
    if(i > j)
      return super.getElem(j,i);
    return super.getElem(i,j);
  }
  public double[] getRow(int i) {
    return getRow(i,0);
  }
  public double[] getColumn(int i) {
    return getRow(i,0);
  }
  public double[] getColumn(int i, int j) {
    return getRow(i,j);
  }
  public double[] getRow(int i, int j) {
    double[] F = new double[rows - j];
    for(int k=0;k<i && k<(rows-j);k++) 
      F[k] = Mat[j+k][i];
    for(int k=i;k<(rows-j);k++) 
      F[k] = Mat[j+k][i];
    return F;
  }
  public void setColumn(int c, double[] F) {
    for(int r=0;r<c;r++) 
      Mat[r][c] = F[r];
    for(int r=c;r<rows;r++) 
      Mat[r][c] = F[r];
  }
  public void setElem(int i, int c, double L) {
    if(i > c)
      super.setElem(c,i,L);
    else
      super.setElem(i,c,L);
  }

  protected double oneNorm() {
    double[] Z = new double[cols];
    
    for(int j=0;j<cols;j++) {
      Z[j] = this.asum(j+1,1,0,j);
      for(int i=0;i<j;i++) {
	Z[i] += Math.abs(this.Mat[i][j]);
      }
    }
    double anorm = 0;
    for(int j=0;j<cols;j++) 
      anorm = Math.max(anorm,Z[j]);
    return anorm;
  }

  /*  Downdate, decompose, update, and exchange are untested 
      Also, need to test both DUtil.rotg and this.rotg

      */

  public void chPDecompose() 
       throws SingularMatrixException
  {
    Pul P = new Pul(0,-1);
    this.chPivot(P);
    this.decom(P);
  }
  public void chDecompose()
       throws SingularMatrixException
  {
    Pul P = new Pul(0,-1);
    this.decom(P);
  }
  private void decom (Pul P) 
       throws SingularMatrixException
  {
    double[] W = new double[cols];

    for(int k=0;k<cols;k++) {
      // Reduction loop //
      
      double maxdia = this.Mat[k][k];
      int kp = k+1;
      int maxln = k;

      if((k >= P.pl) && (k < P.pu)) {
	for(int q=kp;q<=P.pu;q++) {
	  if(this.Mat[q][q] > maxdia) {
	    maxdia = this.Mat[q][q];
	    maxln = q;
	  }
	}
      }
      if(maxdia < 0) 
	throw new SingularMatrixException(k+1);
      
      else {
	if(k != maxln) {
	  this.swap(k,1,0,k,this,1,0,maxln);
	  this.Mat[maxln][maxln] = this.Mat[k][k];
	  this.Mat[k][k] = maxdia;
	  DUtil.swapElems(pivot, k, maxln);
	}
	W[k] = Math.sqrt(this.Mat[k][k]);
	this.Mat[k][k] = W[k];
	for(int j=kp;j<cols;j++) {
	  if(k != maxln) 
	    this.chSwaps(k,j,maxln);
	  this.Mat[k][j] = this.Mat[k][j]/W[k];
	  W[j] = this.Mat[k][j];
	  this.axpy(j-k,-W[j],W,1,kp,1,kp,j);
	}
      }
    }
  }
  private void chPivot(Pul P) {
    
    this.swaperSki(P);

    // change pivots //
    
    P.pu = cols-1;
    for(int k=P.pu;k>=P.pl;k--) {
      if(pivot[k] < 0) {
	pivot[k] = -pivot[k];
	if(P.pu != k) {
	  this.swap(k,1,0,k,this,1,0,P.pu);
	  this.swapElem(k,k,P.pu,P.pu);
	  for(int j=k+1;j<cols;j++) 
	    this.chSwaps(k,j,P.pu);
	  DUtil.swapElems(pivot,k,P.pu);
	}
	P.pu--;
      }
    }
  }
  void swap2(int j, Pul P) {
    if(j != P.pl) {
      this.swap(P.pl,1,0,j,this,1,0,P.pl);
      this.swapElem(j,j,P.pl,P.pl);
      for(int q=P.pl+1;q<cols;q++) 
	this.chSwaps(P.pl,q,j);
      pivot[j] = pivot[P.pl];
      pivot[P.pl] = j+1;
    }
    P.pl++;
  }
  private void chSwaps(int k, int j, int q) {
    if(j < q) 
      this.swapElem(k,j,j,q);
    else if(j > q) 
      this.swapElem(k,j,q,j);
  }



  public void factor() 
       throws SingularMatrixException
  {
    for(int j=0;j<cols;j++) {
      double S = 0;
      for(int k=0;k<j;k++) {
	double T = this.Mat[k][j] - this.dot(k,1,0,k,this,1,0,j);
	T = T/this.Mat[k][k];
	this.Mat[k][j] = T;
	S += T*T;
      }
      S = this.Mat[j][j] - S;
      if(S <= 0) 
	throw new SingularMatrixException(j+1);
      else 
	this.Mat[j][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();

    //  Factor //
    
    this.factor();
    
    /*  RCond = 1/(norm(A)*(estimate of norm(inverse(A)))) .
	Estimate = norm(Z)/norm(Y) where  A*Z = Y  and  A*Y = E .
	The components of  E  are chosen to cause maximum local 
	growth in the elements of W  where  trans(R)*W = E .
	The vectors are frequently rescaled to avoid overflow.
	*/
    //  Solve trans(R)*W = E  //
    
    for(int j=0;j<cols;j++) 
      Z[j] = 0;
      
    this.solveTransUW(Z);
    
    double S = 1/DUtil.asum(cols,Z,1);
    DUtil.scal(cols,S,Z,1);

    //  Solve R*Y = W  //
    
    this.solveUZ(Z);
    
    //  Solve trans(R)*V = Y  //
    
    double ynorm = this.solveTransRV(Z);
    
    //  Solve R*Z = V  //
    
    ynorm = this.solveUZ(Z,ynorm);
    
    LDouble R;
    
    if(anorm != 0) 
      R = new LDouble(ynorm/anorm);
    else
      R = new LDouble(anorm);
    
    return R;
  }
  void solveUZ(double[] Z) {
    for(int k=cols-1;k>=0;k--) {
      if(Math.abs(Z[k]) > Math.abs(this.Mat[k][k])) {
	double S = Math.abs(this.Mat[k][k])/Math.abs(Z[k]);
	DUtil.scal(cols,S,Z,1);
      } 
      if(this.Mat[k][k] != 0) 
	Z[k] = Z[k]/this.Mat[k][k];
      else 
	Z[k] = 1;

      this.axpy(k,-Z[k],1,0,k,Z,1,0);
    }
    double S = 1.0/DUtil.asum(cols,Z,1);
    DUtil.scal(cols,S,Z,1);
  }
  private double solveTransRV (double[] Z) {
    double ynorm = 1;
    for(int k=0;k<cols;k++) {
      Z[k] = Z[k] - this.dot(k,1,0,k,Z,1,0);
      if(Math.abs(Z[k]) > this.Mat[k][k]) {
	double S = this.Mat[k][k]/Math.abs(Z[k]);
	DUtil.scal(cols,S,Z,1);
	ynorm *= S;
      }
      Z[k] = Z[k]/this.Mat[k][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 Job) 
       throws WrongDataTypeException
  {
    this.solve(B);
  }
  public void solve(Vector Be) 
       throws WrongDataTypeException
  {
    
    double[] B = Be.getDoubleArray();

    // Solve Trans(R)*Y = B //

    this.solveTransAX(B);

    // Solve R*X = Y  //

    this.solveUX(B);

    return;
  }
  
  public Vector determ() {
    double[] Det = new double[2];
    
    Det[0] = 1;
    Det[1] = 0;
    for(int i=0;i<cols && Det[0] != 0;i++) {
      Det[0] *= Math.pow(this.Mat[i][i],2);
      if(Det[0] != 0) 
	DUtil.detNorm(Det);
    }
    DVector V = new DVector(Det);
    return V;
  }
  
  public void inverse() {
    
    //  Compute inverse(R)  //

    this.invertU();

    //  Form inv(R)*trans(inv(R))  //
    
    for(int j=0;j<cols;j++) {
      for(int k=0;k<j;k++) 
	this.axpy(k+1,this.Mat[k][j],1,0,j,this,1,0,k);
      this.scal(j+1,this.Mat[j][j],1,0,j);
    }
  }
}