las2.c

求矩阵奇异分解svd算法

{
   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(fp_out1, "%s\n", error[error_index]);
   return(error_index);
}

extern FILE *fp_out1;
/***********************************************************************
 *								       *
 *			  write_data()				       *
 *   Function writes out header of output file containing ritz values  *
 *								       *
 ***********************************************************************/

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

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

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

 fprintf(fp_out1, " ... KAPPA                     =%10.2E\n", kappa);
 fprintf(fp_out1, " %s\n", title);
 fprintf(fp_out1, "           %s\n", name);
 fprintf(fp_out1, " ... NO. OF TERMS     (ROWS)   = %8ld\n", nrow);
 fprintf(fp_out1, " ... NO. OF DOCUMENTS (COLS)   = %8ld\n", ncol);
 fprintf(fp_out1, " ... ORDER OF MATRIX A         = %8ld\n", n);
 fprintf(fp_out1, " ... \n");
 return;
}
#include <stdio.h>
#include <math.h>
#include <errno.h>
#include <fcntl.h>

extern double eps, eps1, reps, eps34, *xv1, *xv2;
extern long nrow, ncol, j;
void   machar(long *, long *, long *, long *, long *);
long   check_parameters(long, long, long, double, double, long,
			long);
double dmax(double, double);
void   lanso(long, long, long, double, double, double *, double *,
	     double *[]);
void   ritvec(long, double, double *, double *, double *, double *,
	      double *, double *);

/***********************************************************************
 *                                                                     *
 *				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 landr(long n, long lanmax, long maxprs, long nnzero, double endl, 
	   double endr, long 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      */

   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(n, kappa, ritz, bnd, wptr[6], wptr[9], wptr[4], wptr[5]);
   }

   free(tptr2);
   return(0);
}

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

extern long nrow, ierr, j, fp_out2, nsig, neig;
extern double *xv1;

void   dscal(long, double, double *,long);
void   dcopy(long, double *, long, double *, long);
void   daxpy(long, double, double *,long, double *, long);
void   store(long, long, long, double *);
void   imtql2(long, long, double *, double *, double *);

/***********************************************************************
 *                                                                     *
 *                        ritvec()                                     *
 * 	    Function computes the singular vectors of matrix A	       *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   This function is invoked by landr() only if eigenvectors of the A'A
   eigenproblem are desired.  When called, ritvec() computes the 
   singular vectors of A and writes the result to an unformatted file.


   Parameters
   ----------

   (input)
   nrow       number of rows of A
   j	      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

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

   BLAS		dscal, dcopy, daxpy
   USER		store
   		imtql2

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

void ritvec(long n, double kappa, double *ritz, double *bnd, double *alf,
	    double *bet, double *w1, double *w2)

{
   long js, jsq, i, k, size, id, id2, tmp;
   double *s;

   js = j + 1;
   jsq = js * js;
   size = sizeof(double) * n;

   if(!(s = (double *) malloc (jsq * sizeof(double)))) {
      perror("MALLOC FAILED in RITVEC()");
      exit(errno);
   }

   /* initialize s to an identity matrix */
   for (i = 0; i < jsq; i++) s[i] = 0.0;
   for (i = 0; i < jsq; i+= (js+1)) s[i] = 1.0;
   dcopy(js, alf, 1, w1, -1);
   dcopy(j, &bet[1], 1, &w2[1], -1);

   /* on return from imtql2(), w1 contains eigenvalues in ascending 
    * order and s contains the corresponding eigenvectors */
   imtql2(js, js, w1, w2, s);
   if (ierr) return;

   write(fp_out2, (char *)&n, sizeof(n));
   write(fp_out2, (char *)&js, sizeof(js));
   write(fp_out2, (char *)&kappa, sizeof(kappa));
   id = 0;
   nsig = 0;
   id2 = jsq - js;
   for (k = 0; k < js; k++) {
      tmp = id2;
      if (bnd[k] <= kappa * fabs(ritz[k]) && k > js-neig-1) {
	 for (i = 0; i < n; i++) w1[i] = 0.0;
         for (i = 0; i < js; i++) {
	    store(n, RETRQ, i, w2);
	    daxpy(n, s[tmp], w2, 1, w1, 1);
	    tmp -= js;
         }
         write(fp_out2, (char *)w1, size);

      /* store the w1 vector row-wise in array xv1;   
       * size of xv1 is (j+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 (j+1) by     
       * (nrow+ncol) and each vector occupies a row  */

         for (i = 0; i < n; i++) xv1[id++] = w1[i];
         id += nrow;
         nsig += 1;
      }
      id2++;
   }
   free(s);
   return;
}

extern long ncol, nrow,mxvcount;
extern long *pointr, *rowind;
extern double *value,*ztemp;

/**************************************************************
 * multiplication of matrix B by vector x, where B = A'A,     *
 * and A is nrow by ncol (nrow >> ncol). Hence, B is of order *
 * n = ncol (y stores product vector).		              *
 **************************************************************/

void opb(long n, double *x, double *y)

{
   long i, j, end;
   
   mxvcount += 2;
   for (i = 0; i < n; i++) y[i] = 0.0;
   for (i = 0; i < nrow; i++) ztemp[i] = 0.0;

   for (i = 0; i < ncol; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++) 
	 ztemp[rowind[j]] += value[j] * (*x); 
      x++;
   }
   for (i = 0; i < ncol; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++) 
	 *y += value[j] * ztemp[rowind[j]];
      y++;
   }
   return;
}

/***********************************************************
 * multiplication of matrix A by vector x, where A is 	   *
 * nrow by ncol (nrow >> ncol).  y stores product vector.  *
 ***********************************************************/

void opa(double *x, double *y)

{
   long end, i, j;
   
   mxvcount += 1;
   for (i = 0; i < nrow; i++) y[i] = 0.0;

   for (i = 0; i < ncol; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)
	 y[rowind[j]] += value[j] * x[i]; 
   }
   return;
}
#include <stdio.h>
#include <math.h>
#define STORQ 1
#define RETRQ 2
#define STORP 3
#define RETRP 4
#define TRUE  1
#define FALSE 0

extern double rnm, anorm, tol, eps, eps1, reps, eps34;
extern long ierr, j, neig;
void   stpone(long, double *[]);
void   error_bound(long *, double, double, double *, double *);
void   lanczos_step(long, long, long, double *[], double *,
                    double *, double *, double *, long *, long *);
long   imin(long, long);
long   imax(long, long);
void   dsort2(long, long, double *, double *);
void   imtqlb(long n, double d[], double e[], double bnd[]);

/***********************************************************************
 *                                                                     *
 *                          lanso()                                    *
 *                                                                     *
 ***********************************************************************/