bls2.c

svd 算法代码 This directory contains instrumented SVDPACKC Version 1.0 (ANSI-C) programs for compiling

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

#define  ONE   1.0

#define  UNIX_CREAT

#ifdef UNIX_CREAT
#define PERMS 0664
#endif


double ddot(long, double *,long, double *, long);
void   opb(long, long, double *, double *);
void   opa(long, long, double *, double *);
void   dscal(long, double, double *,long);
void   daxpy(long, double, double *,long, double *, long);
double  timer();
long   imin(long, long);
long   validate(FILE *, long, long, long, long, long, long, double);
long   blklan2(FILE *, long, long, long, long, double *, double *,
	       long, long, double, double *, long, long *, long *,
	       long *, long*);

/***********************************************************************
 *                                                                     *
 *                        main()                                       *
 *     Sparse SVD Via Hybrid Block Lanczos Procedure for Equivalent    *
 *               A'A Eigensystems  (double precision)                  *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   This is a sample driver that invokes BLKLAN2 to compute singular 
   triplets of a large sparse matrix A.  In this test program, bls2,
   the Harwell-Boeing sparse matrix format is employed for accessing
   elements of the m by n sparse matrix A and its transpose.  Other 
   sparse matrix formats can be used, of course.  Approximate singular
   values of A and their residuals can be found in formatted output
   file blo2 and corresponding right and left singular vectors can be
   found in unformatted output file blv2.

   User supplied routines:  opa, opm, opb, timer
   
   opa(m, n, x, y) takes an n-vector x and returns the m-vector y = A * x.

   opb(m, n, x, y) takes an n-vector x and returns the m-vector y = A'A * x.  

   opm(m, n, nc, x, y) takes a block of nc vectors length n stored in two-
   dimensional array x and returns the block of n-vectors y = A'A * x. 

   User should edit timer() with an appropriate call to an longrinsic
   timing routine that returns elapsed user cpu time.


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

   Defined and documented in bls2.h


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

  (input)
   maxit     maximum number of iterations
   nc        maximum subspace size
   nb        initial block size
   nums      number of singular values desired 
   tol       residual tolerance
   vtf       TRUE or FALSE to indicate whether both singular values and
		vectors are wanted or only singular values are wanted.
		If vtf is set to "TRUE", the unformatted output file blv2
		will contain the approximate singular vectors written in 
		the order u[0], v[0], u[1], v[1], ..., u[iko-1], v[iko-1].
		u[i], v[i] denote the left and right singular vectors,
		respectively, corresponding to the i-th approximate singular
		value sing[i], and iko is the number of singular values 
		approximated.

   (output)
   sing      array of singular values
   res       array of residual norms of the singular values calculated
   memory    total memory allocated in bytes


   Functions called
   --------------

   BLAS         ddot, daxpy, dscal
   USER         opa, opb, timer
   MISC         validate
   BLS2         blklan2


   BLS2 development
   ----------------

   BLS2 is a C translation of the Fortran-77 BLS2 from the SVDPACK
   library written by Michael W. Berry, University of Tennessee,
   Dept. of Computer Science, 107 Ayres Hall, Knoxville, TN, 37996-1301.
   This particular implementation is discussed in "Multiprocessor Sparse
   SVD Algorithms and Applications", Ph.D. thesis by M. Berry, University
   of Illinois at Urbana-Champaign, October 1990.

   Date written:  05 May 1992

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

 ***********************************************************************/
  main() 
{
   double t0, exetime;
   long i, j, k, nrow, ncol, nnzero, nc, nb, nums, maxit, nn;
   long ncold, nbold, nk, size1, size2, indexu, indexv;
   long memory, vectors;
   double tol, *sing, *v, *u, *res, *tmpv, *tptr1;
   double tmp1, tmp2, xnorm;
   char title[73], name[41], vtf[6];
   char *in1, *in2, *out1, *out2;
   FILE *fp_in1, *fp_in2;
   FILE *fp_out1 = NULL; 
   long  fp_out2;
 
   in1 = "blp2";
   in2 = "matrix";
   out1 = "blo2";
   out2 = "blv2";

   /* 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';
    nn=imin(nrow,ncol);
    fscanf (fp_in1,"%s%ld%ld%ld%ld%lf%s", name, &maxit, &nc, &nb, 
	    &nums, &tol, vtf);

    if (!(strcmp(vtf, "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;

    /* 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 (validate(fp_out1, nnzero, nrow, ncol, nc, nb, nums, tol)) {
	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)   *
     * sing   -                                               (nums)   *
     * res    -                                               (nums)   *
     * v      -                                        (ncol * nums)   *
     * ztemp  -                     (nrow * nc; assume nrow >= ncol)   *
     * (ztempp)                                                 (nc)   *
     *******************************************************************/
     size1 = sizeof(double) * (nnzero + nums * (ncol+2) + nrow * nc);
     size2 = sizeof(long) * (ncol + nnzero + 1);
     if (!(rowind = (long *)   malloc(size2))   ||
         !(value  = (double *) malloc(size1))   ||
	 !(ztempp = (double **) malloc(sizeof(double *) * nc)))    {
	 perror("MALLOC FAILED in MAIN()");
	 exit(errno);
     }
     tptr1 = value;

     /* calculate memory allocated for problem */
     memory = size1 + size2 + sizeof(double *) * nc;

     pointr = rowind + nnzero;
     tptr1 += nnzero;
     v      = tptr1;
     tptr1 += ncol * nums;
     sing   = tptr1;
     tptr1 += nums;
     res    = tptr1;
     tptr1 += nums;
     ztemp  = tptr1;
     j = 0;
     k = nc * nrow;
     for (i = 0; i < k; i += nrow) ztempp[j++] = &ztemp[i];
    /* 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]);

    ncold = nc;
    nbold = nb;
    nk    = nums;
    mxvcount = 0;
    mtxvcount = 0;

    exetime = timer();
    if (blklan2(fp_out1, nnzero, nrow, ncol, nums, v, sing, ncold, 
               nbold, tol, res, maxit, &nk, &nc, &nb, &memory)) {

       free(value);
       free(rowind);
       fclose(fp_in1);
       fclose(fp_in2);
       fclose(fp_out1);
       if (vectors) close(fp_out2);
       exit(-1);
    }

    exetime = timer() - exetime;
    if (vectors) {
    /*******************************************************************
     * allocate memory						       *
     * tmpv   -                                               (ncol)   *
     * u      -                                          (nrow * nk)   *
     *******************************************************************/
       size1 = sizeof(double) * (ncol + nk * nrow);
       memory += size1;
       if (!(tmpv = (double *) malloc(size1))) {
	   perror("SECOND MALLOC FAILED in MAIN()");
	   exit(errno);
       }
       u = tmpv + ncol;
       size1 = sizeof(double) * nrow;
       size2 = sizeof(double) * ncol;
       indexu = 0;
       indexv = 0;

       t0 = timer();
       for (i = 0; i < nk; i++) {

	  /* multiply by A'A */
          opb(nrow, ncol, &v[indexv], tmpv);
	  tmp1 = ddot(ncol, &v[indexv], 1, tmpv, 1);
	  daxpy(ncol, -tmp1, &v[indexv], 1, tmpv, 1);
	  tmp1 = sqrt(tmp1);
	  xnorm = sqrt(ddot(ncol, tmpv, 1, tmpv, 1));

	  /* multiply by A to get scaled left singular vectors */
	  opa(nrow, ncol, &v[indexv], &u[indexu]);
	  tmp2 = ONE / tmp1;
          dscal(nrow, tmp2, &u[indexu], 1);
	  res[i] = xnorm * tmp2;
	  sing[i] = tmp1;
          write (fp_out2, (char *)&u[indexu], size1);
          write (fp_out2, (char *)&v[indexv], size2);
	  indexu += nrow;
	  indexv += ncol;
       }
       exetime += timer() - t0;
    }
    /* write results to output file */
    fprintf(fp_out1,"\n");
    fprintf(fp_out1, " ... \n");
    fprintf(fp_out1, " ... HYBRID BLOCK LANCZOS [A^TA]\n");
    fprintf(fp_out1, " ... NO. OF EQUATIONS          =%10ld\n", nn);
    fprintf(fp_out1, " ... MAX. NO. OF ITERATIONS    =%10ld\n", maxit);
    fprintf(fp_out1, " ... NO. OF ITERATIONS TAKEN   =%10ld\n", iter);
    fprintf(fp_out1, " ... NO. OF TRIPLETS SOUGHT    =%10ld\n", nums);
    fprintf(fp_out1, " ... NO. OF TRIPLETS FOUND     =%10ld\n", nk);
    fprintf(fp_out1, " ... INITIAL BLOCKSIZE         =%10ld\n", nbold);
    fprintf(fp_out1, " ... FINAL   BLOCKSIZE         =%10ld\n", nb);
    fprintf(fp_out1, " ... MAXIMUM SUBSPACE BOUND    =%10ld\n", ncold);
    fprintf(fp_out1, " ... FINAL   SUBSPACE BOUND    =%10ld\n", nc);
    fprintf(fp_out1, " ... NO. MULTIPLICATIONS BY A  =%10ld\n", mxvcount);
    fprintf(fp_out1, " ... NO. MULT. BY TRANSPOSE(A) =%10ld\n", mtxvcount);
    fprintf(fp_out1, " ... TOTAL SPARSE MAT-VEC MULT.=%10ld\n", 
            mxvcount + mtxvcount);
    fprintf(fp_out1, " ... MEMORY NEEDED IN BYTES    =%10ld\n", memory);
    if (vectors) 
       fprintf(fp_out1, " ... WANT S-VECTORS ON OUTPUT?           T\n");
    else
       fprintf(fp_out1, " ... WANT S-VECTORS ON OUTPUT?           F\n");
    fprintf(fp_out1, " ... TOLERANCE                 =%10.2E\n\n", tol);
    fprintf(fp_out1, " %s\n", title);
    fprintf(fp_out1, "           %s\n", name);
    fprintf(fp_out1, " ... NO. OF DOC.  (ROWS OF A)   = %8ld\n", nrow);
    fprintf(fp_out1, " ... NO. OF TERMS (COLS OF A)   = %8ld\n", ncol);
    fprintf(fp_out1, " ... \n\n");
    fprintf(fp_out1, " ...... BLSVD EXECUTION TIME=%10.2E\n",exetime);
    fprintf(fp_out1, " ...... \n");
    fprintf(fp_out1, " ......         COMPUTED S-VALUES     (RES. NORMS)\n");
    fprintf(fp_out1, " ...... \n");
    for (i = 0; i < nk; i++)
           fprintf(fp_out1," ...... %3ld   %22.14E  (%11.2E)\n",i+1,sing[i],res[
i]);


    free(value);
    free(rowind);
    fclose(fp_in1);
    fclose(fp_in2);
    fclose(fp_out1);
    if (vectors) {
       free(tmpv);
       close(fp_out2);
    }
    exit(0);
}
#include <stdio.h>
#include <math.h>
#include <errno.h>
#include <fcntl.h>

extern double *tres, *y, *vv, *v0, *uvtmp, *pp, *qq;
extern double *alpha, *beta, *p, *q, *t, *z, *ztemp;
extern double **yp, **vvp, **uvtmpp, **ppp, **qqp, **ztempp;
extern long iconv, nn, iter;

#define   SEED      91827211
#define   MINBLKS   2
#define   TRANSP    1
#define   NTRANSP   0
#define   TRUE      1
#define   FALSE     0
#define   ZERO      0.0
#define   ONE       1.0

long   tql2(long, double *, double *, double **);
void   dgemv(long, long, long, double, double **, double *, double, double *);
void   dgemm2(long, long, long, long, long, double, double **, double **, 
              double, double **);
long   imin(long, long);
double random(long *);
long   validate(FILE *, long, long, long, long, long, long, double);
long   polong2(long *, long, long, double *, double *, double *, double);
void   block2(double **, double **, double **, long, long, long, long, long *);
void   orthg(long, long, long, double **, double **, double *);

/***********************************************************************
 *                                                                     *
 *                        blklan2()                                    *
 *                  Block Lanczos SVD Alogrithm                        *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   blklan2.c is designed to compute singular values and singular vectors
   of a large sparse matrix A.  The singular values (and singular vectors)
   of A are determined by the eigenvalues (and eigenvectors) of the matrix
   B, where B = A'A.

   The eigenvalues of B are the squares of the singular values of A, the
   eigenvectors of B correspond to the right singular vectors of A only.
   The left singular vectors of A are determined by 

		u = 1/sigma A*v,