s_las2.cc

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

// File: s_las2.cc -- Implements the las2 class
// Author: Suvrit Sra
// Date: 14 nov, 2003

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

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

#include "s_las2.h"


#define STORQ 1
#define RETRQ 2
#define STORP 3
#define RETRP 4
#define MAXLL 2

using namespace ssvd;
/***********************************************************************
 *                                                                     *
 *				landr()				       *
 *        Lanczos algorithm with selective orthogonalization           *
 *                    Using Simon's Recurrence                         *
 *                       (double precision)                            *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   landr() is the LAS2 driver routine that, upon entry,
     (1)  checks for the validity of input parameters of the 
	  B-eigenproblem 
     (2)  determines several machine constants
     (3)  makes a Lanczos run
     (4)  calculates B-eigenvectors (singular vectors of A) if requested 
	  by user


   arguments
   ---------

   (input)
   n        dimension of the eigenproblem for A'A
   lanmax   upper limit of desired number of Lanczos steps
   maxprs   upper limit of desired number of eigenpairs
   nnzero   number of nonzeros in matrix A
   endl     left end of interval containing unwanted eigenvalues of B
   endr     right end of interval containing unwanted eigenvalues of B
   vectors  1 indicates both eigenvalues and eigenvectors are wanted
              and they can be found in output file lav2; 
	    0 indicates only eigenvalues are wanted
   kappa    relative accuracy of ritz values acceptable as eigenvalues
	      of B (singular values of A)
   r        work array

   (output)
   j        number of Lanczos steps actually taken                     
   neig     number of ritz values stabilized                           
   ritz     array to hold the ritz values                              
   bnd      array to hold the error bounds


   External parameters
   -------------------

   Defined and documented in las2.h


   local parameters
   -------------------

   ibeta    radix for the floating-point representation
   it       number of base ibeta digits in the floating-point significand
   irnd     floating-point addition rounded or chopped
   machep   machine relative precision or round-off error
   negeps   largest negative integer
   wptr	    array of pointers each pointing to a work space


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

   MISC         dmax, machar, check_parameters
   LAS2         ritvec, lanso

 ***********************************************************************/

long s_las2::landr(long fpo2, FILE* fp, bool ascii, long n, 
		   long lanmax, long maxprs, long nnzero, double endl,  
		   double endr, bool vectors, double kappa, double *ritz, 
		   double *bnd, double *r)

{
   long i, size, ibeta, it, irnd, machep, negep;
   double *wptr[10], *tptr, *tptr2;

   /* data validation */
   if (check_parameters(maxprs, lanmax, n, endl, endr, vectors, nnzero))
      return(-1);

   /* Compute machine precision */ 
   machar(&ibeta, &it, &irnd, &machep, &negep);

   eps1 = eps * sqrt( (double) n );
   reps = sqrt(eps);
   eps34 = reps * sqrt(reps);

   /* allocate work area and initialize pointers         *
    * ptr             symbolic name         size         *
    * wptr[0]             r                  n           *
    * wptr[1]             q		     n           *
    * wptr[2]             q_previous         n           *
    * wptr[3]             p		     n           *
    * wptr[4]             p_previous         n           *
    * wptr[5]             wrk                n           *
    * wptr[6]             alf              lanmax        *
    * wptr[7]             eta              lanmax        *
    * wptr[8]             oldeta           lanmax        *
    * wptr[9]             bet              lanmax+1      */
   //fprintf(stderr, "s_las2::landr() about to start doing computation\n");
   size = 5 * n + (lanmax * 4 + 1);
   tptr = NULL;
   if (!(tptr = (double *) malloc(size * sizeof(double)))){
      perror("FIRST MALLOC FAILED in LANDR()");
      exit(errno);
   }
   tptr2 = tptr;
   wptr[0] = r;
   for (i = 1; i <= 5; i++) {
      wptr[i] = tptr;
      tptr += n;
   }
   for (i = 6; i <= 9; i++) {
      wptr[i] = tptr;
      tptr += lanmax;
   }

   lanso(n, lanmax, maxprs, endl, endr, ritz, bnd, wptr);

   /* compute eigenvectors */
   if (vectors) {
      if (!(xv1 = (double *) malloc((nrow+ncol)*(j+1)*sizeof(double))) ||
	  !(xv2 = (double *) malloc(ncol * sizeof(double)))) {
	  perror("SECOND MALLOC FAILED in LANDR()");
	  exit(errno);
	 }
      kappa = dmax(fabs(kappa), eps34);
      ritvec(fpo2, fp, ascii, n, kappa, ritz, bnd, wptr[6], wptr[9], wptr[4], wptr[5]);
   }

   free(tptr2);
   return(0);
}





/***********************************************************************
 *                                                                     *
 *                         stpone()                                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   Function 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)
   ierr   error flag
   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]
	    wptr[6]             diagonal elements of matrix T 


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

   BLAS		daxpy, datx, dcopy, ddot, dscal
   USER		store, opb
   LAS		startv

 ***********************************************************************/

void s_las2::stpone(long n, double *wrkptr[])

{
   double t, *alf;
   alf = wrkptr[6];

   /* get initial vector; default is random */
   rnm = startv(n, wrkptr);
   if (rnm == 0.0 || ierr != 0) return;

   /* normalize starting vector */
   t = 1.0 / rnm;
   datx(n, t, wrkptr[0], 1, wrkptr[1], 1);
   dscal(n, t, wrkptr[3], 1);

   /* take the first step */
   mopb(n, wrkptr[3], wrkptr[0]);
   alf[0] = ddot(n, wrkptr[0], 1, wrkptr[3], 1);
   daxpy(n, -alf[0], wrkptr[1], 1, wrkptr[0], 1);
   t = ddot(n, wrkptr[0], 1, wrkptr[3], 1);
   daxpy(n, -t, wrkptr[1], 1, wrkptr[0], 1);
   alf[0] += t;
   dcopy(n, wrkptr[0], 1, wrkptr[4], 1);
   rnm = sqrt(ddot(n, wrkptr[0], 1, wrkptr[4], 1));
   anorm = rnm + fabs(alf[0]);
   tol = reps * anorm;
   return;
}


/***********************************************************************
 *                                                                     *
 *                         startv()                                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   Function delivers a starting vector in r and returns |r|; it returns 
   zero if the range is spanned, and ierr is non-zero 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
   j      starting index for a Lanczos run
   eps    machine epsilon (relative precision)

   (output)
   wptr   array of pointers that point to work space that contains
	  r[j], q[j], q[j-1], p[j], p[j-1]
   ierr   error flag (nonzero if no starting vector can be found)

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

   BLAS		ddot, dcopy, daxpy
   USER		opb, store
   MISC		random

 ***********************************************************************/

double s_las2::startv(long n, double *wptr[])

{
   double rnm2, *r, t;
   long irand;
   long id, i;

   /* get initial vector; default is random */
   rnm2 = ddot(n, wptr[0], 1, wptr[0], 1);
   irand = 918273 + j;
   r = wptr[0];
   for (id = 0; id < 3; id++) {
      if (id > 0 || j > 0 || rnm2 == 0) 
	 for (i = 0; i < n; i++) r[i] = mrandom(&irand);
      dcopy(n, wptr[0], 1, wptr[3], 1);

      /* apply operator to put r in range (essential if m singular) */
      mopb(n, wptr[3], wptr[0]);
      dcopy(n, wptr[0], 1, wptr[3], 1);
      rnm2 = ddot(n, wptr[0], 1, wptr[3], 1);
      if (rnm2 > 0.0) break;
   }

   /* fatal error */
   if (rnm2 <= 0.0) {
      ierr = 8192;
      return(-1);
   }
   if (j > 0) {
      for (i = 0; i < j; i++) {
         store(n, RETRQ, i, wptr[5]);
	 t = -ddot(n, wptr[3], 1, wptr[5], 1);
	 daxpy(n, t, wptr[5], 1, wptr[0], 1);
      }

      /* make sure q[j] is orthogonal to q[j-1] */
      t = ddot(n, wptr[4], 1, wptr[0], 1);
      daxpy(n, -t, wptr[2], 1, wptr[0], 1);
      dcopy(n, wptr[0], 1, wptr[3], 1);
      t = ddot(n, wptr[3], 1, wptr[0], 1);
      if (t <= eps * rnm2) t = 0.0;
      rnm2 = t;
   }
   return(sqrt(rnm2));
}



/***********************************************************************
 *                                                                     *
 *				imtqlb()			       *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------

   imtqlb() is a translation of a Fortran version of the Algol
   procedure IMTQL1, Num. Math. 12, 377-383(1968) by Martin and 
   Wilkinson, as modified in Num. Math. 15, 450(1970) by Dubrulle.  
   Handbook for Auto. Comp., vol.II-Linear Algebra, 241-248(1971).  
   See also B. T. Smith et al, Eispack Guide, Lecture Notes in 
   Computer Science, Springer-Verlag, (1976).

   The function 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	             

   (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.					    
   ierr   set to zero for normal return, j if the j-th eigenvalue has
            not been determined after 30 iterations.		    

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

   UTILITY	fsign
   MISC		pythag

 ***********************************************************************/

void s_las2::imtqlb(long n, double d[], double e[], double bnd[])

{
   long last, l, m, i, iteration;

   /* various flags */
   long exchange, convergence, underflow;	

   double b, test, g, r, s, c, p, f;

   if (n == 1) return;
   ierr = 0;
   bnd[0] = 1.0;
   last = n - 1;
   for (i = 1; i < n; i++) {
      bnd[i] = 0.0;
      e[i-1] = e[i];
   }
   e[last] = 0.0;
   for (l = 0; l < n; l++) {
      iteration = 0;
      while (iteration <= 30) {
	 for (m = l; m < n; m++) {
	    convergence = FALSE;
	    if (m == last) break;
	    else {
	       test = fabs(d[m]) + fabs(d[m+1]);
	       if (test + fabs(e[m]) == test) convergence = TRUE;
	    }
	    if (convergence) break;
	 }
	    p = d[l]; 
	    f = bnd[l]; 
	 if (m != l) {
	    if (iteration == 30) {
	       ierr = l;
	       return;
	    }
	    iteration += 1;
	    /*........ form shift ........*/
	    g = (d[l+1] - p) / (2.0 * e[l]);