/*  
 ##  class NPoFull
 ##  
 ##  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.NMatrix;
import edu.rice.linpack.util.*;
import edu.rice.linpack.LNumber.*;

public class NPoFull extends NFull {

  public NPoFull() {
    super();
  }
  public NPoFull(int n, int q) {
    super(n, q);
  }
  public NPoFull(LNumber[][] f) {
    super(f);
  }
  public NPoFull(NPoFull F) {
    super(F);
  }

  public NPoPack pack() {
    NPoPack P = new NPoPack(Mat);
    return P;
  }
  
  public void transpose() {
  }

  public LNumber getElem(int i, int j) {
    if(i > j)
      return super.getElem(j,i);
    return super.getElem(i,j);
  }
  public LNumber[] getRow(int i) {
    return getRow(i,0);
  }
  public LNumber[] getColumn(int i) {
    return getRow(i,0);
  }
  public LNumber[] getColumn(int i, int j) {
    return getRow(i,j);
  }
  public LNumber[] getRow(int i, int j) {
    LNumber[] F = new LNumber[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, LNumber[] F) {
    for(int r=0;r<c;r++) 
      Mat[r][c] = F[r].Clone();
    for(int r=c;r<rows;r++) 
      Mat[r][c] = F[r].Clone();
  }
  public void setElem(int i, int c, LNumber L) {
    if(i > c)
      super.setElem(c,i,L);
    else
      super.setElem(i,c,L);
  }

  protected LNumber oneNorm() {
    LNumber[] Z = new LNumber[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].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 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);
  }
  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++;
  }
  
  void decom(Pul P) 
       throws SingularMatrixException
  {
    LNumber[] W = new LNumber[cols];
    
    for(int k=0;k<cols;k++) {
      // Reduction loop //
      LNumber maxdia = 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(Mat[q][q].greaterThan(maxdia)) {
	    maxdia = this.Mat[q][q];
	    maxln = q;
	  }
	}
      }
      if(!maxdia.greaterOrEqual(0)) 
	throw new SingularMatrixException(k+1);
      else {
	if(k != maxln) {
	  this.swap(k,1,0,k,this,1,0,maxln);
	  Mat[maxln][maxln] = Mat[k][k];
	  Mat[k][k] = maxdia;
	  NUtil.swapElems(pivot, k, maxln);
	}
	W[k] = Mat[k][k].sqrt();
	Mat[k][k] = W[k].Clone();
	for(int j=kp;j<cols;j++) {
	  if(k != maxln) {
	    this.chSwaps(k,j,maxln);
	  }
	  Mat[k][j].divTo(W[k]);
	  W[j] = Mat[k][j];
	  this.axpy(j-k,W[j].negate(),W,1,kp,1,kp,j);
	}
      }
    }
  }
  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);
	  }
	  NUtil.swapElems(pivot,k,P.pu);
	}
	P.pu--;
      }
    }
  }
  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++) {
      LNumber S = Mat[0][0].Clone();
      S.setZero();
      for(int k=0;k<j;k++) {
	LNumber T = Mat[k][j].sub(this.dot(k,1,0,k,this,1,0,j));
	T.divTo(Mat[k][k]);
	Mat[k][j] = T;
	S.addTo(T.square());
      }
      S = Mat[j][j].sub(S);
      if(!S.greaterThan(0)) 
	throw new SingularMatrixException(j+1);
      else 
	Mat[j][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();
    
    //  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 LNumber[]s are frequently rescaled to avoid overflow.
	*/
    //  Solve trans(R)*W = E  //
    
    for(int j=0;j<cols;j++) {
      Z[j] = Mat[0][0].Clone();
      Z[j].setZero();
    }
      
    this.solveTransUW(Z);
    
    LNumber S = (NUtil.asum(cols,Z,1)).invTo();
    NUtil.scal(cols,S,Z,1);

    //  Solve R*Y = W  //
    
    this.solveUZ(Z);
    
    //  Solve trans(R)*V = Y  //
    
    LNumber ynorm = this.solveTransRV(Z);
    
    //  Solve R*Z = V  //
    
    ynorm = this.solveUZ(Z,ynorm);
    
    if(anorm.equals(0))
      return anorm;
   
    ynorm.divTo(anorm);  
    return ynorm;
  }
  void solveUZ(LNumber[] Z) {
    for(int k=cols-1;k>=0;k--) {
      if((Z[k].abs()).greaterThan(Mat[k][k].abs())) {
	LNumber S = (Mat[k][k].abs()).div(Z[k].abs());
	NUtil.scal(cols,S,Z,1);
      } 
      if(!Mat[k][k].equals(0)) 
	Z[k].divTo(Mat[k][k]);
      else 
	Z[k].setOne();
      
      this.axpy(k,Z[k].negate(),1,0,k,Z,1,0);
    }
    LNumber S = (NUtil.asum(cols,Z,1)).invTo();
    NUtil.scal(cols,S,Z,1);
  }
  private LNumber solveTransRV (LNumber[] Z) {
    LNumber ynorm = Z[0].Clone();
    ynorm.setOne();
    for(int k=0;k<cols;k++) {
      Z[k].subTo(this.dot(k,1,0,k,Z,1,0));
      if((Z[k].abs()).greaterThan(Mat[k][k])) {
	LNumber S = Mat[k][k].div(Z[k].abs());
	NUtil.scal(cols,S,Z,1);
	ynorm.multTo(S);
      }
      Z[k].divTo(Mat[k][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 Job) {
    this.solve(B);
  }
  public void solve(LNumber[] B) {
    // Solve Trans(R)*Y = B //
    
    this.solveTransAX(B);

    // Solve R*X = Y  //

    this.solveUX(B);

    return;
  }
  
  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();
    for(int i=0;i<cols;i++) {
      Det[0].multTo(Mat[i][i].square());
      if(!Det[0].equals(0)) 
	NUtil.detNorm(Det);
      else
	return Det;
    }
    return Det;
  }
  
  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);
    }
  }
}