svdpack.h

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

/* svdpack.h -- interface to SVDPACK   -*- c++ -*-
   Copyright (C) 2003 Suvrit Sra ([e-mail]) */


// (c) Suvrit Sra 2004 All Rights Reserved


#ifndef _SVDPACK_H
#define _SVDPACK_H

#include <cstdio>
#include <mat/SparseMatrix.h>

namespace ssvd {
#define DBG(x) std::cerr << x;

/* Class svdpack -- a crude first attempt to provide for myself a simple
 * interface to the various methods in the netlib SVDPACK.
 * Right now the interface is crude with very little stuff factored
 * out...basically i just want to get this thing running so don't have time
 * to make it nice etc....if you want it to be niced out, please let me
 * know, or contribute....thank you.
 * Date: 13 Nov, 2003
 * Author: Suvrit Sra ([e-mail])
 */

class svdpack {
protected:
  /* Shared constant values among all derived classes */
  static const int MINBLKS  =  2;
  static const int TRANSP   =  1;
  static const int NTRANSP  =  0;
  static const int TRUE     =  1;
  static const int FALSE    =  0;
  static const double CONST =  100.0;
  static const double ZERO  =  0.0;
  static const double ONE   =  1.0;
  static const double RDWARF=  3.834e-20;
  static const double RGIANT=  1.304e19;

  /* Standard params restricting input size etc.*/
  long int NZMAX;
  long int NMAX;
  long int SEED;
  long int NCMAX;
  long int LMTNW;

  /* Common variables for all classes      */
  long    *pointr ,      /* pointer to column start array         */
          *rowind ;      /* pointer to row indices array          */
  double  *value  ;      /* pointer to nonzero values array       */

  double  *sing;	 /* singular values are filled into this */
  long    mxvcount,      /* matrix-vector multiplications counter */ 
          mtxvcount,     /* transposed matrix-vector mult. counter*/
          iconv;         /* converged vector counter              */


  long	  ierr,		 /* error flag			       */
          j,		 /* number of lanczos steps taken      */
          neig,		 /* number of ritz values stabilized   */
          nsig,		 /* number of accepted ritz values     *
			  * based on kappa (relative accuracy) */
          ncol,		 /* number of columns of A	       */
          nrow,		 /* number of rows of A		       */
          m_nz;          /* number of non zeros in A           */
  
  double *alpha,         /* diagonal elements of bidiagonal matrix*
			  * from single vector Lanczos (inner)    *
			  * recursion                             */
         *beta,          /* super-diagonals of bidiagonal matrix  *
			  * from single vector Lanczos (inner)    *
			  * recursion                             */
         *p,             /* work array                            */
         *q,             /* work array                            */
         *t,             /* work array                            */
         *z;             /* work array                            */
  double *tres,          /* temporary residual array              */
         *y,             /* work array                            */
         *temp,          /* temporary work space                  */
         *uu,            /* work array                            */
         *vv,            /* work array                            */
         *u0,            /* array of converged left S-vectors     */
         *v0,            /* array of converged right S-vectors    */
         *uvtmp,         /* temporary work space                  */
         *pp,            /* left S-vectors of block upper bi-     *
			  * diagonal matrix from outer recursion  */ 
         *qq;            /* right S-vectors of block upper bi-    *
			  * diagonal matrix from outer recursion  */

  double **uup,          /* corresponding 2-dimensional           */
         **yp,           /* ...array representation of the above  */
         **vvp,          /* ...linear arrays                      */
         **uvtmpp, 
         **ppp, 
         **qqp,
         **ztempp;

  /* From bls1.h */
  double *xv1,           /* temp arrays needed for computing   */
         *xv2,           /* singular vectors		       */
         *a;	         /* pointer to area used by user-      *
			  * supplied procedure store and holds *
			  * lanczos vectors 		       */

  /* Copied from las1.h */
  double  rnm,		 /* norm of the next residual vector   */
          anorm,
	  tol,
	  eps,		 /* positive machine epsilon	      */
	  eps1,		 /* roundoff estimate for dot product  *
			  * of two unit vector		      */
	  reps,
	  eps34;

  double  *ztemp;

  long iter, nn;

  /* Common set of methods used across all derived classes */
  
  /* BLAS routines, replace them by using <clapack.h> */
  double  dasum  (long, double *, long                       );  
  void    datx   (long, double, double *,long, double *, long);
  void    daxpy  (long, double, double *,long, double *, long);

  void    dcopy  (long n,double *dx,long incx,double *dy,
		  long incy                                  );
  double  ddot   (long, double *,long, double *, long        );

  void    dgemm  (long, long, long, long, long, double, 
		  double **, double *, double, double *      );
  void    dgemm2 (long, long, long, long, long, double, 
		  double **, double **, double, double **    );
  void    dgemm3 (long, long, long, long, long, double,
		  double **, long, long, double **, long,
		  long, double, double **, long, long        );

  void    dgemv  (long, long, long, double, double **,
		  double *, double, double *                 );
  void    dgemv2 (long, long, long, double, double **, long,
		  long, double *, double, double *           );

  double  dmax   (double, double                             );
  double  dmin   (double, double                             );

  void    drot   (long, double *, double *, double, double   );
  void    drotg  (double *, double *, double *, double *     );

  void    dscal  (long, double, double *,long                );
  void    dsbmv  (long, long, double, double **, 
		  double *, double, double *                 );
  void    dsort2 (long,long, double *, double *              );
  void    dswap  (long, double *, long, double *, long       );
  void    dtbmv  (long, long, long, double **, double *      );

  double  enorm  (long, double *                             );
  double  fsign  (double , double                            );


  long    imin   (long, long                                 );
  long    imax   (long, long                                 );
  void    imtql2 (long, long, double *, double *, double *   );

  double  mrandom(long *                                     );
  double  norm_1 ();
  void    opa    (double *, double *                         );
  void    opa    (long, long, double *, double *             );
  void    opat   (double *,double *                          );
  void    opat   (long, double *, double *                   );

  void    opb    (long, double *, double *                   );
  void    opb    (long, double *, double *, bool             );
  void    opb    (long, long, double *, double *             );
  //void    opb    (long, double *, double *, double           );
  void    opm    (long, long, long, double **, double **     );

  void    orthg  (long, long, long, double **, double **, 
		  double *                                   );  
  double  pythag (double, double                             );
  void    qriter2(long, double *, double *, double **,
		  double **                                  );

  float   timer  (void);
  long    tql2   (long, long, double *, double *, double **  );
  long    tql2   (long, double*, double*, double**           );
  void    tred2  (long, long, double **, double *,
		  double *, double **                        );

  int dcomp(double* a, double* b)
  {
    if (*a < *b) return -1;
    if (*a > * b) return 1;
    return 0;
  }

  void write_header(long, FILE*, bool, size_t, size_t, int, char*);
public:
  svdpack()
  { mxvcount = mtxvcount = ierr = 0; NMAX = 20000; NCMAX = 20000; NZMAX = 10000000; SEED = 91827211; LMTNW = 10000000;}
  
  svdpack(long, long, long);
  
  virtual int runIt() = 0;

  int    init       (char*, char*);	// init(fname, txx);
  void   extern_init(SparseMatrix* s);
  void   extern_init(long*, long*, double*, long, long, long);
  void   cleanup() { /* deallocate stuff etc. */ }	

  double* getSing() { return sing;}
  virtual ~svdpack() {}
};

}
#endif