las.cc

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

// File: las.cc -- Implements the las base class
// Author: Suvrit Sra
// Date: 14 Nov, 2003

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

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cerrno>

#include <unistd.h>

#include "las.h"

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

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

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

   This function is a partial translation of a Fortran-77 subroutine 
   written by W. J. Cody of Argonne National Laboratory.
   It dynamically determines the listed machine parameters of the
   floating-point arithmetic.  According to the documentation of
   the Fortran code, "the determination of the first three uses an
   extension of an algorithm due to M. Malcolm, ACM 15 (1972), 
   pp. 949-951, incorporating some, but not all, of the improvements
   suggested by M. Gentleman and S. Marovich, CACM 17 (1974), 
   pp. 276-277."  The complete Fortran version of this translation is
   documented in W. J. Cody, "Machar: a Subroutine to Dynamically 
   Determine Determine Machine Parameters," TOMS 14, December, 1988.


   Parameters reported 
   -------------------

   ibeta     the radix for the floating-point representation       
   it        the number of base ibeta digits in the floating-point
               significand					 
   irnd      0 if floating-point addition chops		      
             1 if floating-point addition rounds, but not in the 
                 ieee style					
             2 if floating-point addition rounds in the ieee style
             3 if floating-point addition chops, and there is    
                 partial underflow				
             4 if floating-point addition rounds, but not in the
                 ieee style, and there is partial underflow    
             5 if floating-point addition rounds in the ieee style,
                 and there is partial underflow                   
   machep    the largest negative integer such that              
                 1.0+float(ibeta)**machep .ne. 1.0, except that 
                 machep is bounded below by  -(it+3)          
   negeps    the largest negative integer such that          
                 1.0-float(ibeta)**negeps .ne. 1.0, except that 
                 negeps is bounded below by  -(it+3)	       

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

void las::machar(long *ibeta, long *it, long *irnd, long *machep, long *negep)

{

   double beta, betain, betah, a, b, ZERO, ONE, TWO, temp, tempa, temp1;
   long i, itemp;

   ONE = (double) 1;
   TWO = ONE + ONE;
   ZERO = ONE - ONE;

   a = ONE;
   temp1 = ONE;
   while (temp1 - ONE == ZERO) {
      a = a + a;
      temp = a + ONE;
      temp1 = temp - a;
   }

   b = ONE;
   itemp = 0;
   while (itemp == 0) {
      b = b + b;
      temp = a + b;
      itemp = (long)(temp - a);
   }
   *ibeta = itemp;
   beta = (double) *ibeta;

   *it = 0;
   b = ONE;
   temp1 = ONE;
   while (temp1 - ONE == ZERO) {
      *it = *it + 1;
      b = b * beta;
      temp = b + ONE;
      temp1 = temp - b;
   }
   *irnd = 0; 
   betah = beta / TWO; 
   temp = a + betah;
   if (temp - a != ZERO) *irnd = 1;
   tempa = a + beta;
   temp = tempa + betah;
   if ((*irnd == 0) && (temp - tempa != ZERO)) *irnd = 2;


   *negep = *it + 3;
   betain = ONE / beta;
   a = ONE;
   for (i = 0; i < *negep; i++) a = a * betain;
   b = a;
   temp = ONE - a;
   while (temp-ONE == ZERO) {
      a = a * beta;
      *negep = *negep - 1;
      temp = ONE - a;
   }
   *negep = -(*negep);

   *machep = -(*it) - 3;
   a = b;
   temp = ONE + a;
   while (temp - ONE == ZERO) {
      a = a * beta;
      *machep = *machep + 1;
      temp = ONE + a;
   }
   eps = a;
   return;
}


/***********************************************************************
 *								       *
 *		      check_parameters()			       *
 *								       *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------
   Function validates input parameters and returns error code (long)  

   Parameters 
   ----------
  (input)
   maxprs   upper limit of desired number of eigenpairs of B           
   lanmax   upper limit of desired number of lanczos steps             
   n        dimension of the eigenproblem for matrix B               
   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 lav1; 0 indicates eigenvalues only
   nnzero   number of nonzero elements in input matrix (matrix A)      
                                                                      
 ***********************************************************************/

long las::check_parameters(long maxprs, long lanmax, long n,
		      double endl, double endr, long vectors,
		      long nnzero) 
{
   long error_index, ncells;
   error_index = 0;

   /* assuming that nrow >= ncol... */
   if (ncol >= NMAX || nnzero > NZMAX) error_index = 1;
   else if (endl >= endr)  error_index = 2;
   else if (maxprs > lanmax)  error_index = 3;
   else if (n <= 0)  error_index = 4;
   else if (lanmax <= 0 || lanmax > n)  error_index = 5;
   else if (maxprs <= 0 || maxprs > lanmax)  error_index = 6;
   else {
       if (vectors) {
	  ncells = 6 * n + 4 * lanmax + 1 + lanmax * lanmax;
	  if (ncells > LMTNW) error_index = 7;
       }
       else {
	  ncells = 6 * n + 4 * lanmax + 1;
	  if (ncells > LMTNW) error_index = 8;
       }
   }
   if (error_index)
     fprintf(stderr, "Error: %s\n", error[error_index]);
   //fprintf(stderr, "%f  %f\n", endl, endr);
   return(error_index);
}


/***********************************************************************
 *								       *
 *			  write_data()				       *
 *   Function writes out header of output file containing ritz values  *
 *								       *
 ***********************************************************************/

void las::write_data(FILE* fp, long lanmax, long maxprs, double endl, double endr,
		     bool vectors, double kappa, char *title,
		     char *name, long nrow, long ncol, long n)

{
   fprintf(fp, " ... \n");
   fprintf(fp, " ... SOLVE THE CYCLIC   EIGENPROBLEM\n");
   fprintf(fp, " ... NO. OF EQUATIONS          =%5ld\n", n);
   fprintf(fp, " ... MAX. NO. OF LANCZOS STEPS =%5ld\n", lanmax);
   fprintf(fp, " ... MAX. NO. OF EIGENPAIRS    =%5ld\n", maxprs);
   fprintf(fp, " ... LEFT  END OF THE INTERVAL =%10.2E\n", endl);
   fprintf(fp, " ... RIGHT END OF THE INTERVAL =%10.2E\n", endr);


 if (vectors) 
   fprintf(fp, " ... WANT S-VECTORS?   [T/F]   =   T\n");
 else
   fprintf(fp, " ... WANT S-VECTORS?   [T/F]   =   F\n");

 fprintf(fp, " ... KAPPA                     =%10.2E\n", kappa);
 fprintf(fp, " %s\n", title);
 fprintf(fp, "           %s\n", name);
 fprintf(fp, " ... NO. OF TERMS     (ROWS)   = %8ld\n", nrow);
 fprintf(fp, " ... NO. OF DOCUMENTS (COLS)   = %8ld\n", ncol);
 fprintf(fp, " ... ORDER OF MATRIX A         = %8ld\n", n);
 fprintf(fp, " ... \n");
 return;
}