/*## 
 ##  class FFull
 ##
 ##  This class extends FMatrix to represent all full matrices.   
 ##  It contains the subroutine for singular value decomposition, the 
 ##  QR decompostion routines, the updating routines, and the Gaussian 
 ##  Elimination factoring routines.
 ##
 ##*/

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

public class FFull extends FMatrix {

  public FFull() {
    super();
    cons();
  }
  public FFull(int i, int j) {
    cons(i,j);
  }
  public FFull(float[][] f) {
    cons(f);
  }
  public FFull(FFull F) {
    cons(F);
  }
  public FFull(int ni, int nj, float[][] f, int strow, int stcol) {
    rows = ni;
    cols = nj;
    Mat = new float[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];
  }

  public FBanded band(int ml, int mu) {
    FBanded B = new FBanded(ml,mu,cols);
    int M = ml+mu;
    for(int j=0;j<cols;j++) {
      int i1 = Math.max(0,j-mu);
      int i2 = Math.min(cols-1,j+ml);
      for(int i=i1;i<=i2;i++) {
	int k = i-j+M;
	B.Mat[k][j] = this.Mat[i][j];
      }
    }
    return B;
  }
  public void setPivot(int[] f) {
    for(int i=0;i<cols;i++) 
      pivot[i] = f[i];
  }
  public void setPivot(int i, int p) {
    pivot[i] = p;
  }

  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 {
      float[][] G = new float[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;
    }
  }
  public void transpose(FFull N) {
    for(int i=0;i<rows;i++) {
      for(int j=0;j<cols;j++) {
	N.Mat[j][i] = Mat[i][j];
      }
    }
  }

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

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

  private float oneNorm() {
    float a = 0;
    for(int j=0;j<cols;j++) 
      a = Math.max(a, this.asum(cols,1,0,j));
    return a;
  }
  
  /*##  
   ##  svDecompose is a subroutine to reduce a single precision nXp 
   ##  FFull by orthogonal transformations U and V to diagonal form.  
   ##  The diagonal elements S[i] are the singular values of this 
   ##  FFull.  The columns of U are the corresponding left singular
   ##  Vectors, and the columns of V the right singular vectors.
   ##
   ##  On Entry:
   ##
   ##     this      FFull where this.numofRows() >= n.
   ##               Contains the FFull whose singular value
   ##               decomposition is to be computed.  This FFull
   ##               is destroyed by svDecompose
   ##
   ##     n         int
   ##               n is the number of columns of the section of this
   ##               FFull that is to be used.
   ##
   ##     p         int
   ##               p is the number of rows of the section of this 
   ##               FFull that is to be used.
   ##
   ##     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(n,p) 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         float[MM], where MM = Math.min(n+1,p).    
   ##               the first min(n,p) entries of S contain the 
   ##               singular values of this FFull arranged in 
   ##               descending order of magnitude.
   ##
   ##     E         float[p]
   ##               E ordinarily contains zeros.  However see the 
   ##               discussion of info for exceptions.
   ##
   ##     U         float[n][k], where if job(A) == 1 then k == n
   ##                               else k == Min(n,p)
   ##               U contains the FFull of right singular vectors.
   ##               U is not referenced if job(A) == 0.  
   ##
   ##     V         float[p][p]
   ##               V contains the FFull 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(n,p)).  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(float[] S, float[] E, FFull U, FFull V, int job) {
    
    float[] W = new float[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;
    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++) {
      int qp = q+1;
      if(q < NCT) { 
	S[q] = this.nrm2(rows-q,1,q,q);
	if(S[q] != 0) {
	  if(this.Mat[q][q] != 0) 
	    S[q] = FUtil.signOf(S[q],this.Mat[q][q]);
	  this.scal(rows-q,1/S[q],1,q,q);
	  this.Mat[q][q]+=1;
	}
	S[q] = -S[q];
      }
      for(int j=qp;j<cols;j++) {
	if(q < NCT) 
	  if(S[q] != 0) {
	    float T = -this.dot(rows-q,1,q,q,this,1,q,j)/this.Mat[q][q];
	    this.axpy(rows-q,T,1,q,q,this,1,q,j);
	  }
	E[j] = this.Mat[q][j];
      }
      if(wantU && q<NCT) 
	for(int i=q;i<rows;i++) 
	  U.Mat[i][q] = this.Mat[i][q];
      
      /*##  
       ##  Compute the q-th row transformation and place the q-th
       ##   super-diagonal in E[q]  
       ##*/

      if(q < NRT) {
	E[q] = FUtil.nrm2(cols-qp,E,1,qp);
	if(E[q] != 0) {
	  if(E[qp] != 0) 
	    E[q] = FUtil.signOf(E[q],E[qp]);
	  FUtil.scal(cols-qp,1/E[q],E,1,qp);
	  E[qp]+=1;
	}
	E[q] = -E[q];
	if(qp < rows && E[q] != 0) {
	  for(int i=qp;i<rows;i++) 
	    W[i] = 0;
	  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]/E[qp],W,1,qp,1,qp,j);
	}
	if(wantV) 
	  for(int i=qp;i<cols;i++) 
	    V.Mat[i][q] = E[i];
      }
    }
    int M = Math.min(rows+1,cols);
    if(NCT < cols) 
      S[NCT] = this.Mat[NCT][NCT];
    if(rows < M) 
      S[M-1] = 0;
    if(NRT+1 < M) 
      E[NRT] = this.Mat[NRT][M-1];
    E[M-1] = 0;

    if(wantU) 
      U.generateU(NCT,NCU,S);
    if(wantV) 
      V.generateV(NRT,E);

    int iter = 0;
    int MM = M;
    while(M>0) {
      if(iter >= MAXIT) 
	return M;
      else {
	int q = qgetter(M-2,S,E);
	if(q == M-2) {
	  q++;
	  if(S[q] < 0) {
	    S[q] = -S[q];
	    if(wantV) 
	      V.scal(cols,-1,1,0,q);
	  }
	  boolean notDone = true;
	  while(q < MM-1 && notDone) {
	    if(S[q] < S[q+1]) {
	      FUtil.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 {
	  int qs = this.findNextStep(q,M,S,E);
	  if(qs == q) {
	    q++;
	    FTrig Tg = this.fggetter(q,M,S,E);

	    for(int k=q;k<M-1;k++) {
	      FUtil.rotg(Tg);
	      if(k != q) 
		E[k-1] = Tg.getA();
	      Tg.setA(Tg.getCos()*S[k] + Tg.getSin()*E[k]);
	      E[k] = Tg.getCos()*E[k] - Tg.getSin()*S[k];
	      Tg.setB(Tg.getSin()*S[k+1]);
	      S[k+1] *= Tg.getCos();
	      if(wantV) 
		V.rot(cols,1,0,k,V,1,0,k+1,Tg);
	      FUtil.rotg(Tg);
	      S[k] = Tg.getA();
	      Tg.setA(Tg.getCos()*E[k] + Tg.getSin()*S[k+1]);
	      S[k+1] = -Tg.getSin()*E[k] + Tg.getCos()*S[k+1];
	      Tg.setB(Tg.getSin()*E[k+1]);
	      E[k+1] *= 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) {
	    q++;
	    int MMOne = M-2;
	    FTrig Tg = new FTrig();
	    Tg.setB(E[MMOne]);
	    E[MMOne] = 0;
	    for(int k=MMOne;k>=q;k--) {
	      Tg.setA(S[k]);
	      FUtil.rotg(Tg);
	      S[k] = Tg.getA();
	      if(k != q) {
		Tg.setB(-Tg.getSin()*E[k-1]);
		E[k-1] *= Tg.getCos();
	      }
	      if(wantV) 
		V.rot(cols,1,0,k,V,1,0,MMOne+1,Tg);
	    }
	  }
	  else {
	    q = qs;
	    FTrig Tg = new FTrig();
	    Tg.setB(E[q]);
	    E[q] = 0;
	    q++;
	    for(int k=q;k<M;k++) {
	      Tg.setA(S[k]);
	      FUtil.rotg(Tg);
	      S[k] = Tg.getA();
	      Tg.setB(-Tg.getSin()*E[k]);
	      E[k] *= 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, float[] S) {
    for(int j=NCT;j<NCU;j++) {
      for(int i=0;i<cols;i++) {
	this.Mat[i][j] = 0;
      }
      this.Mat[j][j] = 1;
    }
    if(NCT > 0) {
      for(int q=NCT-1;q>=0;q--) {
	if(S[q] != 0) {
	  for(int j=q+1;j<NCU;j++) {
	    float T = -this.dot(cols-q,1,q,q,this,1,q,j)/this.Mat[q][q];
	    this.axpy(cols-q,T,1,q,q,this,1,q,j);
	  }
	  this.scal(cols-q,-1,1,q,q);
	  this.Mat[q][q]++;
	  for(int i=0;i<q;i++) 
	    this.Mat[i][q] = 0;
	}
	else {
	  for(int i=0;i<cols;i++) 
	    this.Mat[i][q] = 0;
	  this.Mat[q][q] = 1;
	}
      }
    }
  }
  private void generateV(int NRT, float[] E) {
    for(int q=cols-1;q>=0;q--) {
      int qp = q+1;
      if((q < NRT) && (E[q] != 0)) {
	for(int j=qp;j<rows;j++) {
	  float T = -this.dot(rows-qp,1,qp,q,this,1,qp,j)/this.Mat[qp][q];
	  this.axpy(rows-qp,T,1,qp,q,this,1,qp,j);
	}	    
      }
      for(int i=0;i<rows;i++) 
	this.Mat[i][q] = 0;
      this.Mat[q][q] = 1;
    }
  }
  private int qgetter(int st, float[] S, float[] E) {
    for(int q=st;q >= 0;q--) {
      float Test = Math.abs(S[q]) + Math.abs(S[q+1]);
      float ZTest = Test + Math.abs(E[q]);
      if(ZTest == Test) {
	E[q] = 0;
	return q;
      }
    }
    return -1;
  }
  private int findNextStep(int q, int M, float[] S, float[] E) {
    for(int qs=M-1;qs>q;qs--) {
      float Test = 0;
      if(qs != M-1) 
	Test += Math.abs(E[qs]);
      if(qs != q+1) 
	Test += Math.abs(E[qs-1]);
      float Ztest = Test + Math.abs(S[qs]);
      if(Ztest == Test) {
	S[qs] = 0;
	return qs;
      }
    }
    return q;
  }
  private FTrig fggetter(int q, int M, float[] S, float[] E) {
    FTrig Tg = new FTrig();
    float scale = FUtil.aMax(S[M-1],S[M-2],E[M-2],S[q],E[q]);
    float SM = S[M-1]/scale;
    float SMM = S[M-2]/scale;
    float EMM = E[M-2]/scale;
    float Sq = S[q]/scale;
    float Eq = E[q]/scale;
    float B = ((SMM + SM)*(SMM - SM) + (float)Math.pow(EMM,2))/2;
    float C = (float)Math.pow(SM*EMM,2);
    float shift = 0;
    if(B != 0 || C != 0) {
      shift = (float)Math.sqrt(Math.pow(B,2) + C);
      if(B < 0) 
	shift = -shift;
      shift = C/(B+shift);
    }
    Tg.setA((Sq + SM)*(Sq - SM) - shift);
    Tg.setB(Sq*Eq);
  
    return Tg;
  }

  /*##
   ##  qrDecompose uses Householder transformations to compute the
   ##  QR factorization of an n by p FFull.  Column pivoting
   ##  based on the 2-norms of the reduced columns may be
   ##  performed at the user's option.
   ##
   ##  On Entry:
   ##
   ##    this    FFull - Mat[L][p], where L >= n            
   ##            FFull contains the FFull whose decomposition is to 
   ##            be computed.
   ##
   ##
   ##    n       int 
   ##            n is the number of rows of FFull.Mat to be used in
   ##            the calculation
   ##
   ##    p       int 
   ##            p is the number of columns of FFull.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 FFull.Mat[0][k]
   ##            is placed in one of three classes according to the 
   ##            value of pivot[k]:
   ##
   ##               if pivot[k] > 0, then FFull.getColumn(k,0) is an
   ##                                  initial column.
   ##
   ##               if pivot[k] == 0, then FFull.getColumn(k,0) is a
   ##                                  free column.
   ##
   ##               if pivot[k] < 0, then FFull.getColumn(k,0) is a 
   ##                                  final column.
   ##
   ##            Before the decomposition is computed, initial columns
   ##            are moved to the beginning of FFull.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 FFull.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    FFull.Mat contains in its upper triangle the upper
   ##            triangular matrix R of the QR factorization.
   ##            Below its diagonal FFull.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 FFull.Mat, but that of FFull.Mat
   ##            with its columns permuted as described by pivot.
   ##
   ##    QRaux   float[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(float[] 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);
	  FUtil.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(float[] QRaux) {
    Pul P = new Pul(0,-1);
    this.qrDecom(P,QRaux);
  }
  private void qrDecom(Pul P, float[] QRaux) {
    float[] Work = new float[cols];

    for(int j=P.pl;j<=P.pu;j++) {
      QRaux[j] = this.nrm2(rows,1,0,j);
      Work[j] = QRaux[j];
    }
    int qup = Math.min(rows,cols);
    for(int q=0;q<qup;q++) {
      if(q >= P.pl && q<P.pu) {
	float maxnorm = 0;
	int maxj = q;
	for(int j=q;j<=P.pu;j++) {
	  if(QRaux[j] >= maxnorm) {
	    maxnorm = QRaux[j];
	    maxj = j;
	  }
	}
	if(maxj != q) {
	  this.swap(rows,1,0,q,this,1,0,maxj);
	  QRaux[maxj] = QRaux[q];
	  Work[maxj] = Work[q];
	  FUtil.swapElems(pivot,maxj,q);
	}
      }
      QRaux[q] = 0;
      if(q < rows-1) {
	float normxq = this.nrm2(rows-q,1,q,q);
	if(normxq != 0) {
	  if(this.Mat[q][q] != 0) 
	    normxq = FUtil.signOf(normxq,this.Mat[q][q]);
	  this.scal(rows-q,1/normxq,1,q,q);
	  this.Mat[q][q]+=1;
	  for(int j=q+1;j<cols;j++) {
	    float T = -this.dot(rows-q,1,q,q,this,1,q,j)/this.Mat[q][q];
	    this.axpy(rows-q,T,1,q,q,this,1,q,j);
	    if(j >= P.pl && j <= P.pu) {
	      if(QRaux[j] != 0) {
		float TT = 1 - (float)Math.pow(Math.abs(this.Mat[q][j])/
					       QRaux[j],2);
		TT = Math.max(TT,0);
		T = TT;
		TT = 1 + (float).5*TT*(float)Math.pow(QRaux[j]/Work[j],2);
		if(TT != 1) 
		  QRaux[j] *= (float)Math.sqrt(T);
		else {
		  QRaux[j] = this.nrm2(rows-q-1,1,q+1,j);
		  Work[j] = QRaux[j];
		}
	      }
	    }
	  }
	  QRaux[q] = this.Mat[q][q];
	  this.Mat[q][q] = -normxq;
	}
      }
    }
  }

  /*##
   ##  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 FFull.Mat that was input to qrDecompose (if no
   ##  pivoting was done, XK consists of the first k+1 columns of 
   ##  FFull.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 FFull.Mat and QRaux
   ##
   ##  On Entry
   ##
   ##    this   FFull.Mat[n][p]
   ##           FFull.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  float[p]
   ##           QRaux contains the auxiliary output from qrDecompose.
   ##
   ##    Y      float[n]
   ##           Y contains the vector that is to be manipulated
   ##           by qr***.
   ##
   ##  On Return from its respective function
   ##
   ##  qrQy:
   ##    Qy     float[n]
   ##           Qy contains Q*Y
   ##
   ##  qrQty:
   ##    Qty    float[n]
   ##           Qty contains Q.transpose()*Y
   ##
   ##  qrQB:
   ##    B      float[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 FFull.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    float[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     float[n]
   ##           Xb contains the least squares approximation XK*B,
   ##           Xb is also the orthogonal projection of Y onto the 
   ##           column space of FFull.Mat.
   ##*/
  
  public void qrQy(float[] QRaux, float[] Y, float[] Qy) {
    
    int k = Math.min(rows,cols);
    int ju = Math.min(k,rows-1) - 1;
    
    if(ju < 0) 
      Qy[0] = Y[0];
    else {
      FUtil.copy(rows,Y,1,0,Qy,1,0);
      for(int j=ju;j>=0;j--) 
	this.qrsolver(j,Qy,QRaux);
    }
  }
  
  public void qrQty(float[] QRaux, float[] Y,float[] Qty) {
    
    int k = Math.min(rows,cols);
    int ju = Math.min(k,rows-1) - 1;
    
    if(ju < 0) 
      Qty[0] = Y[0];
    else {
      FUtil.copy(rows,Y,1,0,Qty,1,0);
      for(int j=0;j<=ju;j++) 
	this.qrsolver(j,Qty,QRaux);
    }
  }

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

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

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

  public void qrXB(float[] QRaux, float[] Y,
		  float[] Qty, float[] 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];
    else {
      FUtil.copy(k,Qty,1,0,Xb,1,0);
      for(int i=k;i<rows;i++) 
	Xb[i] = 0;
      for(int j=ju;j>=0;j--) 
	this.qrsolver(j,Xb,QRaux);
    }
  }
  public void qrRSD(float[] QRaux, float[] Y,
		    float[] Qty, float[] Rsd) {

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

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

    for(int j=0;j<cols;j++) {
      float Xj = X[j];
      
      for(int i=0;i<j;i++) {
	float T = C[i]*this.Mat[i][j] + S[i]*Xj;
	Xj = C[i]*Xj - S[i]*this.Mat[i][j];
	this.Mat[i][j] = T;
      }
      FTrig Tg = new FTrig(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++) {
      float Zeta = Y[j];
      for(int i=0;i<cols;i++) {
	float T = C[i]*Z[i][j] + S[i]*Zeta;
	Zeta = C[i]*Zeta - S[i]*Z[i][j];
	Z[i][j] = T;
      }
      float azeta = Math.abs(Zeta);
      if(azeta != 0 && Rho[j] >= 0) {
	float scale = azeta*Rho[j];
	Rho[j] = scale*(float)Math.sqrt(Math.pow(azeta/scale,2)+
					Math.pow(Rho[j]/scale,2));
      }
    }
  }
  public int chDowndate(float[] X,float[][] Z,
		      float[] Y,float[] Rho,float[] C,float[] S) {
    
    int info = 0;
    float Zeta;

    S[0] = X[0]/this.Mat[0][0];
    for(int j=1;j<cols;j++) {
      S[j] = X[j] - this.dot(j,1,0,j,S,1,0);
      S[j] = S[j]/this.Mat[j][j];
    }
    float norm = FUtil.nrm2(cols,S,1,0);
    if(norm >= 1) 
      info = -1;
    else {
      float alpha = (float)Math.sqrt(1 - (float)Math.pow(norm,2));
      for(int i=cols-1;i>=0;i--) {
	float scale = alpha + Math.abs(S[i]);
	float A = alpha/scale;
	float B = S[i]/scale;
	norm = (float)Math.sqrt(Math.pow(A,2) + Math.pow(B,2));
	C[i] = A/norm;
	S[i] = B/norm;
	alpha = scale*norm;
      }
      for(int j=0;j<cols;j++) {
	float XX = 0;
	for(int i=j;i>=0;i--) {
	  float T = C[i]*XX + S[i]*this.Mat[i][j];
	  this.Mat[i][j] = C[i]*this.Mat[i][j] - S[i]*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] - S[i]*Zeta)/C[i];
	  Zeta = C[i]*Zeta - S[i]*Z[i][j];
	}
	float azeta = Math.abs(Zeta);
	if(azeta > Rho[j]) {
	  info = 1;
	  Rho[j] = -1;
	}
	else 
	  Rho[j] *= (float)Math.sqrt(1-(float)Math.pow(azeta/Rho[j],2));
      }
    }
    return info;
  }

  public void chExchange(int k, int q, FFull Z,
		       float[] C, float[] S, int job) {
    
    /* k and l 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++) {
	int ii = q-i;
	S[i] = this.Mat[ii][q];
      }
      for(int j=qm;j>=k;j--) {
	for(int i=0;i<=j;i++) 
	  this.Mat[i][j+1] = this.Mat[i][j];
	this.Mat[j+1][j+1] = 0;
      }
      for(int i=0;i<k;i++) {
	int ii = q-i;
	this.Mat[i][k] = S[ii];
      }
      float T = S[0];
      FTrig Tg = new FTrig();
      for(int i=0;i<qmkp;i++) {
	Tg.setA(S[i+1]);
	Tg.setB(T);
	FUtil.rotg(Tg);
	S[i+1] = Tg.getA();
	T = S[i+1];
	C[i] = Tg.getCos();
	S[i] = Tg.getSin();
      }
      this.Mat[k][k] = T;
      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.cols;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++) {
	int ii = qmkp + i;
	S[ii] = this.Mat[i][k];
      }
      for(int j=k;j<q;j++) {
	for(int i=0;i<=j;i++) 
	  this.Mat[i][j] = this.Mat[i][j+1];
	S[j-k] = this.Mat[j+1][j+1];
      }
      for(int i=0;i<=k;i++) 
	this.Mat[i][q] = S[qmkp+i];

      for(int i=kp;i<=q;i++) 
	this.Mat[i][q] = 0;

      for(int j=k;j<cols;j++) {
	if(j != k) {
	  int iu = Math.min(j-1,qm);
	  for(int i=k;i<=iu;i++) {
	    int ii = i-k;
	    this.applyTrans(j,i,ii,C,S);
	  }
	}
	if(j < q) {
	  int jj = j-k;
	  FTrig Tg = new FTrig(S[jj],C[jj],S[jj]);
	  this.rotg(j,j,Tg);
	  C[jj] = Tg.getCos();
	  S[jj] = Tg.getSin();
	}
      }
      for(int j=0;j<Z.cols;j++) 
	for(int i=k;i<q;i++) {
	  int ii = i-k;
	  Z.applyTrans(j,i,ii,C,S);
	}
    }
  } 	   
  private void applyTrans(int j, int i, int ii, float[] C, float[] S) {
    float T = C[ii]*this.Mat[i][j] + S[ii]*this.Mat[i+1][j];
    this.Mat[i+1][j] = C[ii]*this.Mat[i+1][j] - S[ii]*this.Mat[i][j];
    this.Mat[i][j] = T;
  }
	
  /*## 
   ##  FFull.factor factors a single precision floating point matrix by 
   ##  Gaussian elimination.  
   ##   
   ##  FFull.factor is usually called by FFull.condition, but it can be
   ##  called directly with a saving in time if Rcond is not needed.
   ##
   ##  On Entry:
   ##
   ##     this      FFull.Mat[n][n]
   ##               FFull.Mat is the matrix to be factored
   ##
   ##     n         int
   ##               n is the order of FFull.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;
    float 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(this.Mat[L][k] != 0) {
	  /*##
	   ##  Interchange if necessary 
	   ##*/
	    
	  if(L != k) 
	    this.swapElem(L,k,k,k);

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

	  T = (float)-1.0/this.Mat[k][k];
	  this.scal(nm-k,T,1,kp,k);

	  /*##
	   ##  Row elimination with column indexing 
	   ##*/
	  
	  for(int j=kp;j<cols;j++) {
	    T = this.Mat[L][j];
	    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(this.Mat[cols-1][cols-1] == 0) 
      info = cols;
    if(info != 0)
      throw new SingularMatrixException(info);
  }
  
  /*##
   ##  FFull.condition factors by Gaussian Elimination and estimates the
   ##  condition of FFull.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     FFull.Mat[n][n]
   ##            The matrix to be factored.
   ##
   ##    n       int
   ##            The order of FFull.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, float 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       float[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 WrongDataTypeException, SingularMatrixException
  {
    Vector Z = new FVector(cols);
    return this.condition(Z);
  }
  public LNumber condition(Vector Ze) 
       throws WrongDataTypeException, SingularMatrixException
  {
    
    float[] Z = Ze.getFloatArray();
    
    /*##
     ##  Compute 1-Norm of A 
     ##*/
    
    float anorm = this.oneNorm();
    float 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.  Vectors are frequently rescaled   ##
	##  to avoid overflow.                              ##
	##*/
      
      /*##
	##  Solve Trans(U)*W = E    
	##*/
      
      for(int j=0;j<cols;j++) 
	Z[j] = 0;
      this.solveTransUW(Z);
      
      float S = (float)1.0/FUtil.asum(cols,Z,1);
      FUtil.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  
	##*/
    }
    LFloat R;
    
    if(anorm != 0) 
      R = new LFloat(ynorm/anorm);
    else
      R = new LFloat(anorm);

    return R;
  }
  void solveTransUW(float[] Z) {

    float ek = 1;

    for(int k=0;k<cols;k++) {
      
      if(Z[k] != 0) 
	ek = FUtil.signOf(ek,-Z[k]);
      
      if(Math.abs(ek-Z[k]) > Math.abs(this.Mat[k][k])) {
	float S = Math.abs(this.Mat[k][k])/Math.abs(ek-Z[k]);
	FUtil.scal(cols,S,Z,1);
	ek *= S;
      }
      float wk = ek - Z[k];
      float wkm = -ek - Z[k];
      float S = Math.abs(wk);
      float SM = Math.abs(wkm);
      if(this.Mat[k][k] != 0) {
	wk /= this.Mat[k][k];
	wkm /= this.Mat[k][k];
      }      
      else {
	wk = 1;
	wkm = 1;
      }
      int kp = k + 1;
      if(kp < cols) {
	for(int j=kp;j<cols;j++) {
	  SM += Math.abs(Z[j] + wkm*this.Mat[k][j]);
	  Z[j] += wk*this.Mat[k][j];
	  S += Math.abs(Z[j]);
	}
	if(S < SM) {
	  float T = wkm - wk;
	  wk = wkm;
	  for(int j=kp;j<cols;j++) {
	    Z[j] += T*this.Mat[k][j];
	  }
	}
      }  
      Z[k] = wk;
    }
  }
  private float solveLV(float[] Z) { 

    float ynorm = 1;
    
    for(int k=0;k<cols;k++) {
      int L = pivot[k] - 1;
      float T = Z[L];
      FUtil.swapElems(Z,L,k);
      
      if(k < cols-1) 
	this.axpy(cols-k-1,T,1,k+1,k,Z,1,k+1);
      
      if(Math.abs(Z[k]) > 1) {
	float S = (float)1.0/Math.abs(Z[k]);
	FUtil.scal(cols,S,Z,1);
	ynorm *= S;
      }
    }
    float S = (float)1.0/FUtil.asum(cols,Z,1);
    FUtil.scal(cols,S,Z,1);
    ynorm *= S;

    return ynorm;
  }
  private void solveTransLY(float[] Z) {
    for(int k=cols-1;k>=0;k--) {
      if(k < cols-1) {
	Z[k] += this.dot(cols-k-1,1,k+1,k,Z,1,k+1);
      }
      if(Math.abs(Z[k]) > 1) {
	float S = (float)1.0/Math.abs(Z[k]);
	FUtil.scal(cols,S,Z,1);
      }
      int L = pivot[k] - 1;
      FUtil.swapElems(Z,L,k);
    }
    
    float S = (float)1.0/FUtil.asum(cols,Z,1);
    FUtil.scal(cols,S,Z,1);
  }
  float solveUZ(float[] Z, float ynorm) {
    for(int k=cols-1;k>=0;k--) {
      if(Math.abs(Z[k]) > Math.abs(this.Mat[k][k])) {
	float S = Math.abs(this.Mat[k][k])/Math.abs(Z[k]);
	FUtil.scal(cols,S,Z,1);
	ynorm *= S;
      } 
      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);
    }
    float S = (float)1.0/FUtil.asum(cols,Z,1);
    FUtil.scal(cols,S,Z,1);
    ynorm *= 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   float[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          float[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 
   ##     LFloat.getFloat > 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(FFull B, int c, int Job) 
       throws WrongDataTypeException
  {
    Vector F = new FVector(B.getColumn(c,0));
    this.solve(F, Job);
    B.setColumn(c,F.getFloatArray());
  }
  public void solve(Vector B) 
       throws WrongDataTypeException
  {
    this.solve(B,0);
  }
  public void solve(Vector Be, int Job) 
       throws WrongDataTypeException
  {
    float[] B = Be.getFloatArray();

    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;
	float T = B[L];
	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] += this.dot(nm-k,1,k+1,k,B,1,k+1);
	int L = pivot[k] - 1;
	if(L != k) 
	  FUtil.swapElems(B,L,k);
      }
    }
    return;
  }
  
  /*##
   ##  Solve U*X = Y  
   ##  This function is used by FPoFull
   ##*/
  
  void solveUX(float[] B) {
    for(int k=cols-1;k>=0;k--) {
      B[k] /= this.Mat[k][k];
      this.axpy(k,-B[k],1,0,k,B,1,0);
    }
  }

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

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

  /*##
   ##  this.determ computes the determinant of a FFull using the
   ##  factors computed by this.condition or this.factor.
   ##
   ##  On Entry
   ##
   ##    this.Mat   float[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        float[2]
   ##               Determinant of original FFull
   ##               Determinant = Det[0] * 10.0**Det[1]
   ##               with  1  <=  Math.abs(Det[0])  <  10
   ##               or  Det[0]  ==  0.
   ##*/

  public Vector determ() {

    float[] Det = new float[2];
    
    /*##
     ##  Compute Determinant  
     ##*/
    
    Det[0] = 1;
    Det[1] = 0;

    for(int i=0;i<cols;i++) {
      if(pivot[i] != i+1)
	Det[0] = -Det[0];
      Det[0] *= this.Mat[i][i];
      if(Det[0] == 0) {
	FVector V = new FVector(Det);
	return V;
      }      
      else 
	FUtil.detNorm(Det);
    }
    FVector V2 = new FVector(Det);
    return V2;
  }
  
  /*##
   ##  this.inverse computes the inverse of a FFull using the
   ##  factors computed by this.condition or this.factor.
   ##
   ##  On Entry
   ##
   ##    this.Mat   float[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 FFull.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 
   ##     FRcond.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++) {
      this.Mat[k][k] = (float)1.0/this.Mat[k][k];
      float T = -this.Mat[k][k];
      this.scal(k,T,1,0,k);
      for(int j=k+1;j<cols;j++) {
	T = this.Mat[k][j];
	this.Mat[k][j] = 0;
	this.axpy(k+1,T,1,0,k,this,1,0,j);
      }
    }
  }
  private void invUinvL() {
    
    float[] Work = new float[cols];
    
    for(int k=cols-1;k>=0;k--) {
      int kp = k+1;
      for(int i=kp;i<cols;i++) {
	Work[i] = this.Mat[i][k];
	this.Mat[i][k] = 0;
      }
      for(int j=kp;j<cols;j++) {
	float T = Work[j];
	this.axpy(cols,T,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);
    }
  }
}