s_sis2.cc

矩阵奇异分解(svd)最新c++版本

// File: s_sis2.cc -- Implements sis2 class
// Author: Suvrit Sra
// Date: 15, Nov, 2003

/*************************************************************************
 THE ORIGINAL SVDPAKC COPYRIGHT
                           (c) Copyright 1993
                        University of Tennessee
                          All Rights Reserved                          
 *************************************************************************/

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cerrno>
#include <cstring>
#include <unistd.h>
#include <fcntl.h>


#include "s_sis2.h"

using namespace ssvd;
/****************************************************************
 *                                                              *
 *                        ritzit()                              *
 *                                                              *
 ****************************************************************/
/****************************************************************

 Description:
 ------------

 This subroutine is a translation of the Fortran-77 procedure
 RITZIT from the SVDPACK library written by Michael W. Berry,
 University of Tennessee, Dept. of Computer Science, 107 Ayres
 Hall, Knoxville TN 37919-1301

 This subroutine determines the absolutely largest eigenvalues
 and corresponding eigenvectors of a real symmetric matrix by
 simultaneous iteration.

 External parameters (constants)
 -------------------------------

 Defined and documented in sisg.h (sis1c.h)

 local parameters: 
 -----------------
 
 (input)

 n      the order of the matrix  whose eigenpairs are sought
        (matrix B for the SVD problem)
 kp     the number of simultaneous iteration vectors
 km     the maximum number of iteration steps to be performed. If
        starting values for the iteration vectors are available,
        km should be prefixed with a minus sign.
 eps    the tolerance for accepting eigenvectors
 opb    the name of the subroutine that defines the matrix B. opb
        is called with the parameters (n, u, w) and must 
        compute w=Bu without altering u for the SVD problem
 inf    the name of the subroutine that may be used to obtain 
        information or exert control during execution. inf is
        called with parameters (ks, kg, kh, f, m), where ks
        is the number of the next iteration step, kg is the
        number of already accepted eigenvectors, kh is the number
        already accepted eigenvalues, and f is the array of error
        quantities for the vectors of x. An element of f has the
        value 4.0 until the corresponding eigenvalue of the matrix
        B  has been accepted.
 kem    the number of eigenvalues and corresponding eigenvectors
         of matrix B desired. kem must be less than kp.
 x      contains, if km is negative, the starting values for the
        iteration vectors.

 (output)
 
 km     is unchanged
 kem    is reset to the number of eigenvalues and eigenvectors
        of matrix B actually accepted within the limit of km steps
 imem   the number of bytes needed for this invocation
 x      contains in its first kem columns orthonormal
        eigenvectors of the matrix B,  corresponding to the 
        eigenvalues in array d.  The remaining columns contain 
        approximations to further eigenvectors.
 d      contains in its first kem positions the absolutely
        largest eigenvalues of the matrix B. The remaining positions
        contain approximations to smaller eigenvalues.

 u, w, b, f, cx, x, s and  rq are temporary storage arrays.

 Functions used
 --------------

 BLAS:  ddot, dscal, daxpy
 EISP:  tred2, tql2, pythag
 MISC:  opb, inf, dmax, imax
 
*****************************************************************/
void s_sis2::ritzit(FILE* fptr, long n, long kp, long km, double eps,long kem, 
					double **x, double *d, double *f, 
					double *cx, double *u, long *imem)

{
  long i, j, l, i1, l1, size, kg, kh, kz, kz1, kz2, ks, m, m1;
  long  jj, ii, ig, ip, ik, jp, k, flag1 ;  
  double  *w[NMAX], *rq, *b[NSIG], *s, *tptr, TmpRes;
  double xks, ee2, e, e2, e1, xkm, xk, xm1, t;

  /*******************************************************************
   * allocate memory for temporary storage arrays used in ritzit only*
   *                                                                 *
   *            rq - (kp)                                            *
   *             w - (NSIG * n)                                      *
   *             b - (NSIG * kp)                                     *
   *             s - (kp)                                            *
   *******************************************************************/

   size =    kp * (NSIG +2) + NSIG * n;

   rq = (double*) malloc (size*sizeof(double));
   if (rq == 0) {
	 perror("MALLOC FAILED IN RITZIT()");
	 exit(errno);
   }

   tptr= rq + kp;
   for (i=0; i< NSIG; i++){
    w[i]= tptr;
    tptr+=n;
   }

   /*w=tptr; tptr+=(n * kp);*/
   for (i=0; i<NSIG; i++) {
    b[i]=tptr;
    tptr+=kp;
   };
   s=tptr;

   /*finished allocating memory*/


  *imem += size;
  *imem = sizeof(double) * (*imem);

  ee2 = ONE + 1.0e-1 * eps;
  e   = ZERO;
  kg  = -1;
  kh  = -1;
  kz  = 1367;
  kz1 = 0;
  kz2 = 0;
  ks  = 0;
  m   = 1;

  for (l=0; l<kp; l++){
    f[l] = 4.0e+0;
   cx[l] = ZERO;
   rq[l] = ZERO;
 } 
 
  if (km >= 0)    /*generate random initial iteration vectors*/
      for (j=0; j<kp; j++)
       for (l=0; l<n; l++){
         kz = (3125 * kz) % 65536;  
        x[j][l] = (double) (kz - 32768);
      }
  km =abs(km);
  l1 = 1;
  i1 = 1;
  jp=kp;
flag1=0;
#define NRow kp
#define NCol n

  /* extend orthonormalization to all kp rows of x 

     Variables used here:
 
     NCol                     : number of columns of x
     NRow                     : number of rows of x
     ii, ik, jp, k, l1        : indexing integers
     t, TmpRes                : double precision storage.


     at the end of this section of code, 

     x   : transpose(Q)
     b   : R

     where X = QR is the orthogonal factorization of the matrix X
     stored in array x

     invokes: ddot, daxpy, dscal (BLAS)
              sqrt               (math.h)
  */

    ik = l1 -1; 
  for (i=0; i< jp; i++){
    for (k=l1-1; k<jp; k++) b[i][k]= ZERO;
   
    if (i >= l1-1){
      ik=i+1;
      b[i][i]=sqrt(ddot(NCol, &x[i][0], 1 ,  &x[i][0], 1));
      t=ZERO;
      if (b[i][i] != ZERO) t=ONE / b[i][i];
      dscal(NCol, t, &(x[i][0]),  1);
     }

   for (ii=ik; ii< NRow; ii++){
     TmpRes=ZERO;
     for (jj=0; jj<NCol; jj++)
       TmpRes+= x[ii][jj] * x[i][jj];
     b[i][ii]=TmpRes;
  }
  
    for (k=ik; k<jp;k++) 
        daxpy(NCol, -b[i][k], &x[i][0], 1, &x[k][0], 1);

  } /* end for i */

  ig = 0;
  ip = kp - 1;
  
 /*statement 70  from original RITZIT begins here*/
while (1){

  flag1=1;  /* so that jacobi step is done at least once */

  while (flag1){
        /* jacobi step modified */
        for(k=ig; k< kp; k++){
          mopb(n, &x[k][0], &w[0][0]);
          for (j=0; j<n; j++) x[k][j]=w[0][j];
         }
        l1=ig + 1;
        jp=kp;
        flag1=0; /*flag1 is set only if re-orthog needs to be done*/

     /* extend orthonormalization to all kp rows of x 

     Variables used here:
 
     NCol                     : number of columns of x
     NRow                     : number of rows of x
     ii, ik, jp, k, l1        : indexing integers
     t, TmpRes                : double precision storage.


     at the end of this section of code, 

     x   : transpose(Q)
     b   : R

     where X = QR is the orthogonal factorization of the matrix X
     stored in array x

     invokes: ddot, daxpy, dscal (BLAS)
              sqrt               (math.h)
  */

    ik = l1 -1; 
  for (i=0; i< jp; i++){
    for (k=l1-1; k<jp; k++) b[i][k]= ZERO;
   
    if (i >= l1-1){
      ik=i+1;
      b[i][i]=sqrt(ddot(NCol, &x[i][0], 1 ,  &x[i][0], 1));
      t=ZERO;
      if (b[i][i] != ZERO) t=ONE / b[i][i];
      dscal(NCol, t, &(x[i][0]),  1);
     }

   for (ii=ik; ii< NRow; ii++){
     TmpRes=ZERO;
     for (jj=0; jj<NCol; jj++)
       TmpRes+= x[ii][jj] * x[i][jj];
     b[i][ii]=TmpRes;
  }
  
    for (k=ik; k<jp;k++) 
        daxpy(NCol, -b[i][k], &x[i][0], 1, &x[k][0], 1);

  } /* end for i */

       if ( ks<=0 ) {
         /* measures against unhappy choice of initial vectors*/
         for (k=0;k<kp; k++)
           if (b[k][k] ==ZERO)
              for (j=0;j<n;j++){
                 kz= (3125 * kz) % 65536;
                 x[k][j] = (double)(kz-32768);
                 flag1=1;
               }
       }
     if (flag1){
       l1=1;
       ks=1; /*we dont want to re-initialize x[][] again*/

   /* extend orthonormalization to all kp rows of x 

     Variables used here:
 
     NCol                     : number of columns of x
     NRow                     : number of rows of x
     ii, ik, jp, k, l1        : indexing integers
     t, TmpRes                : double precision storage.


     at the end of this section of code, 

     x   : transpose(Q)
     b   : R

     where X = QR is the orthogonal factorization of the matrix X
     stored in array x

     invokes: ddot, daxpy, dscal (BLAS)
              sqrt               (math.h)
  */

    ik = l1 -1; 
  for (i=0; i< jp; i++){
    for (k=l1-1; k<jp; k++) b[i][k]= ZERO;
   
    if (i >= l1-1){
      ik=i+1;
      b[i][i]=sqrt(ddot(NCol, &x[i][0], 1 ,  &x[i][0], 1));
      t=ZERO;
      if (b[i][i] != ZERO) t=ONE / b[i][i];
      dscal(NCol, t, &(x[i][0]),  1);
     }

   for (ii=ik; ii< NRow; ii++){
     TmpRes=ZERO;
     for (jj=0; jj<NCol; jj++)
       TmpRes+= x[ii][jj] * x[i][jj];
     b[i][ii]=TmpRes;
  }
  
    for (k=ik; k<jp;k++) 
        daxpy(NCol, -b[i][k], &x[i][0], 1, &x[k][0], 1);

  } /* end for i */
      } 
    } /*end while flag1 */

  for (k= ig; k<kp; k++)
    for (l=k; l<kp; l++){
      t = ZERO;

      for (i=l; i<kp; i++)
        t+= b[k][i] * b[l][i];

      /* negate matrix to reverse eigenvalue ordering */
      b[k][l] = -t;
    }
  j=kp - kg - 1;

  tred2(ig, kp, b, d, u, b);
  ii=tql2(ig, kp, d, u, b);


  for (k=ig; k< kp; k++)
   d[k]=sqrt(dmax(-d[k], ZERO));
  

 for (j=ig; j<kp; j++)
  for (k=0; k<n; k++) 
   w[j][k]=ZERO;

   for (j=ig; j<kp; j++)
     for (l=ig; l<kp; l++){
       TmpRes=b[j][l];
       for (k=0; k<n; k++)
         w[j][k] += TmpRes  * x[l][k];}   /* w is going to be transposed as 
                                         compared to the fortran version */
     

  for (j=ig; j<kp; j++)
   for (k=0; k<n; k++)
      x[j][k]=w[j][k]; 

  xks=(double)(++ks);
  if (d[kp-1] > e) e=d[kp-1 ];

  /* randomization */
  if (kz1<3) {
    for (j=0; j<n; j++) {
       kz= (3125 * kz) % 65536;
       x[kp-1][j] = (double) (kz -32768);
    }
    l1=kp;
    jp=kp;

  /* extend orthonormalization to all kp rows of x 

     Variables used here:
 
     NCol                     : number of columns of x
     NRow                     : number of rows of x
     ii, ik, jp, k, l1        : indexing integers
     t, TmpRes                : double precision storage.


     at the end of this section of code, 

     x   : transpose(Q)
     b   : R

     where X = QR is the orthogonal factorization of the matrix X
     stored in array x

     invokes: ddot, daxpy, dscal (BLAS)
              sqrt               (math.h)
  */

    ik = l1 -1; 
  for (i=0; i< jp; i++){
    for (k=l1-1; k<jp; k++) b[i][k]= ZERO;
   
    if (i >= l1-1){
      ik=i+1;
      b[i][i]=sqrt(ddot(NCol, &x[i][0], 1 ,  &x[i][0], 1));
      t=ZERO;
      if (b[i][i] != ZERO) t=ONE / b[i][i];
      dscal(NCol, t, &(x[i][0]),  1);
     }

   for (ii=ik; ii< NRow; ii++){
     TmpRes=ZERO;
     for (jj=0; jj<NCol; jj++)
       TmpRes+= x[ii][jj] * x[i][jj];
     b[i][ii]=TmpRes;
  }
  
    for (k=ik; k<jp;k++) 
        daxpy(NCol, -b[i][k], &x[i][0], 1, &x[k][0], 1);

  } /* end for i */

  }
  
  /*compute control quantities cx */
  for (k=ig; k<ip;k++){
    t=(d[k] - e)*(d[k] + e);
    if (t <= ZERO) cx[k]=ZERO; 
    else 
           if (e ==ZERO) cx[k] = 1.0e+3 + log(d[k]);
           else  cx[k]= log( (d[k]+sqrt(t)) / e);
  }
    
    /*acceptance test for eigenvalues including adjustment of 
      kem and kh such that