las1.c

求矩阵奇异分解svd算法

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

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

{
   fprintf(fp_out1, " ... \n");
   fprintf(fp_out1, " ... SOLVE THE CYCLIC   EIGENPROBLEM\n");
   fprintf(fp_out1, " ... NO. OF EQUATIONS          =%5d\n", n);
   fprintf(fp_out1, " ... MAX. NO. OF LANCZOS STEPS =%5d\n", lanmax);
   fprintf(fp_out1, " ... MAX. NO. OF EIGENPAIRS    =%5d\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 DOCUMENTS (ROWS)   = %8d\n", nrow);
 fprintf(fp_out1, " ... NO. OF TERMS     (COLS)   = %8d\n", ncol);
 fprintf(fp_out1, " ... ORDER OF MATRIX A         = %8d\n", n);
 fprintf(fp_out1, " ... \n");
 return;
}
#include <stdio.h>
#include <math.h>
#include <errno.h>
#include <fcntl.h>
extern float eps, eps1, reps, eps34, *xv1, *xv2;
extern int j;
void   machar(int *, int *, int *, int *, int *);
int   check_parameters(int, int, int, float, float, int,
			int);
float dmax(float, float);
void   lanso(int, int, int, float, float, float *, float *,
	     float *[]);
void   ritvec(int, float, float *, float *, float *, float *,
              float *, float *);

/***********************************************************************
 *                                                                     *
 *				landr()				       *
 *        Lanczos algorithm with selective orthogonalization           *
 *                    Using Simon's Recurrence                         *
 *                       (float precision)                            *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   landr() is the LAS1 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 matrix B                 
   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 lav1; 
	    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 las1.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
   LAS1         ritvec, lanso

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

int landr(int n, int lanmax, int maxprs, int nnzero, float endl,
	   float endr, int vectors, float kappa, float *ritz,
	   float *bnd, float *r)

{
   int i, size, ibeta, it, irnd, machep, negep;
   float *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( (float) n );
   reps = sqrt(eps);
   eps34 = reps * sqrt(reps);


   size = 5*n + (lanmax*4 + 1);
   tptr = NULL;

   /* allocate work area and initialize pointers         *
    * pointer            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	 */

   if (!(tptr = (float *) malloc(size * sizeof(float)))){
      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 = (float *) malloc(n * (j+1) * sizeof(float))) ||
	  !(xv2 = (float *) malloc(n * sizeof(float)))) {
	  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 int ierr, j, nsig, fp_out2;
extern float *xv1;
void   dscal(int, float, float *,int);
void   dcopy(int, float *, int, float *, int);
void   daxpy(int, float, float *,int, float *, int);
void   store(int, int, int, float *);
void   imtql2(int, int, float *, float *, float *);
/***********************************************************************
 *                                                                     *
 *                        ritvec()                                     *
 * 	    Function computes the singular vectors of matrix A	       *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

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


   Parameters
   ----------

   (input)
   n	      dimension of the eigenproblem for matrix B
   kappa      relative accuracy of ritz values acceptable as 
		eigenvalues of B (singular values of 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
   fp_out2    pointer to unformatted output file
   j	      number of Lanczos iterations performed

   (output)
   xv1        array of eigenvectors of B (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(int n, float kappa, float *ritz, float *bnd,
	    float *alf, float *bet, float *w1, float *w2)

{
   int js, jsq, i, k, size, id, id2, tmp;
   float *s;

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

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

   /* initialize s as 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);

   /* compute eigenvalues of T */
   imtql2(js, js, w1, w2, s);
   if (ierr) return;

   /* on return w1 contains eigenvalues in ascending order and s 
    * contains the corresponding eigenvectors */

   write(fp_out2, (char *)&n, sizeof(n));
   write(fp_out2, (char *)&js, sizeof(js));
   write(fp_out2, (char *)&kappa, sizeof(kappa));
   nsig = 0;
   id = 0;
   id2 = jsq-js;
   for (k = 0; k < js; k++) {
      tmp = id2;
      if (bnd[k] <= kappa * fabs(ritz[k]) ) {
         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 sequentially in array xv1 */
         for (i = 0; i < n; i++) xv1[id++] = w1[i];
         nsig += 1;
      }
      id2++;
   }
   free(s);
   return;
}
extern int ncol, nrow, mxvcount;
extern int *pointr, *rowind;
extern float *value;

/************************************************************** 
 * multiplication of 2-cyclic matrix B by a vector x, where   *
 *							      *
 * B = [0  A]						      *
 *     [A' 0] , where A is nrow by ncol (nrow >> ncol). Hence,*
 * B is of order n = nrow+ncol (y stores product vector).     *
 **************************************************************/ 
void opb(int n,float *x, float *y)

{
   int i, j, end;
   float *tmp;
   
   mxvcount += 2;
   for (i = 0; i < n; i++) y[i] = 0.0;

   tmp = &x[nrow]; 
   for (i = 0; i < ncol; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++) 
	 y[rowind[j]] += value[j] * (*tmp); 
      tmp++; 
   }
   for (i = nrow; i < n; i++) {
      end = pointr[i-nrow+1];
      for (j = pointr[i-nrow]; j < end; j++) 
	 y[i] += value[j] * x[rowind[j]];
   }
   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 float rnm, anorm, tol, eps, eps1, reps, eps34;
extern int ierr, j, neig;
void   stpone(int, float *[]);
void   error_bound(int *, float, float, float *, float *);
void   lanczos_step(int, int, int, float *[], float *,
                    float *, float *, float *, int *, int *);
int   imin(int, int);
int   imax(int, int);
void   dsort2(int, int, float *, float *);
void   imtqlb(int n, float d[], float e[], float bnd[]);

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

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

   Function determines when the restart of the Lanczos algorithm should 
   occur and when it should terminate.

   Arguments 
   ---------

   (input)
   n         dimension of the eigenproblem for matrix B
   lanmax    upper limit of desired number of lanczos steps           
   maxprs    upper limit of desired number of eigenpairs             
   endl      left end of interval containing unwanted eigenvalues
   endr      right end of interval containing unwanted eigenvalues
   ritz      array to hold the ritz values                       
   bnd       array to hold the error bounds                          
   wptr      array of pointers that point to work space:            
  	       wptr[0]-wptr[5]  six vectors of length n		
  	       wptr[6] array to hold diagonal of the tridiagonal matrix T
  	       wptr[9] array to hold off-diagonal of T	
  	       wptr[7] orthogonality estimate of Lanczos vectors at 
		 step j
 	       wptr[8] orthogonality estimate of Lanczos vectors at 
		 step j-1

   (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
   ierr      (globally declared) error flag
	     ierr = 8192 if stpone() fails to find a starting vector
	     ierr = k if convergence did not occur for k-th eigenvalue
		    in imtqlb()
	     ierr = 0 otherwise


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

   LAS		stpone, error_bound, lanczos_step
   MISC		dsort2
   UTILITY	imin, imax

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

void lanso(int n, int lanmax, int maxprs, float endl,
	   float endr, float *ritz, float *bnd, float *wptr[])

{
   float *r, *alf, *eta, *oldeta, *bet, *wrk;
   int ll, first, last, ENOUGH, id1, id2, id3, i, l;