las2.c

求矩阵奇异分解svd算法

/*************************************************************************
                           (c) Copyright 1993
                        University of Tennessee
                          All Rights Reserved                          
 *************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <fcntl.h>
#include "las2.h"

/* #define  UNIX_CREAT */

#ifdef UNIX_CREAT
#define PERMS 0664
#endif

long   landr           (long, long, long, long, double, double, long,
	                double, double *, double *, double *);
void   dscal           (long, double, double *,long);
double ddot            (long, double *,long, double *, long);
void   daxpy           (long, double, double *,long, double *, long);
void   opb             (long, double *, double *);
void   opa             (double *, double *);
void   write_data      (long, long, double, double, long, double,
		        char *,char *, long, long, long);
long   check_parameters(long, long, long, double, double, long,
		        long);
float  timer           (void);

/***********************************************************************
 *                                                                     *
 *                        main()                                       *
 * Sparse SVD(A) via Eigensystem of A'A symmetric Matrix 	       *
 *                  (double precision)                                 *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   This sample program uses landr to compute singular triplets of A via
   the equivalent symmetric eigenvalue problem                         

   B x = lambda x, where x' = (u',v'), lambda = sigma**2,
   where sigma is a singular value of A,
                                                                     
   B = A'A , and A is m (nrow) by n (ncol) (nrow >> ncol),                
                                                                 
   so that {u,sqrt(lambda),v} is a singular triplet of A.        
   (A' = transpose of A)                                      
                                                            
   User supplied routines: opa, opb, store, timer              
                                                        
   opa(     x,y) takes an n-vector x and returns A*x in y.
   opb(ncol,x,y) takes an n-vector x and returns B*x in y.
                                                                  
   Based on operation flag isw, store(n,isw,j,s) stores/retrieves 
   to/from storage a vector of length n in s.                   
                                                               
   User should edit timer() with an appropriate call to an intrinsic
   timing routine that returns elapsed user time.                      


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

   Defined and documented in las2.h


   Local parameters 
   ----------------

  (input)
   endl     left end of interval containing unwanted eigenvalues of B
   endr     right end of interval containing unwanted eigenvalues of B
   kappa    relative accuracy of ritz values acceptable as eigenvalues
              of B
	      vectors is not equal to 1
   r        work array
   n	    dimension of the eigenproblem for matrix B (ncol)
   maxprs   upper limit of desired number of singular triplets of A
   lanmax   upper limit of desired number of Lanczos steps
   nnzero   number of nonzeros in A
   vectors  1 indicates both singular values and singular vectors are 
	      wanted and they can be found in output file lav2;
	      0 indicates only singular values are wanted 
   		
  (output)
   ritz	    array of ritz values
   bnd      array of error bounds
   d        array of singular values
   memory   total memory allocated in bytes to solve the B-eigenproblem


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

   BLAS		daxpy, dscal, ddot
   USER		opa, opb, timer
   MISC		write_data, check_parameters
   LAS2		landr


   Precision
   ---------

   All floating-point calculations are done in double precision;
   variables are declared as long and double.


   LAS2 development
   ----------------

   LAS2 is a C translation of the Fortran-77 LAS2 from the SVDPACK
   library written by Michael W. Berry, University of Tennessee,
   Dept. of Computer Science, 107 Ayres Hall, Knoxville, TN, 37996-1301

   31 Jan 1992:  Date written 

   Theresa H. Do
   University of Tennessee
   Dept. of Computer Science
   107 Ayres Hall
   Knoxville, TN, 37996-1301
   internet: [e-mail]

 ***********************************************************************/
		    
main()

{
   float t0,exetime;
   double endl, endr, kappa, tmp0, tmp1, xnorm;
   double *r, *ritz, *bnd, *d, *tptr1;
   long nn, k, i, id, ida, n, lanmax, maxprs, nnzero;
   long memory, vectors, size1, size2, count1, count2;
   char title[73], name[41], v[6];
   char *in1, *in2, *out1, *out2;
   FILE *fp_in1, *fp_in2;

   in1 = "lap2";
   in2 = "matrix";
   out1 = "lao2";
   out2 = "lav2";

   /* open files for input/output */
    if (!(fp_in1 = fopen(in1, "r"))) { 
       printf("cannot open file %s for reading\n", in1);
       exit(-1);
    }
    if (!(fp_in2 = fopen(in2, "r"))) {
       printf("cannot open file %s for reading\n", in2);
       exit(-1);
    }
    if (!(fp_out1 = fopen(out1, "w"))) { 
       printf("cannot open output file %s \n", out1);
       exit(-1);
    }

    /* read data */
    fscanf (fp_in2, "%72c%*s%*s%*s%ld%ld%ld%*d", 
	    title, &nrow, &ncol, &nnzero);
    title[73] = '\0';
    fscanf (fp_in1, "%s %ld %ld %lf %lf %s %lf", name, &lanmax, 
            &maxprs, &endl, &endr, v, &kappa);

    if (!(strcmp(v, "TRUE"))) {
       vectors = 1;

#if  !defined UNIX_CREAT

       if ((fp_out2 = open(out2, O_CREAT | O_RDWR)) == -1) {
          printf("cannot open output file %s \n", out2);
          exit(-1);
       }
#else
       if ((fp_out2 = creat(out2, PERMS)) == -1) {
          printf("cannot open output file %s \n", out2);
          exit(-1);
       }
#endif

    }
    else vectors = 0;

    n = ncol;
    nn = ncol + nrow;

    /* write header of output file */
    write_data(lanmax, maxprs, endl, endr, vectors, kappa, title, 
	       name, nrow, ncol, n);

    /* even though the validity of the parameters will be checked in the
     * SVD code, some parameter checking should also be done before 
     * allocating memory to ensure that they are nonnegative */ 
    if (check_parameters(maxprs, lanmax, n, endl, endr, vectors, nnzero)) {
       fclose(fp_in1);
       fclose(fp_in2);
       fclose(fp_out1);
       if (vectors) close(fp_out2);
       exit(-1);
    }

    /*******************************************************************
     * allocate memory                                                 *
     * pointr - column start array of harwell-boeing sparse matrix     *
     *          format                                       (ncol+1)  *
     * rowind - row indices array of harwell-boeing format   (nnzero)  *
     * value  - nonzero values array of harwell-boeing sparse matrix   *
     *          format                                       (nnzero)  *
     * r      - work array                                        (n)  *
     * ritz   - array of ritz values                              (n)  *
     * bnd    - array of error bounds                             (n)  *
     * d      - array of approximate singular values of matrix A  (n)  *
     * ztemp  - work array for user function opb	       (nrow)  *
     * a      - storage area for Lanczos vectors     (n * (lanmax+2))  *
     *******************************************************************/

    size1 = sizeof(double) * (6 * n + nrow + nnzero + (n * lanmax));
    size2 = sizeof(long) * (ncol + 1 + nnzero);
    if (!(pointr = (long *)   malloc(size2))   ||
        !(value  = (double *) malloc(size1))){
	perror("MALLOC FAILED in MAIN()");
	exit(errno);
    }
    tptr1 = value;
    /* allocated memory including work array used in landr */
    memory = size1 + size2 + sizeof(double) * (5 * n + lanmax * 4 + 1);

    rowind = pointr + ncol + 1;
    tptr1 += nnzero;
    r      = tptr1;
    tptr1 += n;
    ritz   = tptr1;
    tptr1 += n;
    bnd    = tptr1;
    tptr1 += n;
    d      = tptr1;
    tptr1 += n;
    ztemp  = tptr1;
    tptr1 += nrow;
    a      = tptr1;

    /* skip data format line */
    fscanf(fp_in2, "%*s %*s %*s %*s");

    /* read data */
    for (i = 0; i <= ncol; i++) fscanf(fp_in2, "%ld", &pointr[i]);
    for (i = 0; i < ncol; i++) pointr[i] -= 1;

    /* define last element of pointr in case it is not */
    pointr[i] = nnzero;

    for (i = 0; i < nnzero; i++) fscanf(fp_in2, "%ld", &rowind[i]);
    for (i = 0; i < nnzero; i++) rowind[i] -= 1;
    for (i = 0; i < nnzero; i++) fscanf(fp_in2, "%lf", &value[i]);

    /* to get a random starting vector, the first n cells must be
     * initialize to zero */
    for (i = 0; i < n; i++) r[i] = 0.;

    exetime = timer();

    /* make a lanczos run; see landr for meaning of parameters */
    if(landr(n, lanmax, maxprs, nnzero, endl, endr, vectors, kappa, 
       ritz, bnd, r)){

       free(value);
       free(pointr);
       fclose(fp_in1);
       fclose(fp_in2);
       fclose(fp_out1);
       if (vectors) {
	  close(fp_out2);
	  free(xv1);
	  free(xv2);
       }
       exit(-1);
    }

    exetime = timer() - exetime;

    /* memory allocated for xv1, xv2 and s in landr() */
    if (vectors) memory += sizeof(double) * (nn * (j+1) + n + 
					     (j+1) * (j+1));

    /* print error code if not zero */
    if (ierr)fprintf(fp_out1, " ... RETURN FLAG = %9ld ...\n", ierr);

    /* print ritz values and error bounds */
    fprintf(fp_out1, "\n");
    fprintf(fp_out1, " ...... ALLOCATED MEMORY (BYTES)= %10.2E\n", (float)memory);
    fprintf(fp_out1, " ...... LANSO EXECUTION TIME=%10.2E\n", exetime);
    fprintf(fp_out1, " ...... \n");
	fprintf(fp_out1, " ...... NUMBER OF LANCZOS STEPS = %3ld       NEIG = %3ld\n", j+1, neig);
    fprintf(fp_out1, " ...... \n");
	fprintf(fp_out1, " ......         COMPUTED RITZ VALUES  (ERROR BNDS)\n");
    fprintf(fp_out1, " ...... \n");
    for (i = 0; i <= j; i++)
       fprintf(fp_out1, " ...... %3ld   %22.14E  (%11.2E)\n",
	       i + 1, ritz[i], bnd[i]);

    /* compute residual error when singular values and vectors are 
     * computed.  This is done only for singular values that are
     * within the relative accuracy (kappa) */
    if (vectors) {
       size1 = sizeof(double) * nrow;
       t0 = timer();
       id = 0;

       for (i = 0; i < nsig; i++) {

	  /* multiply by matrix B first */
	  opb(n, &xv1[id], xv2);
	  tmp0 = ddot(n, &xv1[id], 1, xv2, 1);
	  daxpy(n, -tmp0, &xv1[id], 1, xv2, 1);
	  tmp0 = sqrt(tmp0);
	  xnorm = sqrt(ddot(n, xv2, 1, xv2, 1));
	  ida = id + ncol;

	  /* multiply by matrix A to get (scaled) left s-vector */
	  opa(&xv1[id], &xv1[ida]);
	  tmp1 = 1.0 / tmp0;
	  dscal(nrow, tmp1, &xv1[ida], 1);
	  xnorm *= tmp1;
	  bnd[i] = xnorm;
	  d[i] = tmp0;

	  /* write left s-vector to output file */
	  write(fp_out2, (char *)&xv1[ida], size1);
	  id += nn;
       }
       exetime += (timer() - t0);
       count1=(mxvcount-nsig)/2 + nsig;
       count2=(mxvcount-nsig)/2;
       fprintf(fp_out1, " ...... \n");
       fprintf(fp_out1, " ...... NO. MULTIPLICATIONS BY A  =%10ld\n", count1);
       fprintf(fp_out1, " ...... NO. MULT. BY TRANSPOSE(A) =%10ld\n", count2);
       fprintf(fp_out1, "\n");
       fprintf(fp_out1, " ...... LASVD EXECUTION TIME=%10.2E\n", exetime);
       fprintf(fp_out1, " ...... \n");
       fprintf(fp_out1, " ......        NSIG = %4ld\n", nsig);
       fprintf(fp_out1, " ...... \n");
       fprintf(fp_out1, " ......         COMPUTED S-VALUES     (RES. NORMS)\n");
       fprintf(fp_out1, " ...... \n");
       for (i = 0; i < nsig; i++)
          fprintf(fp_out1, " ...... %3ld   %22.14E  (%11.2E)\n",
		  i + 1, d[i], bnd[i]);
    }
    else {
       for (i = j; i >= 0; i--)
	  if (bnd[i] > kappa * fabs(ritz[i])) break;
       nsig = j - i;

       count1=(mxvcount-nsig)/2 + nsig;
       count2=(mxvcount-nsig)/2;
       fprintf(fp_out1," ...... \n");
       fprintf(fp_out1," ...... NO. MULTIPLICATIONS BY A  =%10ld\n", count1);
       fprintf(fp_out1," ...... NO. MULT. BY TRANSPOSE(A) =%10ld\n", count2);
       fprintf(fp_out1, "\n");
       fprintf(fp_out1," ...... LASVD EXECUTION TIME = %10.2E\n", exetime);
       fprintf(fp_out1," ...... \n");
       fprintf(fp_out1," ......         NSIG = %4ld\n" , nsig);
       fprintf(fp_out1," ...... \n");
       fprintf(fp_out1," ......         COMPUTED S-VALUES   (ERROR BNDS)\n");
       fprintf(fp_out1," ...... \n");

       k = j + 1 - nsig;
       for (i = 1 ; i <= nsig; i++) {
	  fprintf(fp_out1," ...... %3ld   %22.14E  (%11.2E)\n", 
	           i, sqrt(ritz[k]), bnd[k]);
          k++;
          }
   }
    free(value);
    free(pointr);
    fclose(fp_in1);
    fclose(fp_in2);
    fclose(fp_out1);
    if (vectors) {
       free(xv1);
       free(xv2);
       close(fp_out2);
    }
    exit(0);
}

extern long ncol,nrow;
extern char *error[];
extern FILE *fp_out1;
/***********************************************************************
 *								       *
 *		      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 lav2; 0 indicates eigenvalues only
   nnzero   number of nonzero elements in input matrix (matrix A)      
                                                                      
 ***********************************************************************/

long check_parameters(long maxprs, long lanmax, long n, 
		      double endl, double endr, long vectors, 
		      long nnzero)