las1.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 "las1.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);
void   daxpy(long, double, double *,long, double *, long);
double ddot(long, double *,long, double *, long);
void   opb(long, 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 Via Eigensystem of Equivalent 2-Cyclic 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, 
   where sigma is a singular value of A,
                                                                     
       [o  A]	
   B = [    ] , and A is m (nrow) by n (ncol) (nrow >> ncol),
       [A' o]
                                                                 
   so that {u, abs(lambda), v} is a singular triplet of A.        
   (A' = transpose of A)                                      
                                                            
   User supplied routines:  opb, store, timer              

   opb(nrow+ncol,x,y) takes an (m+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 cpu time.


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

   Defined and documented in las1.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
   r        work array
   n	    dimension of the eigenproblem for matrix B (nrow + 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 lav1; 
	    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 B-eigenproblem


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

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


   Precision
   ---------

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


   LAS1 development
   ----------------

   LAS1 is a C translation of the Fortran-77 LAS1 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, tmp2, tmp3, xnorm;
   double *r, *ritz, *bnd, *d, *tptr1;
   long k, i, id, n, lanmax, maxprs, nnzero, size1, size2, memory, vectors;
   long *tptr3, hcount;
   char title[73], name[41], v[6];
   char *in1, *in2, *out1, *out2;
   FILE *fp_in1, *fp_in2;
 
   in1 = "lap1";
   in2 = "matrix";
   out1 = "lao1";
   out2 = "lav1";

   /* 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 = nrow + ncol;

    /* 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)  *
     * a      - storage area for Lanczos vectors     (n * (lanmax+2))  *
     *******************************************************************/
	 size1 = sizeof(double) * (n * (lanmax + 6) + nnzero);
	 size2 = sizeof(long) * (ncol + nnzero + 1);
	 if (!(pointr = (long *)   malloc(size2))   ||
	     !(value  = (double *) malloc(size1))){
	     perror("MALLOC FAILED in MAIN()");
	     exit(errno);
         }
	 tptr1 = value;
	 
    /* calculate memory allocated for problem (including work space in landr */
    memory = size1 + size2 + sizeof(double) * (5 * n + 4 * lanmax + 1);	
     
     rowind = pointr + (ncol + 1);
     tptr1 += nnzero;
     r      = tptr1;
     tptr1 += n;
     ritz   = tptr1;
     tptr1 += n;
     bnd    = tptr1;
     tptr1 += n;
     d      = tptr1;
     tptr1 += n;
     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; exit upon error */
    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;

    /* calculating singular vectors requires extra work space
     * (xv1, xv2, s) in landr() */
    if (vectors) memory += sizeof(double) * 
			   (n*(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 */
    if (vectors) {
       t0 = timer();
       id = 0;

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

	  /* multiply by 2-cyclic indefinite matrix */
	  opb(n, &xv1[id], xv2);
	  tmp0 = ddot(n, &xv1[id], 1, xv2, 1);
	  tmp1 = ddot(nrow, &xv1[id], 1, &xv1[id], 1);
	  tmp2 = ddot(ncol, &xv1[id + nrow], 1, &xv1[id + nrow], 1);
	  tmp3 = sqrt(tmp1 + tmp2);
	  daxpy(n, -tmp0, &xv1[id], 1, xv2, 1);
	  xnorm = ddot(n, xv2, 1, xv2, 1);
	  xnorm = sqrt(xnorm);
	  bnd[i] = xnorm / tmp3;
	  d[i] = fabs(tmp0);
	  id += n;
       }  
       exetime += (timer() - t0);
       hcount=mxvcount/2;
       fprintf(fp_out1," ...... \n");
       fprintf(fp_out1," ...... NO. MULTIPLICATIONS BY A  =%10ld\n", hcount);
       fprintf(fp_out1," ...... NO. MULT. BY TRANSPOSE(A) =%10ld\n", hcount);
       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;
       hcount=mxvcount/2;
       fprintf(fp_out1," ...... \n");
       fprintf(fp_out1," ...... NO. MULTIPLICATIONS BY A  =%10ld\n", hcount);
       fprintf(fp_out1," ...... NO. MULT. BY TRANSPOSE(A) =%10ld\n", hcount);
       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, 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 lav1; 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) 
{
   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]);