Laboratoire ICAR » Plateforme TXM

/*

 * To change this template, choose Tools | Templates

 * and open the template in the editor.

 */

package JamaPlus;

import JamaPlus.util.*;

import java.util.*;

/**

 *

 * @author Bernard

 */

public class SparseSVD {  // based on code adapted from SVDLIBC - svdLAS2 (adapted from SVDPACKC)

  private boolean transpose; // si transposition de la matrice en entrée

  private int rows, cols, dim; // nb rows, nb cols, nb de val sing. désirées

  // (après transposition éventuelle)

  private double[][] Ut; // (d x m) transpose of left singular vectors (rows of Ut)

  private double[][] Vt; // (d x n) transpose of right singular vectors (rows of Vt)

  private double[] S; // d singular values

  /* Constantes */

  private final static int MAXLL = 2; // ??

  private final static double[] END = {-1e-30, 1e-30};

  private final static double EPS = Machine.eps(); // precision de la machine

  private final static double RAD_EPS = Math.sqrt(EPS);

  private final static double RAD34_EPS = RAD_EPS*Math.sqrt(RAD_EPS);

  private final static double KAPPA = Math.max(RAD34_EPS, 1e-6);

  private double EPS1;  // = precision * sqrt(nbcols)

  /* Variables de travail : mémoire à libérer en fin de calcul */

  private double[] alph;  // diagonale de la matrice tridiagonale T

  private double[] bet; // éléments hors-diagonale de T

  private double[] ritz; // valeurs propres de la  projection de A'A sur l'espace de Krylov

  private double[] bnd; // a vérifier : initialisation (127) nécessaire ?? 

  private double[] rj ;

  private  double[] qj;

  private  double[] qj1;

  private  double[] pj ;

  private  double[] pj1;

  private  double[] wk;

   private   double[] eta; 

  private  double[] oldeta; 

//    private static int irandom = 0;  // deboguage : à supprimer

  private static class Storage {  // store Lanczos vectors q(j) 

    private double[][] P, Q;  // P supplementary storage for j < MAXLL (j= 0 or 1)

    private Storage(int nb, int maxll) {

      P = new double[maxll][];

      Q = new double[nb][];

    }

    private void storeP(int n, int j, double[] p) {

      if (P[j]==null) P[j] = new double[n];

      System.arraycopy(p, 0, P[j], 0, n);

    }

    private void storeQ(int n, int j, double[] q) {

      if (Q[j]==null) Q[j] = new double[n];

      System.arraycopy(q, 0, Q[j], 0, n);

    }

    private void retrieveP(int n, int j, double[] p) {

      System.arraycopy(P[j], 0, p, 0, n);

    }

    private void retrieveQ(int n, int j, double[] q) {

      System.arraycopy(Q[j], 0, q, 0, n);

    }

  }

  private Storage storage;

 public enum CodeSVDException {

   DIM_OUT_OF_BOUNDS ("Nombre de valeurs singulières désirées invalide"),

   FIRST_STEP_FAILURE("Impossible d'initialiser le premier pas de calcul"),

   CONVERGENCE_FAILURE("Echec de convergence ");

   String message;

  CodeSVDException(String message) {

   this.message = message;

 }

 }

  public static class SparseSVDException extends Exception { // Erreurs

    private CodeSVDException code;

    private int k;

    private SparseSVDException(CodeSVDException code, int k) {

      super();

      this.code = code;

      this.k = k;

    }

    public CodeSVDException getCode() {

      return code;

    }

    public int getK() {

      return k;

    }

    @Override

    public String getMessage() {

      return super.getMessage()+"\ncode = "+getCode().message +", n = "+getK();

    }

  }

  public SparseSVD(SparseMatrix spm, int d) throws SparseSVDException {

    int iter = Math.min(spm.getRowDimension(), spm.getColumnDimension()); // nb d'itérations 

    if (d<=0 || d>iter)

      throw new SparseSVDException(CodeSVDException.DIM_OUT_OF_BOUNDS, d);

    dim = d;

    transpose = spm.getColumnDimension()>=spm.getRowDimension()*1.2;

    SparseMatrix A = transpose ? spm.transpose() : spm.copy();

    rows = A.getRowDimension();

    cols = A.getColumnDimension();

    EPS1 = EPS*Math.sqrt(cols);

    alph = new double[iter];  // diagonale de la matrice tridiagonale T

    bet = new double[iter+1]; // éléments hors-diagonale de T

    ritz = new double[iter+1];

    bnd = new double[iter+1]; // a vérifier : initialisation (127) nécessaire ?? 

    storage = new Storage(iter, MAXLL);

    // Calcul de la matrice tridiagonale 

    int[] steps_neig_ = lanso(A, iter);

    Ut = new double[dim][rows];

    Vt = new double[dim][cols];

    S = new double[dim];

    // Calcul des  vecteurs propres

    //  nsig  : nb de val. propres "kappa-acceptables"

    int nsig = ritvec(A, steps_neig_);

    if (nsig<dim) dim = nsig;  // moins de valeurs obtenues que demandées

    // libère la mémoire

    storage = null;

    alph =  bet = ritz = bnd = null;

  }

  private int[] lanso(SparseMatrix A, int iter) throws SparseSVDException {

    // Calcul de la matrice tridiagonale

    int neig = 0;

    // espace de travail

    rj = new double[cols];

    qj = new double[cols];

     qj1 = new double[cols];

     pj = new double[cols];

     pj1 = new double[cols];

     wk = new double[cols];

    // step one 

    double[] rnm_tol_ = stpone(A, cols); // initialise rmn et tol

     eta = new double[iter]; // orthogonality estimate of Lanczos vectors

    oldeta = new double[iter];

    if (rnm_tol_[0]==0) {

      ritz[0] = alph[0];

      storage.storeQ(cols, 0, qj);

      return new int[]{0, 1}; // steps = 0 et neig = 1

    }

    eta[0] = oldeta[0] = EPS1;

    boolean[] enough_ = {false}; // flag de sortie de boucle

    int steps = 0; // j = steps

    int[] ll_ = {0};  // 0, 1 or 2 initial Lanczos vectors in local orthog.

    int first = 1;

    int last = Math.min(dim+Math.max(8, dim), iter);  // = iter dans notre version

    int i, l;

    int intro = 0;

    while (!enough_[0]) {

      steps = lanczos_step(A, first, last, 

              ll_, enough_, rnm_tol_, cols);

      steps = enough_[0] ? steps-1 : last-1;

      first = steps+1;

      bet[first] = rnm_tol_[0];

      // analyse T

      l = 0;

      for (int id2 = 0; id2<steps; id2++) {

        if (l>steps) break;

        for (i = l; i<=steps; i++) if (bet[i+1]==0) break;

        if (i>steps) i = steps;

        // now i is at the end of an unreduced submatrix

        DoubleArraysPlus.copyReverse(i-l+1, alph, l, ritz, l);

        DoubleArraysPlus.copyReverse(i-l, bet, l+1, wk, l+1);

        imtqlb(i-l+1, l, ritz, wk, bnd);

        for (int id3 = l; id3<=i; id3++)

          bnd[id3] = rnm_tol_[0]*Math.abs(bnd[id3]);

        l = i+1;

      }

      // sort eigenvalues into increasing order

      DoubleArraysPlus.sort2(ritz, bnd, steps+1, 0);

      // massage error bounds for very close ritz values 

      neig = error_bound(enough_, END[0], END[1], steps, rnm_tol_[1]);

      // should we stop? 

      if (neig<dim) {

        if (neig==0) {

          last = first+9;

          intro = first;

        } else {

          last = first+Math.max(3, 1+((steps-intro)*(dim-neig))/neig);

        }

        last = Math.min(last, iter);

      } else {

        enough_[0] = true;

      }

      enough_[0] = enough_[0]||first>=iter;

    }

    storage.storeQ(cols, steps, qj);

    // libère la mémoire

     rj = qj =  qj1 = pj =  pj1 =  wk = null;

     eta = oldeta = null;

    return new int[]{steps, neig};

  }

  private double[] stpone(SparseMatrix A,  int n)

          throws SparseSVDException {

    /*  stp-one performs the first step of the Lanczos algorithm.  It also

    does a step of extended local re-orthogonalization.

    Arguments 

    (input)

    n      dimension of the eigenproblem for matrix B

    (output)

    wptr   array of pointers that point to work space that contains

    wptr[0]             r[j]

    wptr[1]             q[j]

    wptr[2]             q[j-1]

    wptr[3]             p

    wptr[4]             p[j-1]

    alph             diagonal elements of matrix T 

    throw a SparseSVDException(-1) if the choice of initial vector (startv) failed

     */

    // get initial normalized vector; default is random 

    startv0(A, n);

    System.arraycopy(pj, 0, qj, 0, n);

    // take the first step   

    A.opb(pj, rj);

    alph[0] = DoubleArraysPlus.dot(rj, pj);

    for (int i = 0; i<n; i++) rj[i] -= alph[0]*qj[i];

    /* extended local reorthogonalization */

    double t = DoubleArraysPlus.dot(rj, pj);

    for (int i = 0; i<n; i++) rj[i] -= t*qj[i];

    alph[0] += t;

    System.arraycopy(rj, 0, pj1, 0, n);

    double rnm = Math.sqrt(DoubleArraysPlus.dot(pj1, pj1));

    double tol = RAD_EPS*(rnm+Math.abs(alph[0]));

    return new double[]{rnm, tol};

  }

//  private static double random2() {

//    double x = 0.11*(++irandom);

//    return x-Math.floor(x);

//  }

//

private void startv0(SparseMatrix A, int n)

          throws SparseSVDException {

    /* startv delivers a starting normalized vector in r, 

     * and throws a SparseSVDException(-1) if no starting vector can be found.

    Parameters 

    (input)

    n      dimension of the eigenproblem matrix B

    wptr   array of pointers that point to work space

    step      starting index for a Lanczos run

    (output)

    wptr   array of pointers that point to work space that contains

    r[j], q[j], q[j-1], p[j], p[j-1]

     */

// get initial vector; default is random

    double rnm2 = 0;

    Random random = new Random(918273);

    for (int id = 0; id<3; id++) {

      for (int i = 0; i<n; i++) wk[i]  = random.nextDouble(); 

      // apply operator to put r in range (essential if m singular)

      A.opb(wk, pj);

      rnm2 = DoubleArraysPlus.dot(pj, pj);

      if (rnm2>0) break;

    }

    // normalize

    if (rnm2<=0) throw new SparseSVDException(CodeSVDException.FIRST_STEP_FAILURE, 0);

    rnm2 = 1/Math.sqrt(rnm2);

    for (int i = 0; i<n; i++) pj[i] *= rnm2;

  }

  private int lanczos_step(SparseMatrix A, int first, int last,int[] ll_,   

          boolean[] enough_, double[] rnm_tol_, int n) throws SparseSVDException {

    /*   embodies a single Lanczos step

    Arguments 

    (input)

    n        dimension of the eigenproblem for matrix B

    first    start of index through loop                                      

    last     end of index through loop                                     

    wptr            array of pointers pointing to work space                    

    alf            array to hold diagonal of the tridiagonal matrix T

    eta      orthogonality estimate of Lanczos vectors at step j   

    oldeta   orthogonality estimate of Lanczos vectors at step j-1

    bet      array to hold off-diagonal of T                     

    ll       number of intitial Lanczos vectors in local orthog. 

    (has value of 0, 1 or 2)                        

    enough   stop flag                        

     */

    int j;

    double t;

    double mid[];

    for (j = first; j<last; j++) {

      mid = qj1;

      qj1 = qj;

      qj = mid;

      mid = pj;

      pj = pj1;

      pj1 = mid;

      storage.storeQ(n, j-1, qj1);

      if (j-1<MAXLL) storage.storeP(n, j-1, pj1);

      bet[j] = rnm_tol_[0];

      // restart if invariant subspace is found 

      if (bet[j]==0) {

        rnm_tol_[0] = startv(A, j, n);

        if (rnm_tol_[0]==0) enough_[0] = true;

      }

      if (enough_[0]) {

        /* added by Doug... */

        /* These lines fix a bug that occurs with low-rank matrices */

        mid = qj1;

        qj1 = qj;

        qj = mid;

        /* ...added by Doug */

        break;

      }

      /* take a lanczos step */

      t = 1.0/rnm_tol_[0];

      for (int i = 0; i<n; i++) qj[i] = t*rj[i];

      for (int i = 0; i<n; i++) pj[i] *= t;

      A.opb(pj, rj);

      for (int i = 0; i<n; i++) rj[i] -= rnm_tol_[0]*qj1[i];

      alph[j] = DoubleArraysPlus.dot(rj, pj);

      for (int i = 0; i<n; i++) rj[i] -= alph[j]*qj[i];

      /* orthogonalize against initial lanczos vectors */

      if (j<=MAXLL&&(Math.abs(alph[j-1])>4.0*Math.abs(alph[j])))

        ll_[0] = j;

      for (int i = 0; i<Math.min(ll_[0], j-1); i++) {

        storage.retrieveP(n, i, wk);

        t = DoubleArraysPlus.dot(wk, rj);

        storage.retrieveQ(n, i, wk);

        for (int ii = 0; ii<n; ii++) rj[ii] -= t*wk[ii];

        oldeta[i] = eta[i] = EPS1;

      }

      /* extended local reorthogonalization */

      t = DoubleArraysPlus.dot(rj, pj1);

      for (int i = 0; i<n; i++) rj[i] -= t*qj1[i];

      if (bet[j]>0) bet[j] += t;

      t = DoubleArraysPlus.dot(rj, pj);

      for (int i = 0; i<n; i++) rj[i] -= t*qj[i];

      alph[j] += t;

      System.arraycopy(rj, 0, pj1, 0, n);

      rnm_tol_[0] = Math.sqrt(DoubleArraysPlus.dot(rj, pj1));

      rnm_tol_[1] = RAD_EPS*(bet[j]+Math.abs(alph[j])+rnm_tol_[0]);

      /* update the orthogonality bounds */

      ortbnd(j, rnm_tol_[0]);

      /* restore the orthogonality state when needed */

      purge(n, ll_[0], j, rnm_tol_);

      if (rnm_tol_[0]<=rnm_tol_[1]) rnm_tol_[0] = 0;

    }

    return j;

  }

  private double startv(SparseMatrix A, int step, int n)

          throws SparseSVDException {

    /* startv delivers a starting vector in r and returns |r|; it returns 

    zero if the range is spanned, and throws a SparseSVDException(-1) if no starting 

    vector within range of operator can be found.

    Parameters 

    (input)

    n      dimension of the eigenproblem matrix B

    wptr   array of pointers that point to work space

    step      starting index for a Lanczos run

    (output)

    wptr   array of pointers that point to work space that contains

    r[j], q[j], q[j-1], p[j], p[j-1]

     */

// get initial vector; default is random

    double rnm2 = DoubleArraysPlus.dot(rj, rj);

    Random random = new Random(918273+step);

    for (int id = 0; id<3; id++) {

      for (int i = 0; i<n; i++) rj[i] = random.nextDouble();

      System.arraycopy(rj, 0, pj, 0, n);

      // apply operator to put r in range (essential if m singular)

      A.opb(pj, rj);

      System.arraycopy(rj, 0, pj, 0, n);

      rnm2 = DoubleArraysPlus.dot(rj, pj);

      if (rnm2>0) break;

    }

    if (rnm2<=0) throw new SparseSVDException(CodeSVDException.FIRST_STEP_FAILURE, 0);

    double t;

    for (int j = 0; j<step; j++) {

      storage.retrieveQ(n, j, wk);

      t = -DoubleArraysPlus.dot(pj, wk);

      for (int i = 0; i<n; i++) rj[i] += t*wk[i];

    }

    // make sure q[step] is orthogonal to q[step-1] 

    t = DoubleArraysPlus.dot(pj1, rj);

    for (int i = 0; i<n; i++) rj[i] -= t*qj1[i];

    System.arraycopy(rj, 0, pj, 0, n);

    t = DoubleArraysPlus.dot(pj, rj);

    if (t<=EPS*rnm2) rnm2 = 0;

    else rnm2 = t;

    return Math.sqrt(rnm2);

  }

  private void ortbnd(int step, double rnm) {

    /*  updates the eta recurrence

    Arguments 

    (input)

    alf      array to hold diagonal of the tridiagonal matrix T         

    eta      orthogonality estimate of Lanczos vectors at step j        

    oldeta   orthogonality estimate of Lanczos vectors at step j-1     

    bet      array to hold off-diagonal of T                          

    n        dimension of the eigenproblem for matrix B                    

    j        dimension of T                                          

    rnm            norm of the next residual vector                         

    eps1            roundoff estimate for dot product of two unit vectors

    (output)

    eta      orthogonality estimate of Lanczos vectors at step j+1     

    oldeta   orthogonality estimate of Lanczos vectors at step j        

     */

    if (step<1) return;

    if (rnm!=0) {

      if (step>1)

        oldeta[0] = (bet[1]*eta[1]+(alph[0]-alph[step])*eta[0]

                -bet[step]*oldeta[0])/rnm+EPS1;

      for (int i = 1; i<=step-2; i++)

        oldeta[i] = (bet[i+1]*eta[i+1]+(alph[i]-alph[step])*eta[i]

                +bet[i]*eta[i-1]-bet[step]*oldeta[i])/rnm+EPS1;

    }

    oldeta[step-1] = EPS1;

    double t;

    for (int i = 0; i<step; i++) {

      t = oldeta[i];

      oldeta[i] = eta[i];

      eta[i] = t;

    }

    eta[step] = EPS1;

  }

  private void purge(int n, int ll, int step, double[] rnm_tol_) {

    /*

    purge examines the state of orthogonality between the new Lanczos

    vector and the previous ones to decide whether re-orthogonalization 

    should be performed

    Arguments 

    (input)

    n        dimension of the eigenproblem for matrix B                       

    ll       number of intitial Lanczos vectors in local orthog.       

    r        residual vector to become next Lanczos vector            

    q        current Lanczos vector                                   

    ra       previous Lanczos vector

    qa       previous Lanczos vector

    wrk      temporary vector to hold the previous Lanczos vector

    eta      state of orthogonality between r and prev. Lanczos vectors 

    oldeta   state of orthogonality between q and prev. Lanczos vectors

    j        current Lanczos step                                     

    (output)

    r            residual vector orthogonalized against previous Lanczos 

    vectors

    q        current Lanczos vector orthogonalized against previous ones    

     */

    if (step<ll+2) return;

    int k = DoubleArraysPlus.indAbsMax(step-(ll+1), eta, ll);

    if (Math.abs(eta[k])<=RAD_EPS) return;

    double REPS1 = EPS1/RAD_EPS;

    double t;

    for (int iteration = 0; iteration<2; iteration++) {

      if (rnm_tol_[0]<=rnm_tol_[1]) break;

      /* bring in a lanczos vector t and orthogonalize both 

       * r and q against it */

      double tq = 0;

      double tr = 0;

      for (int i = ll; i<step; i++) {

        storage.retrieveQ(n, i, wk);

        t = -DoubleArraysPlus.dot(pj, wk);

        tq += Math.abs(t);

        for (int ii = 0; ii<n; ii++) qj[ii] += t*wk[ii];

        t = -DoubleArraysPlus.dot(pj1, wk);

        tr += Math.abs(t);

        for (int ii = 0; ii<n; ii++) rj[ii] += t*wk[ii];

      }

      System.arraycopy(qj, 0, pj, 0, n);

      t = -DoubleArraysPlus.dot(rj, pj);

      tr += Math.abs(t);

      for (int i = 0; i<n; i++) rj[i] += t*qj[i];

      System.arraycopy(rj, 0, pj1, 0, n);

      rnm_tol_[0] = Math.sqrt(DoubleArraysPlus.dot(pj1, rj));

      if (tq<=REPS1&&tr<=REPS1*rnm_tol_[0]) break;

    }

    for (int i = ll; i<=step; i++) {

      eta[i] = EPS1;

      oldeta[i] = EPS1;

    }

  }

  private static void imtqlb(int n, int i0, double[] d, double[] e, double[] bnd)

          throws SparseSVDException {

    /*   imtqlb() finds the eigenvalues of a symmetric tridiagonal

    matrix by the implicit QL method.

    Arguments 

    (input)

    n      order of the symmetric tridiagonal matrix   

    d      contains the diagonal elements of the input matrix           

    e      contains the subdiagonal elements of the input matrix in its

    last n-1 positions.  e[0] is arbitrary

    i0   starting index for the 3 vectors

    (elements with index < i0 or >= i0+n are preserved)

    (output)

    d      contains the eigenvalues in ascending order.  if an error

    exit is made, the eigenvalues are correct and ordered for

    indices 0,1,...ierr, but may not be the smallest eigenvalues.

    e      has been destroyed.                                            

    throws a SparseSVDException(j) if the j-th eigenvalue has

    not been determined after 30 iterations.                    

     */

    if (n==1) return;

    bnd[i0] = 1;

    for (int i = i0+1; i<i0+n; i++) {

      bnd[i] = 0;

      e[i-1] = e[i];

    }

    e[i0+n-1] = 0;

    for (int l = i0; l<i0+n; l++) {

      for (int iteration = 0; iteration<=30; iteration++) {

        double p = d[l];

        double f = bnd[l];

        int m;

        for (m = l; m<i0+n-1; m++) { // test de convergence

          double test = Math.abs(d[m])+Math.abs(d[m+1]);

          if (test+Math.abs(e[m])==test) break;

        }

        if (m==l) {

          d[l] = p;

          bnd[l] = f;

          /* order the eigenvalues */

          DoubleArraysPlus.sort2(d, bnd, l-i0+1, i0);

          break;

        } else {

          if (iteration==30) throw new SparseSVDException(CodeSVDException.CONVERGENCE_FAILURE, l);

          /*........ form shift ........*/

          double g = (d[l+1]-p)/(2.0*e[l]);

          double r = Math.hypot(g, 1);

          g = d[m]-p+e[l]/(g+Math.copySign(r, g));

          double s = 1;

          double c = 1;

          p = 0;

          boolean underflow = false;

          for (int i = m-1; i>=l; i--) {

            f = s*e[i];

            double b = c*e[i];

            r = Math.hypot(f, g);

            e[i+1] = r;

            if (r==0) {

              /*........ recover from underflow .........*/

              underflow = true;

              d[i+1] -= p;

              e[m] = 0.0;

              break;

            }

            s = f/r;

            c = g/r;

            g = d[i+1]-p;

            r = (d[i]-g)*s+2*c*b;

            p = s*r;

            d[i+1] = g+p;

            g = c*r-b;

            f = bnd[i+1];

            bnd[i+1] = s*bnd[i]+c*f;

            bnd[i] = c*bnd[i]-s*f;

          }

          if (!underflow) {

            d[l] -= p;

            e[l] = g;

            e[m] = 0.0;

          }

        }

      }

    }

  }

  private int error_bound(boolean[] enough_, double endl, double endr,

          int step, double tol) {

    /*

    massages error bounds for very close ritz values by placing 

    a gap between them.  The error bounds are then refined to reflect this.

    Arguments 

    (input)

    endl     left end of interval containing unwanted eigenvalues

    endr     right end of interval containing unwanted eigenvalues

    ritz     array to store the ritz values

    bnd      array to store the error bounds

    enough   stop flag

     */

    int mid = DoubleArraysPlus.indAbsMax(step+1, bnd, 1);

    for (int i = step; i>mid; i--)

      if (Math.abs(ritz[i-1]-ritz[i])<RAD34_EPS*Math.abs(ritz[i]))

        if (bnd[i]>tol&&bnd[i-1]>tol) {

          bnd[i-1] = Math.hypot(bnd[i], bnd[i-1]);

          bnd[i] = 0;

        }

    for (int i = 1; i<mid; i++)

      if (Math.abs(ritz[i+1]-ritz[i])<RAD34_EPS*Math.abs(ritz[i]))

        if (bnd[i]>tol&&bnd[i+1]>tol) {

          bnd[i+1] = Math.hypot(bnd[i], bnd[i+1]);

          bnd[i] = 0.0;

        }

    /* refine the error bounds */

    int neig = 0;

    double gapl = ritz[step]-ritz[0];

    for (int i = 0; i<=step; i++) {

      double gap = gapl;

      if (i<step) gapl = ritz[i+1]-ritz[i];

      gap = Math.min(gap, gapl);

      if (gap>bnd[i]) bnd[i] *= (bnd[i]/gap);

      if (bnd[i]<=16.0*EPS*Math.abs(ritz[i])) {

        neig++;

        if (!enough_[0]) enough_[0] = endl<ritz[i]&&ritz[i]<endr;

      }

    }

    return neig;

  }

  private int ritvec0(SparseMatrix A, int[] steps_neig_) throws SparseSVDException {

    int js, jsq, i, k, /*size,*/ id2, tmp, nsig, x;

    double[] s, xv2;

    double tmp0, tmp1, xnorm;

    double[] w1 = Vt[0];

    double[] w2 = new double[cols];

    js = steps_neig_[0]+1;

    jsq = js*js;

    /*size = sizeof(double) * n;*/

    s = new double[jsq];

    xv2 = new double[cols];

    /* initialize s to an identity matrix */

    for (i = 0; i<jsq; i += (js+1)) s[i] = 1.0;

    DoubleArraysPlus.copyReverse(js, alph, 0, w1, 0);

    DoubleArraysPlus.copyReverse(steps_neig_[0], bet, 1, w2, 1);

    /* on return from imtql2(), w1 contains eigenvalues in ascending 

     * order and s contains the corresponding eigenvectors */

    double[][] z = new double[js][js];

    for (int r = 0; r<js; r++) for (int c = 0; c<js; c++) z[r][c] = s[r*js+c];

    imtql2(js, js, w1, w2, z);

    for (int r = 0; r<js; r++) for (int c = 0; c<js; c++) s[r*js+c] = z[r][c];

    /*fwrite((char *)&n, sizeof(n), 1, fp_out2);

    fwrite((char *)&js, sizeof(js), 1, fp_out2);

    fwrite((char *)&kappa, sizeof(kappa), 1, fp_out2);*/

    /*id = 0;*/

    nsig = 0;

    x = 0;

    id2 = jsq-js;

    for (k = 0; k<js; k++) {

      tmp = id2;

      if (bnd[k]<=KAPPA*Math.abs(ritz[k])&&k>js-steps_neig_[1]-1) {

        if (--x<0) x = dim-1;

        w1 = Vt[x];

        for (i = 0; i<cols; i++) w1[i] = 0.0;

        for (i = 0; i<js; i++) {

          storage.retrieveQ(cols, i, w2);

          for (int ii = 0; ii<cols; ii++) w1[ii] += s[tmp]*w2[ii];

          tmp -= js;

        }

        /*fwrite((char *)w1, size, 1, fp_out2);*/

        /* store the w1 vector row-wise in array xv1;   

         * size of xv1 is (steps+1) * (nrow+ncol) elements 

         * and each vector, even though only ncol long,

         * will have (nrow+ncol) elements in xv1.      

         * It is as if xv1 is a 2-d array (steps+1) by     

         * (nrow+ncol) and each vector occupies a row  */

        /* j is the index in the R arrays, which are sorted by high to low 

        singular values. */

        /*for (i = 0; i < n; i++) R->Vt->value[x]xv1[id++] = w1[i];*/

        /*id += nrow;*/

        nsig++;

      }

      id2++;

    }

    s = null;

    /* Rotate the singular vectors and values. */

    /* x is now the location of the highest singular value. */

    DoubleArraysPlus.rotateArray(Vt, x);

    dim = Math.min(dim, nsig);

    for (x = 0; x<dim; x++) {

      /* multiply by matrix B first */

      A.opb(Vt[x], xv2);

      tmp0 = DoubleArraysPlus.dot(Vt[x], xv2);

      for (int ii = 0; ii<cols; ii++) xv2[ii] -= tmp0*Vt[x][ii];

      tmp0 = Math.sqrt(tmp0);

      xnorm = Math.sqrt(DoubleArraysPlus.dot(xv2, xv2));

      /* multiply by matrix A to get (scaled) left s-vector */

      A.opa(Vt[x], Ut[x]);

      tmp1 = 1.0/tmp0;

      for (int ii = 0; ii<rows; ii++) Ut[x][ii] *= tmp1;

      xnorm *= tmp1;

      bnd[i] = xnorm;

      S[x] = tmp0;

    }

    xv2 = w2 = null;

    return nsig;

    /* multiply by matrix A to get (scaled) left s-vector */

  }

  private int ritvec(SparseMatrix A, int[] steps_neig_) throws SparseSVDException {

    //  computes the singular vectors of matrix A               *

/*  Parameters

    (input)

    nrow       number of rows of A

    steps      number of Lanczos iterations performed

    fp_out2    pointer to unformatted output file

    n              dimension of matrix A

    kappa      relative accuracy of ritz values acceptable as 

    eigenvalues of A'A

    ritz       array of ritz values

    bnd        array of error bounds

    alf        array of diagonal elements of the tridiagonal matrix T

    bet        array of off-diagonal elements of T

    w1, w2     work space

    (output)

    xv1        array of eigenvectors of A'A (right singular vectors of A)

    ierr              error code

    0 for normal return from imtql2()

    k if convergence did not occur for k-th eigenvalue in

    imtql2()

    nsig       number of accepted ritz values based on kappa

    (local)

    s              work array which is initialized to the identity matrix

    of order (j + 1) upon calling imtql2().  After the call,

    s contains the orthonormal eigenvectors of the symmetric 

    tridiagonal matrix T

     */

    // int n = cols;

    double[] w1 = Vt[0];

    double[] w2 = new double[cols];

    int js = steps_neig_[0]+1;

    double[][] s = new double[js][js];

    double[] xv2 = new double[cols];

    /* initialize s to an identity matrix */

    for (int i = 0; i<js; i++) s[i][i] = 1;

    DoubleArraysPlus.copyReverse(js, alph, 0, w1, 0);

    DoubleArraysPlus.copyReverse(steps_neig_[0], bet, 1, w2, 1);

    /* on return from imtql2(), w1 contains eigenvalues in ascending 

     * order and s contains the corresponding eigenvectors */

    imtql2(js, js, w1, w2, s);

    int nsig = 0;

    int x = 0;

    for (int k = 0; k<js; k++) {

      if (bnd[k]<=KAPPA*Math.abs(ritz[k])&&k>js-steps_neig_[1]-1) {

        if (--x<0) x = dim-1;

        w1 = Vt[x];

        for (int i = 0; i<cols; i++) w1[i] = 0;

        for (int i = 0; i<js; i++) {

          storage.retrieveQ(cols, i, w2);

          for (int ii = 0; ii<cols; ii++) w1[ii] += s[k][js-1-i]*w2[ii];

        }

        /*fwrite((char *)w1, size, 1, fp_out2);*/

        /* store the w1 vector row-wise in array xv1;   

         * size of xv1 is (steps+1) * (nrow+ncol) elements 

         * and each vector, even though only ncol long,

         * will have (nrow+ncol) elements in xv1.      

         * It is as if xv1 is a 2-d array (steps+1) by     

         * (nrow+ncol) and each vector occupies a row  */

        /* j is the index in the R arrays, which are sorted by high to low 

        singular values. */

        /*for (i = 0; i < n; i++) R->Vt->value[x]xv1[id++] = w1[i];*/

        /*id += nrow;*/

        nsig++;

      }

    }

    s = null;  // libère la mémoire

  /* Rotate the singular vectors and values. */

    /* x is now the location of the highest singular value. */

    DoubleArraysPlus.rotateArray(Vt, x);

    int d = Math.min(dim, nsig);

    for (x = 0; x<d; x++) {

      /* multiply by matrix B first */

      A.opb(Vt[x], xv2);

      double tmp0 = DoubleArraysPlus.dot(Vt[x], xv2);

      for (int i = 0; i<cols; i++) xv2[i] -= tmp0*Vt[x][i];

      tmp0 = Math.sqrt(tmp0);

      double xnorm = Math.sqrt(DoubleArraysPlus.dot(xv2, xv2));

      /* multiply by matrix A to get (scaled) left s-vector */

      A.opa(Vt[x], Ut[x]);

      double tmp1 = 1/tmp0;

      for (int i = 0; i<rows; i++) Ut[x][i] *= tmp1;

      xnorm *= tmp1;

      //   bnd[i] = xnorm; bizarre ??

      S[x] = tmp0;

    }

    xv2 = w2 = null;

    return nsig;

  }

  private static void imtql2(int nm, int n, double d[], double[] e, double[][] z)

          throws SparseSVDException {

    /*   imtqlb() finds the eigenvalues  and eigenvectors of a symmetric tridiagonal

    matrix by the implicit QL method.

    Arguments 

    (input)

    nm     row dimension of the symmetric tridiagonal matrix   

    n      order of the symmetric tridiagonal matrix   

    d      contains the diagonal elements of the input matrix           

    e      contains the subdiagonal elements of the input matrix in its

    last n-1 positions.  e[0] is arbitrary

    z      contains the identity matrix

    (output)

    d      contains the eigenvalues in ascending order.  if an error

    exit is made, the eigenvalues are correct but unordered for

    indices 0,1,...ierr.

    e      has been destroyed.                                            

    z      contains orthonormal eigenvectors of the symmetric   

    tridiagonal (or full) matrix.  if an error exit is made,

    z contains the eigenvectors associated with the stored 

    eigenvalues.                                        

    throws a SparseSVDException(j) if the j-th eigenvalue has

    not been determined after 30 iterations.                    

     */

    if (n==1) return;

    for (int i = 1; i<n; i++) e[i-1] = e[i];

    e[n-1] = 0.0;

    int nnm = n*nm;

    for (int l = 0; l<n; l++) {

      /* look for small sub-diagonal element */

      for (int iteration = 0; iteration<=30; iteration++) {

        double p = d[l];

        int m;

        for (m = l; m<n-1; m++) {// test de convergence

          double test = Math.abs(d[m])+Math.abs(d[m+1]);

          if (test+Math.abs(e[m])==test) break;

        }

        if (m==l) break;

        /* set error -- no convergence to an eigenvalue after

         * 30 iterations. */

        if (iteration==30) throw new SparseSVDException(CodeSVDException.CONVERGENCE_FAILURE, l);

        /* form shift */

        double g = (d[l+1]-p)/(2*e[l]);

        double r = Math.hypot(g, 1);

        g = d[m]-p+e[l]/(g+Math.copySign(r, g));

        double s = 1;

        double c = 1;

        p = 0;

        boolean underflow = false;

        for (int i = m-1; i>=l; i--) {

          double f = s*e[i];

          double b = c*e[i];

          r = Math.hypot(f, g);

          e[i+1] = r;

          if (r==0) {

            /*........ recover from underflow .........*/

            underflow = true;

            d[i+1] -= p;

            e[m] = 0.0;

            break;

          }

          s = f/r;

          c = g/r;

          g = d[i+1]-p;

          r = (d[i]-g)*s+2*c*b;

          p = s*r;

          d[i+1] = g+p;

          g = c*r-b;

          /* form vector */

          for (int k = 0; k<nm; k++) {

            f = z[i+1][k];

            z[i+1][k] = s*z[i][k]+c*f;

            z[i][k] = c*z[i][k]-s*f;

          }

        }

        if (!underflow) {

          d[l] -= p;

          e[l] = g;

          e[m] = 0.0;

        }

      }

    }

    /* order the eigenvalues and corresponding eigenvectors */

    DoubleArraysPlus.sort2(d, z, n, 0);

  }

  /* ------------------------

  Public Methods

   * ------------------------ */

  /** Return the left singular vectors

  @return     U

   */

  public Matrix getU() {

    return transpose ? (new Matrix(Vt, dim, cols)).transpose()

            : (new Matrix(Ut, dim, rows)).transpose();

  }

  /** Return the right singular vectors

  @return     V

   */

  public Matrix getV() {

    return transpose ? (new Matrix(Ut, dim, rows)).transpose()

            : (new Matrix(Vt, dim, cols)).transpose();

  }

  /** Return the one-dimensional array of singular values

  @return     diagonal of S.

   */

  public double[] getSingularValues() {

    return Arrays.copyOfRange(S, 0, dim);

  }

}