/*##
 ##  class NFull
 ## 
 ##  This class contains the factor, determ, inverse, conditon, and solve 
 ##   for full matrices.  It also contains the Singular Value Decompostion,
 ##   QR-decomposition, and Cholesky updating routines. 
 ##  This class has methods to compute the transpose of a matrix and several
 ##   matrix-matrix and matrix-vector multiplication routines.  
 ##
 ##*/

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

public class NFull extends NMatrix {

  public NFull() {
    super();
    cons();
  }
  public NFull(int i, int j) {
    cons(i,j);
  }
  public NFull(LNumber[][] f) {
    cons(f);
  }

  /*##
   ##  This is the copy constructor.
   ##*/

  public NFull(NFull F) {
    cons(F);
  }

  /*##
   ##  This is the constructor that allows for the manipulation of a 
   ##   subsection of a matrix by creating a new NFull from the subsection.
   ##*/

  public NFull(int ni, int nj, LNumber[][] f, int strow, int stcol) {
    rows = ni;
    cols = nj;
    Mat = new LNumber[rows][cols];
    for(int i=0;i<rows;i++) {
      for(int j=0;j<cols;j++) {
	Mat[i][j] = f[i+strow][j+stcol];
      }
    }
    pivot = new int[cols];
  }

  /*##
   ##  This bands an NFull. 
   ##   
   ##  ml is the number of rows below the diagonal.
   ##  mu is the number of rows above the diagonal.
   ##*/

  public NBanded band(int ml, int mu) {
    NBanded B = new NBanded(ml,mu,cols);
    int M = ml+mu;
    for(int j=0;j<cols;j++) {
      int i1 = Math.max(0,M-j);
      int i2 = Math.min(cols,cols+M-j);
      for(int i=0;i<rows;i++) {
	int k = i-j+M;
	if(k>=i1 && k<i2 && k<B.numofRows()) {
	  if(k>=ml) 
	    B.setElem(k,j,Mat[i][j]);	
	  else {
	    LNumber T = newOfType();
	    T.setZero();
	    B.setElem(k,j,T);
	  }
	}
      }
    }
    return B;
  }

  /*## 
   ##  These routines set the pivot array of an NFull.
   ##  This is necessary for qrPDecompose.
   ##*/

  public void setPivot(int i, int p) {
    pivot[i] = p;
  }
  public void setPivots(int[] p) {
    int q = Math.min(p.length,cols);
    for(int i=0;i<q;i++) 
      pivot[i] = p[i];
  }

  /*## 
   ##  This routine sets an NFull to be its own transpose()
   ##*/

  public void transpose() {
    if(rows == cols) 
      for(int i=0;i<rows;i++) {
	for(int j=0;j<i;j++) {
	  this.swapElem(i,j,j,i);
	}
      }
    else {
      LNumber[][] G = new LNumber[cols][rows];
      for(int i=0;i<rows;i++) {
	for(int j=0;j<cols;j++) {
	  G[j][i] = Mat[i][j];
	}
      }
      Mat = G;
      cols = rows;
      rows = G.length;
    }
  }

  /*##
   ##  This creates a new NFull and sets it as the transpose.
   ##*/

  public void transpose(NFull N) {
    for(int i=0;i<rows;i++) {
      for(int j=0;j<cols;j++) {
	N.setElem(j,i,Mat[i][j]);
      }
    }
  }

  /*##  
   ##  NFull times a NFull multiplication.
   ##*/

  public NFull matMult(NFull N) {
    NFull G = new NFull(rows,N.numofCols());
    for(int i=0;i<rows;i++) {
      LNumber[] L = this.getRow(i,0);
      for(int j=0;j<N.numofCols();j++) {
	G.setElem(i,j,N.dot(cols,1,0,j,L,1,0));
      }
    }
    return G;
  }

  /*##
   ##  NFull-NFull multiplication where calling NFull is set to be 
   ##   the solution.
   ##*/

  public void matMultTo(NFull N) {
    LNumber[][] G = new LNumber[rows][N.numofCols()];
    for(int i=0;i<rows;i++) {
      LNumber[] L = this.getRow(i,0);
      for(int j=0;j<N.numofCols();j++) {
	G[i][j] = N.dot(cols,1,0,j,L,1,0);
      }
    }
    Mat = G;
    cols = N.numofCols();
  }
  public NFull matMultRev(NFull N) {
    NFull G = new NFull(N.numofRows(),cols);
    for(int i=0;i<N.numofRows();i++) {
      LNumber[] L = N.getRow(i,0);
      for(int j=0;j<cols;j++) {
	G.setElem(i,j,this.dot(rows,1,0,j,L,1,0));
      }
    }
    return G;
  }

  /*##
   ##  NFull times a vector.
   ##*/

  public LNumber[] matVec(LNumber[] N) {
    LNumber[] G = new LNumber[rows];
    for(int i=0;i<rows;i++) {
      LNumber[] L = this.getRow(i,0);
      G[i] = NUtil.dot(cols,L,1,0,N,1,0);
    }
    return G;
  }

  /*##
   ##  This takes the 1-Norm of an NFull.  This is necessary for condition.
   ##*/

  private LNumber oneNorm() {
    LNumber a = this.asum(cols,1,0,0);
    for(int j=1;j<cols;j++) 
      a.maxTo(this.asum(cols,1,0,j));
    return a;
  }
  
  /*##  
   ##  svDecompose is a subroutine to reduce a single precision rows X cols 
   ##  NFull by orthogonal transformations U and V to diagonal form.  
   ##  The diagonal elements S[i] are the singular values of this 
   ##  NFull.  The columns of U are the corresponding left singular
   ##  vectors, and the columns of V the right singular vectors.
   ##
   ##  On Entry:
   ##
   ##     this      The NFull whose singular value decomposition is to be 
   ##               computed.  This NFull is destroyed by svDecompose.
   ##
   ##     job       int
   ##               job controls the computation of the singular
   ##               vectors.  It has the decimal expansion AB
   ##               with the following meaning:
   ##
   ##                    A == 0    Do not compute the left singular
   ##                              vectors.
   ##                    A == 1    Return the n left singular vectors
   ##                              in U.
   ##                    A >= 2    Return the first min(rows,cols) singluar
   ##                              vectors in U.
   ##                    B == 0    Do not compute the right singular
   ##                              vectors.
   ##                    B == 1    Return the right singular vectors
   ##                              in V.
   ##
   ##  On Return:
   ##
   ##     S         LNumber[MM], where MM = Math.min(rows+1,cols).    
   ##               The first min(rows,cols) entries of S contain the 
   ##               singular values of this NFull arranged in 
   ##               descending order of magnitude.
   ##
   ##     E         LNumber[cols]
   ##               E ordinarily contains zeros.  However see the 
   ##               discussion of info for exceptions.
   ##
   ##     U         LNumber[rows][k], where if job(A) == 1 then k == rows
   ##                                       else k == Min(rows,cols).
   ##               U contains the NFull of right singular vectors.
   ##               U is not referenced if job(A) == 0.  
   ##
   ##     V         LNumber[cols][cols]
   ##               V contains the NFull of right singular vectors.
   ##               V is not referenced if job == 0.  
   ##
   ##     info      int
   ##               The singular values (and their corresponding
   ##               singular vectors) S(info),S(info+1),...,S(M-1)
   ##               are correct (here M=min(rows,cols)).  Thus if 
   ##               info == 0, all the singular values and their 
   ##               vectors are correct.  In any event, the matrix
   ##               B = U.transpose()*this*V is the bidiagonal matrix
   ##               with the elements of S on its diagonal and the 
   ##               elements of E on its super-diagonal.  Thus the 
   ##               singular values of this and B are the same.
   ##*/


  public int svDecompose(LNumber[] S, LNumber[] E, NFull U, NFull V, 
			 int job) {
    
    LNumber[] W = new LNumber[rows];
    int info = 0;
    
    /*##  
      ##  This is the maximum number of iterations.
      ##*/

    int MAXIT = 30;

    /*##  
      ##  Determine what is to be computed
      ##*/

    boolean wantU = false;
    boolean wantV = false;
    int jobu = (job%100)/10;
    int NCU = rows;
    if(jobu > 1) 
      NCU = Math.min(cols,rows);
    if(jobu != 0) 
      wantU = true;
    if(job%10 != 0) 
      wantV = true;

    /*##
     ##  Reduce Mat to bidiagonal form, storing the diagonal elements
     ##   in S and the super-diagonal elements in E.
     ##*/

    int NCT = Math.min(rows-1,cols);
    int NRT = Math.max(0,Math.min(cols-2,rows));
    int lu = Math.max(NCT,NRT);
    for(int q=0;q<lu;q++) {
      //System.out.println("++++++++++++++++++++++++++++++++++++++++");
      //NUtil.Print(S);
      int qp = q+1;
      if(q < NCT) { 
	
	/*##  
	 ##  Compute the transformation for the q-th column and place the
	 ##   q-th diagonal in S[q].
	 ##*/
	//this.Print();
	S[q] = nrm2(rows-q,1,q,q);
	//System.out.println(S[q]+" norm2");
	if(!S[q].equals(0)) {
	  if(!Mat[q][q].equals(0)) 
	    S[q] = NUtil.signOfM(S[q],Mat[q][q]);
	  this.scal(rows-q,S[q].inv(),1,q,q);
	  Mat[q][q].pplus();
	}
	S[q].negateTo();
	//System.out.println(S[q]+" negated");
	//this.Print();
      }
      for(int j=qp;j<cols;j++) {
	if(q < NCT) {
	  if(!S[q].equals(0)) {

	    /*##
	     ##  Apply the transformation
	     ##*/
	    if(q==0 /*&& j==qp*/) {
	      //this.Print();
	    }
	    //System.out.println((rows-q)+" "+q+" "+j+" "+Mat[q][q]);
	    LNumber T = (dot(rows-q,1,q,q,this,1,q,j).div(Mat[q][q]))
	      .negateTo();
	    //System.out.println((dot(rows-q,1,q,q,this,1,q,j)));
	   //System.out.println((dot(rows-q,1,q,q,this,1,q,j).div(Mat[q][q])));
	    //System.out.println("T "+T);
	    this.axpy(rows-q,T,1,q,q,this,1,q,j);
	    if(q==0 && false) {
	      //this.Print();
	      //System.out.println("VARS "+T+" "+S[q]+" "+j);
	    }
	  }
	}
	
	/*##
	 ##  Place the q-th row of Mat into E for the subsequent
	 ##   calculation of the row transformation.
	 ##*/

	// This should be the CONJ

	E[j] = Mat[q][j].conj();
      }
      //if(q==0)
      //this.Print();
      if(wantU && q<NCT) {

	/*##
	 ##  Place the transformation in U for subsequent back 
	 ##   multiplication.
	 ##*/
	
	for(int i=q;i<rows;i++) 
	  U.setElem(i,q,Mat[i][q]);
      }
      if(q < NRT) {

	/*##  
	 ##  Compute the q-th row transformation and place the q-th
	 ##   super-diagonal in E[q]  
	 ##*/

	E[q] = NUtil.nrm2(cols-qp,E,1,qp);
	if(!E[q].equals(0)) {
	  if(!E[qp].equals(0)) 
	    E[q] = NUtil.signOfM(E[q],E[qp]);
	  NUtil.scal(cols-qp,E[q].inv(),E,1,qp);
	  E[qp].pplus();
	}

	// This should be the negative CONJ

	E[q] = (E[q].conj()).negateTo();
	if(qp < rows && !E[q].equals(0)) {
	  
	  /*## 
	   ##  Apply the transformation.
	   ##*/

	  for(int i=qp;i<rows;i++) {
	    W[i] = E[q].Clone();
	    W[i].setZero();
	  }
	  for(int j=qp;j<cols;j++) 
	    this.axpy(rows-qp,E[j],1,qp,j,W,1,qp);
	  for(int j=qp;j<cols;j++) 
	    this.axpy(rows-qp,((E[j].div(E[qp])).negateTo()).conj()
		      ,W,1,qp,1,qp,j);
	}
	if(wantV) {

	  /*##
	   ##  Place the transformation in V for subsequent back
	   ##   multiplication.
	   ##*/
	  
	  for(int i=qp;i<cols;i++) 
	    V.setElem(i,q,E[i].Clone());
	}

      }
      //System.out.println(S[q]);
    }

    /*##  
     ##  Set up the final bidiagonal matrix of order M.
     ##*/

    int M = Math.min(rows+1,cols);
    if(NCT < cols) 
      S[NCT] = Mat[NCT][NCT];
    if(rows < M) {
      S[M-1] = newOfType();
      S[M-1].setZero();
    }
    if(NRT+1 < M) 
      E[NRT] = Mat[NRT][M-1].Clone();
    E[M-1].setZero();

    /*##
     ##  If required generate U and/or V.
     ##*/

    //NUtil.Print(S);
    //NUtil.Print(E);
    /*    
    NUtil.Print(S);
    System.out.println("++++++++++++++++++++++++++++");
    NUtil.Print(E);
    System.out.println("++++++++++++++++++++++++++++");
    this.Print();
    */

    if(wantU) 
      U.generateU(NCT,NCU,S);
    if(wantV) 
      V.generateV(NRT,E);
    /*
    NUtil.Print(S);
    System.out.println("++++++++++++++++++++++++++++");
    NUtil.Print(E);
    System.out.println("++++++++++++++++++++++++++++");
    this.Print();
    */

    //NUtil.Print(S);
    //NUtil.Print(E);

    // Major problem still right here - code for complex, but not real
    // This should work for real as well, but it adds computation.

    //System.out.println("DUMBASS you forgot to fix this!!!");

    for(int i=0;i<M;i++) {
      if(!(S[i].abs()).equals(0)) {
	LNumber T = S[i].modulus();
	LNumber R = S[i].div(T);
	S[i] = T;
	if(i < M-1)
	  E[i].divTo(R);
	if(wantU)
	  U.scal(rows,R,1,0,i);
      }
      if(i < M-1) {
	if(!(E[i].abs()).equals(0)) {
	  LNumber T = E[i].modulus();
	  LNumber R = T.div(E[i]);
	  E[i] = T;
	  S[i+1].multTo(R);
	  if(wantV)
	    V.scal(cols,R,1,0,i+1);
	}
      }
    }
    /*    
    NUtil.Print(S);
    System.out.println("++++++++++++++++++++++++++++");
    NUtil.Print(E);
    System.out.println("++++++++++++++++++++++++++++");
    this.Print();
    */
    
    /*##  
     ##  Main iteration loop for the singular values.
     ##*/
    
    int iter = 0;
    int MM = M;
    while(M>0) {

      /*##  
       ##  If too many iterations have been performed, return error flag.
       ##*/

      if(iter >= MAXIT) 
	return M;
      else {

	/*##
	 ##  qgetter determines what to do next by searching for negligible 
	 ##   elements in S and E.
	 ##*/

	int q = qgetter(M-2,S,E);
	if(q == M-2) {

	  /*##
	   ##  E[M-2] is negligible (convergence).
	   ##*/
	  
	  q++;
	  
	  /*## 
	   ##  Make the singular value positive.
	   ##*/
	  
	  if(S[q].lessThan(0)) {
	    S[q].negateTo();
	    if(wantV) 
	      V.scal(cols,-1,1,0,q);
	  }
	  
	  /*##
	   ##  Order the singular value.
	   ##*/

	  boolean notDone = true;
	  while(q < MM-1 && notDone) {
	    if(S[q].lessThan(S[q+1])) {
	      NUtil.swapElems(S,q,q+1);
	      if(wantV && q < cols-1) 
		V.swap(cols,1,0,q,V,1,0,q+1);
	      if(wantU && q < rows-1) 
		U.swap(rows,1,0,q,U,1,0,q+1);
	      q++;
	    }
	    else 
	      notDone = false;
	  }
	  iter = 0;
	  M--;
	}
	else {
	  
	  /*## 
	   ##  If there is not convergence, findNextStep helps qgetter 
	   ##   determine what to do next.
	   ##*/

	  int qs = this.findNextStep(q,M,S,E);
	  if(qs == q) {
	    
	    /*##
	     ##  E[q-1] is negligible, q < M, and S[q], ... , S[M-1] are
	     ##   not negligible (qr step).
	     ##  Perform one qr shift.
	     ##*/

	    q++;

	    /*##
	     ##  Calculate the shift.
	     ##*/

	    NTrig Tg = this.fggetter(q,M,S,E);

	    /*##
	     ##  Chase zeros.
	     ##*/
	    
	    for(int k=q;k<M-1;k++) {
	      NUtil.rotg(Tg);
	      if(k != q) 
		E[k-1] = Tg.getA();

	      // Problem right here

	      Tg.setA(((Tg.getCos()).mult(S[k].realPt()))
		      .add((Tg.getSin()).mult(E[k].realPt())));
	      E[k] = ((Tg.getCos()).mult(E[k])).sub((Tg.getSin()).mult(S[k]));
	      Tg.setB((Tg.getSin()).mult(S[k+1].realPt()));
	      S[k+1].multTo(Tg.getCos());
	      if(wantV) 
		V.rot(cols,1,0,k,V,1,0,k+1,Tg);
	      NUtil.rotg(Tg);
	      S[k] = Tg.getA();
	      Tg.setA(((E[k].realPt()).mult(Tg.getCos()))
		      .add((S[k+1].realPt()).mult(Tg.getSin())));
	      S[k+1] = (E[k].mult((Tg.getSin()).negate()))
		.add(S[k+1].mult(Tg.getCos()));
	      Tg.setB((Tg.getSin()).mult(E[k+1].realPt()));
	      E[k+1].multTo(Tg.getCos());
	      if(wantU && k < rows-1) 
		U.rot(rows,1,0,k,U,1,0,k+1,Tg);
	    }
	    E[M-2] = Tg.getA();
	    iter++;
	  }
	  else if(qs == M-1) {

	    /*## 
	     ##  S[M] and E[q-1] are negligible and q < M.
	     ##  Deflate negligible S[M].
	     ##*/

	    q++;
	    int MMOne = M-2;
	    NTrig Tg = new NTrig();
	    Tg.setB(E[MMOne].realPt());
	    E[MMOne].setZero();
	    for(int k=MMOne;k>=q;k--) {

	      // S[k] should only be real - is this a prob?

	      Tg.setA(S[k].realPt());
	      NUtil.rotg(Tg);
	      S[k] = Tg.getA();
	      if(k != q) {

		// Need it here, too

		Tg.setB(((Tg.getSin()).mult(E[k-1].realPt())).negateTo());
		E[k-1].multTo(Tg.getCos());
	      }
	      if(wantV) 
		V.rot(cols,1,0,k,V,1,0,MMOne+1,Tg);
	    }
	  }
	  else {
	    
	    /*## 
	     ##  S[q] is negligible and q < M.
	     ##  Split at negligible S[q].
	     ##*/

	    q = qs;
	    NTrig Tg = new NTrig();

	    // Real part of E[q] wanted here, too
	    
	    Tg.setB(E[q].realPt());
	    E[q].setZero();
	    q++;
	    for(int k=q;k<M;k++) {

	      // Here problem, too

	      Tg.setA(S[k].realPt());
	      NUtil.rotg(Tg);
	      S[k] = Tg.getA();
	      Tg.setB(((Tg.getSin()).mult(E[k].realPt())).negateTo());
	      E[k].multTo(Tg.getCos());
	      if(wantU) 
		U.rot(rows,1,0,k,U,1,0,q-1,Tg);
	    }
	  }
	}
      }
    }
    return info;
  }
  private void generateU(int NCT, int NCU, LNumber[] S) {
    for(int j=NCT;j<NCU;j++) {
      for(int i=0;i<cols;i++) {
	Mat[i][j] = S[0].Clone();
	Mat[i][j].setZero();
      }
      Mat[j][j].setOne();
    }
    if(NCT > 0) {
      for(int q=NCT-1;q>=0;q--) {
	if(!S[q].equals(0)) {
	  for(int j=q+1;j<NCU;j++) {
	    LNumber T = (dot(cols-q,1,q,q,this,1,q,j).div(Mat[q][q]))
	      .negateTo();
	    this.axpy(cols-q,T,1,q,q,this,1,q,j);
	  }
	  this.scal(cols-q,-1,1,q,q);
	  this.Mat[q][q].pplus();
	  for(int i=0;i<q;i++) { 
	    Mat[i][q] = S[0].Clone();
	    this.Mat[i][q].setZero();
	  }
	}
	else {
	  for(int i=0;i<cols;i++) {
	    Mat[i][q] = S[0].Clone();
	    this.Mat[i][q].setZero();
	  }
	  this.Mat[q][q].setOne();
	}
      }
    }
  }
  private void generateV(int NRT, LNumber[] E) {
    for(int q=cols-1;q>=0;q--) {
      int qp = q+1;
      if((q < NRT) && (!E[q].equals(0))) {
	for(int j=qp;j<rows;j++) {
	  LNumber T = (dot(rows-qp,1,qp,q,this,1,qp,j).div(Mat[qp][q]))
		       .negateTo();
	  this.axpy(rows-qp,T,1,qp,q,this,1,qp,j);
	}	    
      }
      for(int i=0;i<rows;i++) {
	Mat[i][q] = E[0].Clone();
	Mat[i][q].setZero();
      }
      Mat[q][q].setOne();
    }
  }
  private int qgetter(int st, LNumber[] S, LNumber[] E) {
    for(int q=st;q >= 0;q--) {
      LNumber Test = (S[q].modulus()).add(S[q+1].modulus());
      LNumber ZTest = Test.add(E[q].modulus());
      if(ZTest.equals(Test)) {
	E[q].setZero();
	return q;
      }
    }
    return -1;
  }
  private int findNextStep(int q, int M, LNumber[] S, LNumber[] E) {
    for(int qs=M-1;qs>q;qs--) {
      LNumber Test = E[qs].modulus();
      if(qs == M-1) 
	Test.setZero();
      if(qs != q+1) 
	Test.addTo(E[qs-1].modulus());
      LNumber Ztest = Test.add(S[qs].modulus());
      if(Ztest.equals(Test)) {
	S[qs].setZero();
	return qs;
      }
    }
    return q;
  }
  private NTrig fggetter(int q, int M, LNumber[] S, LNumber[] E) {
    NTrig Tg = new NTrig();
    LNumber scale = NUtil.mMax(S[M-1],S[M-2],E[M-2],S[q],E[q]);

    // These are all supposed to be dreal

    LNumber SM = (S[M-1].realPt()).div(scale);
    LNumber SMM = (S[M-2].realPt()).div(scale);
    LNumber EMM = (E[M-2].realPt()).div(scale);
    LNumber Sq = (S[q].realPt()).div(scale);
    LNumber Eq = (E[q].realPt()).div(scale);

    // It isn't a problem until here

    LNumber B = (((SMM.add(SM)).mult(SMM.sub(SM))).add(EMM.square())).div(2);
    LNumber C = (SM.mult(EMM)).square();
    LNumber shift = C.Clone();
    shift.setZero();
    if(!(B.equals(0) && C.equals(0))) {
      shift = ((B.square()).add(C)).sqrt();
      if(!B.greaterOrEqual(0)) 
	shift.negateTo();
      shift = C.div(B.add(shift));
    }
    Tg.setA(((Sq.add(SM)).mult(Sq.sub(SM))).sub(shift));
    Tg.setB(Sq.mult(Eq));
    
    return Tg;
  }

  /*##
   ##  qrDecompose uses Householder transformations to compute the
   ##  QR factorization of an n by p NFull.  Column pivoting
   ##  based on the 2-norms of the reduced columns may be
   ##  performed at the user's option.
   ##
   ##  On Entry:
   ##
   ##    this    NFull - Mat[L][p], where L >= n            
   ##            NFull contains the NFull whose decomposition is to 
   ##            be computed.
   ##
   ##
   ##    n       int 
   ##            n is the number of rows of NFull.Mat to be used in
   ##            the calculation
   ##
   ##    p       int 
   ##            p is the number of columns of NFull.Mat to be used
   ##            in the calculation
   ##
   ##    pivot   int[p]
   ##            pivot contains ints that control the selection
   ##            of the pivot columns.  The k-th column NFull.Mat[0][k]
   ##            is placed in one of three classes according to the 
   ##            value of pivot[k]:
   ##
   ##               if pivot[k] > 0, then NFull.getColumn(k,0) is an
   ##                                  initial column.
   ##
   ##               if pivot[k] == 0, then NFull.getColumn(k,0) is a
   ##                                  free column.
   ##
   ##               if pivot[k] < 0, then NFull.getColumn(k,0) is a 
   ##                                  final column.
   ##
   ##            Before the decomposition is computed, initial columns
   ##            are moved to the beginning of NFull.Mat and final
   ##            columns to the end.  Both initial and final columns
   ##            are frozen in place during the computation and only 
   ##            free columns are moved.  At the k-th stage of the 
   ##            reduction, if NFull.getColumn(k,0) is occupied by a
   ##            free column it is interchanged with the free column of
   ##            largest reduced norm.  pivot is not referenced if 
   ##            job == 0.
   ##
   ##
   ##  On Return
   ##
   ##    this    NFull.Mat contains in its upper triangle the upper
   ##            triangular matrix R of the QR factorization.
   ##            Below its diagonal NFull.Mat contains information 
   ##            from which the orthogonal part of the decomposition
   ##            can be recovered.  Note that if pivoting has
   ##            been requested, the decomposition is not that
   ##            of the original NFull.Mat, but that of NFull.Mat
   ##            with its columns permuted as described by pivot.
   ##
   ##    QRaux   LNumber[p]
   ##            QRaux contains further information required to recover 
   ##            the orthogonal part of the decomposition.
   ##
   ##    pivot   pivot[k-1] contains the index of the column of the
   ##            original matrix that has been interchanged into 
   ##            the k-th column, if pivoting was requested.
   ##*/

  public void qrPDecompose(LNumber[] QRaux) {
    Pul P = new Pul(0,-1);
    this.qrPivot(P);
    this.qrDecom(P, QRaux);
  }
  private void qrPivot(Pul P) {
    this.swaperSki(P);
    P.pu = cols-1;
    for(int j=P.pu;j>=0;j--) {
      if(pivot[j] < 0) {
	pivot[j] = -pivot[j];
	if(j != P.pu) {
	  this.swap(rows,1,0,P.pu,this,1,0,j);
	  NUtil.swapElems(pivot,P.pu,j);
	}
	P.pu--;
      }
    }
  }
  void swaperSki(Pul P) {
    for(int j=0;j<cols;j++) {
      boolean swapj = (pivot[j] > 0);
      if(pivot[j] < 0) 
	pivot[j] = -j-1;
      else 
	pivot[j] = j+1;

      if(swapj) 
	this.swap2(j,P);
    }
  }
  void swap2(int j, Pul P) {
    if(j != P.pl) 
      this.swap(rows,1,0,P.pl,this,1,0,j);
    pivot[j] = pivot[P.pl];
    pivot[P.pl] = j+1;
    P.pl++;
  }
   
  public void qrDecompose(LNumber[] QRaux) {
    Pul P = new Pul(0,-1);
    qrDecom(P,QRaux);
  }
  private void qrDecom(Pul P, LNumber[] QRaux) {
    LNumber[] Work = new LNumber[cols];

    for(int j=P.pl;j<=P.pu;j++) {
      QRaux[j] = this.nrm2(rows,1,0,j);
      Work[j] = QRaux[j].Clone();
    }
    int qup = Math.min(rows,cols);
    for(int q=0;q<qup;q++) {
      if(q >= P.pl && q<P.pu) {
	LNumber maxnorm = QRaux[q];
	int maxj = q;
	for(int j=q;j<=P.pu;j++) {
	  if(QRaux[j].greaterThan(maxnorm)) {
	    maxnorm = QRaux[j].realPt();
	    maxj = j;
	  }
	}
	if(maxj != q) {
	  this.swap(rows,1,0,q,this,1,0,maxj);
	  QRaux[maxj] = QRaux[q].Clone();
	  Work[maxj] = Work[q].Clone();
	  NUtil.swapElems(pivot,maxj,q);
	}
      }
      QRaux[q] = newOfType();
      QRaux[q].setZero();
      if(q < rows-1) {
	LNumber normxq = this.nrm2(rows-q,1,q,q);
	if(!normxq.equals(0)) {
	  if(!this.Mat[q][q].equals(0)) 
	    normxq = NUtil.signOfM(normxq,Mat[q][q]);
	  this.scal(rows-q,normxq.inv(),1,q,q);
	  this.Mat[q][q].pplus();
	  for(int j=q+1;j<cols;j++) {
	    LNumber T = ((dot(rows-q,1,q,q,this,1,q,j)).div(Mat[q][q]))
	      .negateTo();
	    this.axpy(rows-q,T,1,q,q,this,1,q,j);
	    if(j >= P.pl && j <= P.pu) {
	      if(!QRaux[j].equals(0)) {

		// In this area are 3 dreals that I do not use

		LNumber TT = (((Mat[q][j].modulus()).div(QRaux[j].realPt()))
			      .square()).subFrom(1);
		if(!TT.greaterThan(0)) 
		  TT.setZero();
		T = TT.Clone();
		TT = ((TT.mult(((QRaux[j].realPt()).div(Work[j].realPt()))
			       .square())).div(2)).add(1);
		if(!TT.equals(1)) {
		  //System.out.println("FUNKY SQRT");
		  QRaux[j].multTo(T.sqrt());
		}
		else {
		  QRaux[j] = this.nrm2(rows-q-1,1,q+1,j);
		  Work[j] = QRaux[j].Clone();
		}
	      }
	    }
	  }
	  QRaux[q] = Mat[q][q];
	  Mat[q][q] = normxq.negate();
	}
      }
    }
  }

  /*##
   ##  The following qr*** routines apply the output of qrDecompose to 
   ##  compute coordinate transformations, projections, and least squares
   ##  solutions.
   ##  for k < Min(n,p), let XK be the matrix
   ##
   ##        XK = (this.getColumn(pivot(0)-1,0),
   ##                 this.getColumn(pivot(1)-1,0),
   ##                 ... ,this.getColumn(pivot(k)-1,0));
   ##
   ##  formed from columns pivot[0], ... , pivot[k] of the original
   ##  n X p matrix NFull.Mat that was input to qrDecompose (if no
   ##  pivoting was done, XK consists of the first k+1 columns of 
   ##  NFull.Mat in their original order).  qrDecompose produces a 
   ##  factored orthogonal matrix Q and an upper triang. matrix R s.t.
   ##
   ##          XK = Q * (R)
   ##                   (0)
   ##
   ##  this information is contained in coded form in NFull.Mat and QRaux
   ##
   ##  On Entry
   ##
   ##    this   NFull.Mat[n][p]
   ##           NFull.Mat contains the output of qrDecompose
   ##
   ##    n      int 
   ##           n is the number of rows of the matrix XK.  It must
   ##           have the same value as n in qrDecompose.
   ##
   ##    k      int 
   ##           k is the number of columns of the matrix XK.  k
   ##           must not be greater than min(n,p), where p is the 
   ##           same as in qrDecompose.
   ##
   ##    QRaux  LNumber[p]
   ##           QRaux contains the auxiliary output from qrDecompose.
   ##
   ##    Y      LNumber[n]
   ##           Y contains the vector that is to be manipulated
   ##           by qr***.
   ##
   ##  On Return from its respective function
   ##
   ##  qrQy:
   ##    Qy     LNumber[n]
   ##           Qy contains Q*Y
   ##
   ##  qrQty:
   ##    Qty    LNumber[n]
   ##           Qty contains Q.transpose()*Y
   ##
   ##  qrQB:
   ##    B      LNumber[k]
   ##           B contains the solution of the least squares problem
   ##
   ##                minimize norm2(Y - XK*B),
   ##
   ##           (Note that if pivoting was requested in qrDecompose, 
   ##           the j-th component of B will be associated with column
   ##           pivot[j-1]-1 of the original NFull.Mat that was input
   ##           into qrDecompose.)
   ##
   ##    info   int
   ##           info is zero unless R is exactly singular in the 
   ##           computation of B.  In this case, info is the index plus
   ##           one of the first zero diagonal element of R and B is 
   ##           left unaltered.
   ##
   ##  qrQRSD:
   ##    Rsd    LNumber[n]
   ##           Rsd contains the least squares residual Y - XK*B,
   ##           Rsd is also the orthogonal projection of Y onto the
   ##           orthogonal complement of the column space of XK.
   ##
   ##  qrXB:
   ##    Xb     LNumber[n]
   ##           Xb contains the least squares approximation XK*B,
   ##           Xb is also the orthogonal projection of Y onto the 
   ##           column space of NFull.Mat.
   ##*/
  
  public void qrQy(LNumber[] QRaux, LNumber[] Y, LNumber[] Qy) {
    
    int k = Math.min(rows,cols);
    int ju = Math.min(k,rows-1) - 1;
    
    if(ju < 0) 
      Qy[0] = Y[0].Clone();
    else {
      NUtil.copy(rows,Y,1,0,Qy,1,0);
      for(int j=ju;j>=0;j--) 
	this.qrsolver(j,Qy,QRaux);
    }
  }
  
  public void qrQty(LNumber[] QRaux, LNumber[] Y,LNumber[] Qty) {
    
    int k = Math.min(rows,cols);
    int ju = Math.min(k,rows-1) - 1;
    
    if(ju < 0) 
      Qty[0] = Y[0].Clone();
    else {
      NUtil.copy(rows,Y,1,0,Qty,1,0);
      for(int j=0;j<=ju;j++) 
	this.qrsolver(j,Qty,QRaux);
    }
  }

  public int qrB(LNumber[] QRaux, LNumber[] Y,
		 LNumber[] Qty, LNumber[] B) {

    int k = Math.min(cols,rows);
    int ju = Math.min(k,rows-1) - 1;
    
    qrQty(QRaux,Y,Qty);

    if(ju < 0) {
      if(Mat[0][0].equals(0)) 
	return 1;
      else 
	B[0] = Y[0].div(Mat[0][0]);
    }
    else {
      NUtil.copy(k,Qty,1,0,B,1,0);
      for(int j=k-1;j>=0;j--) {
	if(Mat[j][j].equals(0)) 
	  return j+1;
	else {
	  B[j].divTo(Mat[j][j]);
	  if(j > 0) 
	    this.axpy(j,B[j].negate(),1,0,j,B,1,0);
	}
      }
    }
    return 0;
  }

  public void qrXB(LNumber[] QRaux, LNumber[] Y,
		  LNumber[] Qty, LNumber[] Xb) {
    
    int k = Math.min(cols,rows);
    int ju = Math.min(k,rows-1) - 1;
    
    qrQty(QRaux,Y,Qty);

    if(ju < 0) 
      Xb[0] = Y[0].Clone();
    else {
      NUtil.copy(k,Qty,1,0,Xb,1,0);
      for(int i=k;i<rows;i++) {
	Xb[i] = Y[0].Clone();
	Xb[i].setZero();
      }
      for(int j=ju;j>=0;j--) 
	this.qrsolver(j,Xb,QRaux);
    }
  }
  public void qrRSD(LNumber[] QRaux, LNumber[] Y,
		    LNumber[] Qty, LNumber[] Rsd) {

    int k = Math.min(cols,rows);
    int ju = Math.min(k,rows-1) - 1;
    
    qrQty(QRaux,Y,Qty);
    
    if(ju < 0) {
      Rsd[0] = Y[0].Clone();
      Rsd[0].setZero();
    }
    else {
      int kp = k+1;
      if(k < rows) 
	NUtil.copy(rows-k,Qty,1,k,Rsd,1,k);
      for(int i=0;i<k;i++) {
	Rsd[i] = Y[0].Clone();
	Rsd[i].setZero();
      }
      for(int j=ju;j>=0;j--) 
	this.qrsolver(j,Rsd,QRaux);
    }
  }
  private void qrsolver(int j, LNumber[] R, LNumber[] QRaux) {
    if(!QRaux[j].equals(0)) {
      LNumber Temp = Mat[j][j].Clone();
      Mat[j][j] = QRaux[j];
      LNumber T = (dot(rows-j,1,j,j,R,1,j).div(Mat[j][j])).negateTo();
      this.axpy(rows-j,T,1,j,j,R,1,j);
      Mat[j][j] = Temp;
    }
  }


  public void chUpdate(LNumber[] X,LNumber[][] Z,LNumber[] Y, LNumber[] Rho, 
		       LNumber[] C, LNumber[] S) {
    
    // Update the matrix //

    for(int j=0;j<cols;j++) {
      LNumber Xj = X[j].Clone();
      for(int i=0;i<j;i++) {
	LNumber T = (C[i].mult(Mat[i][j])).add(S[i].mult(Xj));
	Xj = (C[i].mult(Xj)).sub(S[i].mult(Mat[i][j]));
	Mat[i][j] = T;
      }
      NTrig Tg = new NTrig(Xj,C[j],S[j]);
      this.rotg(j,j,Tg);
      C[j] = Tg.getCos();
      S[j] = Tg.getSin();
    }
    for(int j=0;j<Z[0].length;j++) {
      LNumber Zeta = Y[j].Clone();
      for(int i=0;i<cols;i++) {
	LNumber T = (C[i].mult(Z[i][j])).add(S[i].mult(Zeta));
	Zeta = (C[i].mult(Zeta)).sub(S[i].mult(Z[i][j]));
	Z[i][j] = T;
      }
      LNumber azeta = Zeta.modulus();
      if(!azeta.equals(0) && Rho[j].greaterOrEqual(0)) {
	LNumber scale = azeta.mult(Rho[j]);
	Rho[j] = scale.mult((((azeta.div(scale)).square())
			     .add((Rho[j].div(scale)).square())).sqrt());
      }
    }
  }
  public int chDowndate(LNumber[] X,LNumber[][] Z,
		      LNumber[] Y,LNumber[] Rho,LNumber[] C,LNumber[] S) {
    
    int info = 0;
    LNumber Zeta;

    S[0] = X[0].div(Mat[0][0]);
    for(int j=1;j<cols;j++) {
      S[j] = X[j].sub(dot(j,1,0,j,S,1,0));
      S[j].divTo(Mat[j][j]);
    }
    LNumber norm = NUtil.nrm2(cols,S,1,0);
    if(norm.greaterOrEqual(1)) 
      info = -1;
    else {
      LNumber alpha = ((norm.square()).subFrom(1)).sqrt();
      for(int i=cols-1;i>=0;i--) {
	LNumber scale = alpha.add(S[i].modulus());
	LNumber A = alpha.div(scale);
	LNumber B = S[i].div(scale);
	norm = (A.square().add(B.square())).sqrt();
	C[i] = A.div(norm);
	S[i] = B.div(norm);
	alpha = scale.mult(norm);
      }
      for(int j=0;j<cols;j++) {
	LNumber XX = alpha.Clone();
	XX.setZero();
	LNumber T;
	for(int i=j;i>=0;i--) {
	  T = (C[i].mult(XX)).add(S[i].mult(Mat[i][j]));
	  Mat[i][j] = (C[i].mult(Mat[i][j])).sub(S[i].mult(XX));
	  XX = T;
	}
      }
      for(int j=0;j<Z[0].length;j++) {
	Zeta = Y[j];
	for(int i=0;i<cols;i++) {
	  Z[i][j] = (Z[i][j].sub(S[i].mult(Zeta))).div(C[i]);
	  Zeta = (C[i].mult(Zeta)).sub(S[i].mult(Z[i][j]));
	}
	LNumber azeta = Zeta.modulus();
	if(azeta.greaterThan(Rho[j])) {
	  info = 1;
	  Rho[j].setOne();
	  Rho[j].negateTo();
	}
	else 
	  Rho[j].multTo((((azeta.div(Rho[j])).square()).subFrom(1)).sqrt());
      }
    }
    return info;
  }

  public void chExchange(int k, int q, NFull Z,
		       LNumber[] C, LNumber[] S, int job) {
    
    /* k and q are java (0 based) indices */

    int km = k-1;
    int kp = k+1;
    int qmkp = q-k;
    int qm = q-1;
    
    if(job == 1) {
      for(int i=0;i<=q;i++) 
	S[i] = Mat[q-i][q].Clone();
      for(int j=qm;j>=k;j--) {
	for(int i=0;i<=j;i++) 
	  Mat[i][j+1] = Mat[i][j].Clone();
	Mat[j+1][j+1].setZero();
      }
      for(int i=0;i<k;i++) 
	Mat[i][k] = S[q-i].Clone();
      LNumber T = S[0].Clone();
      NTrig Tg = new NTrig();
      for(int i=0;i<qmkp;i++) {
	Tg.setA(S[i+1]);
	Tg.setB(T);
	NUtil.rotg(Tg);
	S[i+1] = Tg.getA();
	T = S[i+1].Clone();
	C[i] = (Tg.getCos()).Clone();
	S[i] = Tg.getSin().Clone();
      }
      Mat[k][k] = T.Clone();
      for(int j=kp;j<cols;j++) {
	int iq = Math.max(0,q-j);
	for(int ii = iq;ii<qmkp;ii++) {
	  int i = qm-ii;
	  this.applyTrans(j,i,ii,C,S);
	}
      }
      for(int j=0;j<Z.numofCols();j++) 
	for(int ii=0;ii<qmkp;ii++) {
	  int i = qm - ii;
	  Z.applyTrans(j,i,ii,C,S);
	}
    }
    else {
      for(int i=0;i<=k;i++) 
	S[qmkp+i] = Mat[i][k].Clone();
      for(int j=k;j<q;j++) {
	for(int i=0;i<=j;i++) 
	  Mat[i][j] = Mat[i][j+1].Clone();
	S[j-k] = Mat[j+1][j+1].Clone();
      }
      for(int i=0;i<=k;i++) 
	Mat[i][q] = S[qmkp+i].Clone();
      for(int i=kp;i<=q;i++) 
	Mat[i][q].setZero();
      for(int j=k;j<cols;j++) {
	if(j != k) {
	  int iu = Math.min(j-1,qm);
	  for(int i=k;i<=iu;i++) 
	    this.applyTrans(j,i,i-k,C,S);
	}
	if(j < q) {
	  int jj = j-k;
	  NTrig Tg = new NTrig(S[jj].Clone(),C[jj],S[jj]);
	  this.rotg(j,j,Tg);
	  C[jj] = Tg.getCos();
	  S[jj] = Tg.getSin();
	}
      }
      for(int j=0;j<Z.numofCols();j++) 
	for(int i=k;i<q;i++) 
	  Z.applyTrans(j,i,i-k,C,S);
    }
  } 	   
  private void applyTrans(int j, int i, int ii, LNumber[] C, LNumber[] S) {
    LNumber T = (C[ii].mult(Mat[i][j])).add(S[ii].mult(Mat[i+1][j]));
    Mat[i+1][j] = (C[ii].mult(Mat[i+1][j])).sub(S[ii].mult(Mat[i][j]));
    Mat[i][j] = T;
  }
	
  /*## 
   ##  NFull.factor factors a single precision LNumbering point matrix by 
   ##  Gaussian elimination.  
   ##   
   ##  NFull.factor is usually called by NFull.condition, but it can be
   ##  called directly with a saving in time if Rcond is not needed.
   ##
   ##  On Entry:
   ##
   ##     this      NFull.Mat[n][n]
   ##               NFull.Mat is the matrix to be factored
   ##
   ##     n         int
   ##               n is the order of NFull.Mat
   ##
   ##  On Return:
   ##
   ##     this      An upper triangular matrix and multipliers used to
   ##               obtain it.  The factorization can be written
   ##               A = L*U where L is a product of permutation and 
   ##               unit lower triangular matrices and U is upper
   ##               triangular.
   ##
   ##     pivot    int[] 
   ##               pivot is a vector of pivot indices
   ##
   ##     info      int
   ##               The returned value where:
   ##               info = 0  normal value
   ##               info = k  if U(k-1,k-1) == 0, indicates solve & 
   ##                determ will divide by 0 if called
   ##*/
  
  public void factor() 
       throws SingularMatrixException
  {
    int info = 0;
    LNumber T;

    int nm = cols - 1;
    if(cols > 1) {
      
      for(int k=0;k < nm; k++) {
	int kp = k+1;
	
	/*##
	 ##  Find L the pivot index  
	 ##*/
	
	int L = this.i_amax(cols-k,1,k,k) + k;

	pivot[k] = L+1;

	/*##
	 ##  Zero pivot implies this column already triangularized  
	 ##*/

	if(!(Mat[L][k].abs()).equals(0)) {
	  /*##
	   ##  Interchange if necessary 
	   ##*/
	    
	  if(L != k) 
	    this.swapElem(L,k,k,k);

	  /*##
	   ##  Compute multipliers 
	   ##*/

	  T = (Mat[k][k].inv()).negateTo();
	  this.scal(nm-k,T,1,kp,k);

	  /*##
	   ##  Row elimination with column indexing 
	   ##*/

	  for(int j=kp;j<cols;j++) {
	    T = Mat[L][j].Clone();
	    if(L != k) 
	      this.swapElem(L,j,k,j);
	    this.axpy(nm-k,T,1,kp,k,this,1,kp,j);
	  }
	}
	else {
	  info = k+1;
	}
      }
    }
    pivot[cols-1] = cols;
    if((Mat[cols-1][cols-1].abs()).equals(0)) 
      info = cols;
    if(info != 0)
      throw new SingularMatrixException(info);
  }
  
  /*##
   ##  NFull.condition factors by Gaussian Elimination and estimates the
   ##  condition of NFull.Mat.
   ##
   ##  If Rcond is not needed, factor is slightly faster.
   ##  To solve  Mat*X = B , follow condition by solve.
   ##  To compute  Inverse(Mat)*C , follow condition by solve.
   ##  To compute  Determinant(Mat) , follow condition by determ.
   ##  To compute  Inverse(Mat) , follow condition by inverse
   ##
   ##  On Entry
   ##
   ##    Mat     NFull.Mat[n][n]
   ##            The matrix to be factored.
   ##
   ##    n       int
   ##            The order of NFull.Mat.
   ##
   ##  On Return
   ##
   ##    Mat     An upper triangular matrix and the multipliers
   ##            used to obtain it.  The factorization can be written
   ##            Mat = L*U  where  L  is a product of permutation and unit 
   ##            lower triangular matrices and  U  is upper trinagular.
   ##
   ##    pivot   int[n]
   ##            An array of pivot indices.
   ##
   ##    R       FRcond(int info, LNumber Rcond)
   ##            info:  The value returned from factor
   ##            Rcond:  An estimate of the reciprocal condition of Mat.
   ##            For the system  Mat*X = B , relative perturbations
   ##            in  Mat  and  B  of size  EPSILON  may cause
   ##            relative perturbations in  X  of size  EPSILON/Rcond .
   ##            If  Rcond  is so small that the logical expression
   ##                       1 + Rcond == 1
   ##            is true, then  Mat  may be singular to working
   ##            precision.  In particular,  Rcond  is zero if
   ##            exact singularity is detected or the estimate 
   ##            underflows.
   ##
   ##    Z       LNumber[n]
   ##            A work vector whose contents are usually unimportant.
   ##            If  Mat  is close to a singular matrix, then  Z  is 
   ##            an approximate Null vector in the sense that
   ##            norm(Mat*Z) = Rcond*norm(Mat)*norm(Z) .
   ##*/

  public LNumber condition() 
       throws SingularMatrixException
  {
    LNumber[] Z = new LNumber[cols];
    return this.condition(Z);
  }
  public LNumber condition(LNumber[] Z) 
       throws SingularMatrixException
  {
    /*##
     ##  Compute 1-Norm of A 
     ##*/
    
    LNumber anorm = this.oneNorm();
    LNumber ynorm;

    /*## 
     ##  Factor 
     ##*/

    try {
      this.factor();
    } finally {

      /*##  
	##  Rcond = 1/A.oneNorm() * est.norm(A.inverse()).  ##
	##  est.norm = norm(Z)/norm(Y) where A*Z = Y and    ##
	##  A.transpose()*Y = E.                            ##
	##  Components of E are chosen to cause max local   ##
	##  growth in W.  Arrays are frequently rescaled    ##
	##  to avoid overflow.                              ##
	##*/
      
      /*##
	##  Solve Trans(U)*W = E    
	##*/
      
      for(int j=0;j<cols;j++) {
	Z[j] = this.newOfType();
	Z[j].setZero();
      }

      this.solveTransUW(Z);
      
      LNumber S = (NUtil.asum(cols,Z,1)).invTo();
      NUtil.scal(cols,S,Z,1);

      /*##
	##  Solve Trans(L)*Y = W 
	##*/

      this.solveTransLY(Z);

      /*##
	##  Solve L*V = Y 
	##*/
      
      ynorm = this.solveLV(Z);

      /*##
	##  Solve U*Z = V 
	##*/
    
      ynorm = this.solveUZ(Z,ynorm);

      /*##
	##  Make Znorm == 1  
	##*/
    }
    LNumber R;
    
    if(!anorm.equals(0)) 
      R = ynorm.div(anorm);
    else 
      R = anorm;
    
    return R;
  }
  void solveTransUW(LNumber[] Z) {

    LNumber ek = this.newOfType();
    ek.setOne();

    for(int k=0;k<cols;k++) {
      
      if(!Z[k].equals(0)) 
	ek = NUtil.signOfA(ek,Z[k].negate());
      
      if(((ek.sub(Z[k])).abs()).greaterThan(Mat[k][k].abs())) {
	LNumber S = (Mat[k][k].abs()).div((ek.sub(Z[k])).abs());
	NUtil.scal(cols,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();
      if(!Mat[k][k].equals(0)) {
	wk.divTo(Mat[k][k]);
	wkm.divTo(Mat[k][k]);
      }      
      else {
	wk.setOne();
	wkm.setOne();
      }
      int kp = k + 1;
      if(kp < cols) {
	for(int j=kp;j<cols;j++) {
	  SM.addTo((Z[j].add(wkm.mult(Mat[k][j]))).abs());
	  Z[j] = Z[j].add(wk.mult(Mat[k][j]));
	  S.addTo(Z[j].abs());
	}
	if(SM.greaterThan(S)) {
	  LNumber T = wkm.sub(wk);
	  wk = wkm;
	  for(int j=kp;j<cols;j++) {
	    Z[j].addTo(T.mult(Mat[k][j]));
	  }
	}
      }  
      Z[k] = wk;
    }
  }
  private LNumber solveLV(LNumber[] Z) { 

    LNumber ynorm = Z[0].abs();
    ynorm.setOne();
    
    for(int k=0;k<cols;k++) {
      int L = pivot[k] - 1;
      LNumber T = Z[L].Clone();
      NUtil.swapElems(Z,L,k);
      
      if(k < cols-1) 
	this.axpy(cols-k-1,T,1,k+1,k,Z,1,k+1);
      
      if((Z[k].abs()).greaterThan(1)) {
	LNumber S = (Z[k].abs()).invTo();
	NUtil.scal(cols,S,Z,1);
	ynorm.multTo(S);
      }
    }
    LNumber S = (NUtil.asum(cols,Z,1)).invTo();
    NUtil.scal(cols,S,Z,1);
    ynorm.multTo(S);

    return ynorm;
  }
  private void solveTransLY(LNumber[] Z) {
    for(int k=cols-1;k>=0;k--) {
      if(k < cols-1) {
	Z[k].addTo(dot(cols-k-1,1,k+1,k,Z,1,k+1));
      }
      if((Z[k].abs()).greaterThan(1)) {
	LNumber S = (Z[k].abs()).invTo();
	NUtil.scal(cols,S,Z,1);
      }
      int L = pivot[k] - 1;
      NUtil.swapElems(Z,L,k);
    }
    
    LNumber S = (NUtil.asum(cols,Z,1)).invTo();
    NUtil.scal(cols,S,Z,1);
  }
  LNumber solveUZ(LNumber[] Z, LNumber ynorm) {
    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);
	ynorm.multTo(S);
      } 
      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);
    ynorm.multTo(S);

    return ynorm;
  }

  /*##
   ##  this.solve solves the single precision system
   ##  A * X = B  or  A.transpose() * X = B 
   ##  using the factors computed by this.condition or this.factor.
   ##
   ##  On Entry
   ##
   ##    this.Mat   LNumber[n][n]
   ##               The output from this.condition or this.factor.
   ##
   ##    n          int    
   ##               The order of this.Mat.      
   ##
   ##    pivot      int(n)    
   ##               The pivot vector from this.condition or this.factor.
   ##
   ##    B          LNumber[n]
   ##               The right hand side vector.
   ##
   ##    Job        int
   ##                = 0         To solve  A*X = B ,
   ##                = non-zero  To solve  A.transpose() * X = B
   ##
   ##  On Return
   ##
   ##     B         The solution vector  X.
   ##
   ##  Error condition
   ##
   ##     A division by zero will occur if the input factor contains a 
   ##     zero on the diagonal.  Technically this indicates singularity
   ##     but it is often caused by improper arguments or improper
   ##     setting of  n.  It will not occur if the subroutines are
   ##     called correctly and if this.condition has set 
   ##     LLNumber.getLNumber > 0  or this.factor threw an exception
   ##
   ##  To compute  A.inverse() * C  where  C  is a FMatrix
   ##  with  p  columns
   ##
   ##        this.condition();
   ##        if (FRcond.Rcond > threshold) {  
   ##          for(int j=0;j<p;j++) 
   ##             this.solve(n,pivot,C,j,0);
   ##        }
   ##*/

  public void solve(NFull B, int c, int Job) {
    LNumber[] F = B.getColumn(c,0);
    this.solve(F, Job);
    B.setColumn(c,F);
  }
  public void solve(LNumber[] B) {
    this.solve(B,0);
  }
  public void solve(LNumber[] B, int Job) {
    int nm = cols - 1;
    if(Job == 0) {
      /*##
       ##  Job == 0 - solve A*X = B  
       ##*/
      /*##
       ##  First solve L*Y = B     
       ##*/

      for(int k=0; k<nm;k++) {
	int L = pivot[k] - 1;
	LNumber T = B[L].Clone();
	if(L != k) {
	  B[L] = B[k];
	  B[k] = T;
	}
	this.axpy(nm-k,T,1,k+1,k,B,1,k+1);
      }

      /*##
       ##  Now solve U*X = Y  
       ##*/

      this.solveUX(B);
    }
    else {

      /*##
       ##  Job = nonzero, solve  A.transpose() * X = B  
       ##*/
      /*##
       ##  First solve  U.transpose() * Y = B
       ##*/

      this.solveTransAX(B);

      /*##
       ##  Solve  L.transpose() * X = Y
       ##*/

      for(int k=cols-2;k>=0;k--) {
	B[k].addTo(dot(nm-k,1,k+1,k,B,1,k+1));
	int L = pivot[k] - 1;
	if(L != k) 
	  NUtil.swapElems(B,L,k);
      }
    }
    return;
  }
  
  /*##
   ##  Solve U*X = Y  
   ##  This function is used by NPoFull
   ##*/
  
  void solveUX(LNumber[] B) {
    for(int k=cols-1;k>=0;k--) {
      B[k].divTo(Mat[k][k]);
      this.axpy(k,B[k].negate(),1,0,k,B,1,0);
    }
  }

  /*##
   ##  Solve  U.transpose() * Y = B
   ##  This function is used by NPoFull
   ##*/

  void solveTransAX (LNumber[] B){
    for(int k=0;k<cols;k++) {
      LNumber T = this.dot(k,1,0,k,B,1,0);
      B[k] = (B[k].sub(T)).div(Mat[k][k]);
    }
  }

  /*##
   ##  this.determ computes the determinant of a NFull using the
   ##  factors computed by this.condition or this.factor.
   ##
   ##  On Entry
   ##
   ##    this.Mat   LNumber[n][n]
   ##               The output from this.condition or this.factor.
   ##
   ##    n          int    
   ##               The order of  this.Mat.
   ##
   ##    pivot      int[n]  
   ##               The pivot vector from this.condition or this.factor.
   ##
   ##  On Return
   ##
   ##    Det        LNumber[2]
   ##               Determinant of original NFull
   ##               Determinant = Det[0] * 10.0**Det[1]
   ##               with  1  <=  Math.abs(Det[0])  <  10
   ##               or  Det[0]  ==  0.
   ##*/

  public LNumber[] determ() {

    LNumber[] Det = new LNumber[2];
    
    /*##
     ##  Compute Determinant  
     ##*/
    
    Det[0] = newOfType();
    Det[1] = newOfType();
    Det[0].setOne();
    Det[1].setZero();

    for(int i=0;i<cols;i++) {
      if(pivot[i] != i+1)
	Det[0].negateTo();
      Det[0].multTo(Mat[i][i]);
      if(Det[0].equals(0)) 
	return Det;
      else 
	NUtil.detNorm(Det);
    }
    return Det;
  }
  
  /*##
   ##  this.inverse computes the inverse of a NFull using the
   ##  factors computed by this.condition or this.factor.
   ##
   ##  On Entry
   ##
   ##    this.Mat   LNumber[n][n]
   ##               The output from this.condition or this.factor.
   ##
   ##    n          int    
   ##               The order of  this.Mat.
   ##
   ##    pivot      int[n]  
   ##               The pivot vector from this.condition or this.factor.
   ##
   ##  On Return
   ##
   ##    this.Mat   Inverse of original NFull.Mat 
   ##
   ##  Error condition
   ##
   ##     A division by zero will occur if the input factor contains
   ##     a zero on the diagonal.  It will not occur if the subroutines
   ##     are called correctly and if this.condition has set 
   ##     NRcond.Rcond > 0 or this.factor has set info == 0.
   ##*/

  public void inverse() {
    /*## 
     ##  Compute Inverse of U 
     ##*/

    this.invertU();
    
    /*##
     ##  Form inverse(U)*inverse(L)  
     ##*/

    if(cols > 1) 
      this.invUinvL();
  }
  void invertU () {
        
    for(int k=0;k<cols;k++) {
      Mat[k][k].invTo();
      LNumber T = Mat[k][k].negate();
      this.scal(k,T,1,0,k);
      for(int j=k+1;j<cols;j++) {
	T = Mat[k][j].Clone();
	Mat[k][j].setZero();
	this.axpy(k+1,T,1,0,k,this,1,0,j);
      }
    }
  }
  private void invUinvL() {
    
    LNumber[] Work = new LNumber[cols];
    
    for(int k=cols-1;k>=0;k--) {
      int kp = k+1;
      for(int i=kp;i<cols;i++) {
	Work[i] = Mat[i][k].Clone();
	Mat[i][k].setZero();
      }
      for(int j=kp;j<cols;j++) 
	this.axpy(cols,Work[j],1,0,j,this,1,0,k);

      int L = pivot[k] - 1;
      if(L != k) 
	this.swap(cols,1,0,k,this,1,0,L);
    }
  }
}