bld2

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

- Introduction

        bls2: sparse svd via hybrid block Lanczos procedure for
              eigensystems of the form A'A.

	bls2.c is an ANSI-C code designed to compute singular values
	and singular vectors of a large sparse matrix A.  This is a     
	modified version of the block Lanczos algorithm first
        published by Golub, Luk, and Overton (ACM TOMS 7(2):149-169, 1981).
        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.

        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 correspond to the right singular vectors only.
        The left singular vectors of A are then determined by

                          u = 1/sigma A*v,

        where {u,sigma,v} is a singular triplet of A.


        This "hybrid" block Lanczos procedure consists of five phases:

        phase 1: Block Lanczos outer iteration to yield a symmetric
                 block tridiagonal matrix S which shares the same 
                 eigenvalues of the matrix B=A'A.  Total or
                 complete re-orthogonalization is used here.

        phase 2: Lanczos method for tridiagonalizing the S matrix from
                 phase 1 to yield the tridiagonal matrix T which preserves
                 the same eigenvalues.  Complete or total re-orthogonal-
                 ization is used for this Lanczos recursion.  A point
                 Lanczos method is used if a blocksize (nb) of 1 is encountered.
                     
        phase 3: Apply an appropriate QL iteration to diagonalize T and 
                 hence produce approximate singular values (array alpha)
                 of the original matrix A.

        phase 4: Convergence test using a user-supplied residual tolerance
                 (tol).

        phase 5: Iteration restart with orthogonal projection with respect
                 to any (all) converged eigenvectors of B=A'A.
	
- Calling sequence

	The calling sequence for procedure blklan2 is

        long   blklan2(FILE *fp, long nnzero, long m, long n, long ik, 
               double *v, double *eig, long ic, long ib, 
               double tol, double *res, long maxit, long *iko, long *ico,
               long *ibo, long *memory)

	The user specifies as part of the parameter list:
	 
        fp              ... a pointer to output file {FILE *}.
	nnzero          ... number of nonzeros in matrix A {long}.
        m               ... row dimension of the sparse matrix A whose SVD
                            is sought {long}.
	n		... column dimension of the sparse matrix a whose SVD
                            is sought {long}.
	ik     	        ... number of singular triplets desired {long}.
	ib     	        ... initial block size for outer iteration {long}.
	ic     	        ... upper bound for dimension of Krylov subspace 
                            generated via outer iteration {long}.  ic is the
                            maximum dimension for the block upper bidiagonal
                            matrix S generated in phase 1 above.
	tol        	... user-specified tolerance for approximate singular
                            triplets {double}.
	maxit           ... maximum number of outer iterations allowed
                            {long}.
	 
	blklan2 returns via its parameter list the following items:
	 
        ik0             ... number of singular triplets approximated {long}.
	ic0    	        ... last bound for dimension of Krylov subspace 
                            used within outer iteration {long}.
	ib0    	        ... final block size used in outer iteration {long}.
	eig             ... one-dimensional array containing the ik0
                            approximate singular values {double}.
	v	        ... two-dimensional array containing the ik0
                            approximate right singular vectors corresponding
                            to the approximate singular values in array
                            eig {double}.
	res             ... one-dimensional array containing the ik0
                            residuals of the approximate eigenpairs of
                            B=A'A {double}.
	memory          ... memory needed in bytes {long}.

- User-supplied routines

        For bls2.c, the user must specify multiplication by matrices
        A and B separately.  Subroutines opb and opm are used when
        multiplying B=A'A by a single vector and a block of vectors,
        respectively.  Subroutine opa is used for multiplying a
        single vector by the matrix A.

        The specification of opa should look something like

           void opa(long m, long n, double *x, double *y)

        so that opa takes an n by 1 vector x and returns the m by 1
        vector y = A*x, where A is m by n (m >> n).

        The specification of opb should look something like

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

        so that opb takes an n by 1 vector x and returns the n by 1
        vector y = A'A*x, where A is m by n (m >> n).

        The specification of opm should look something like

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

        so that opm takes an n by nc block of vectors X and returns 
        the n by nc block of vectors Y = A'A*X, where A is m by n (m >> n).

        This version of bls2.c is designed to approximate the ik-largest
        singular triplets of A.  Users interested in the ik-smallest
        singular triplets need only sort the alpha array in increasing
        (as opposed to the default ascending order) following the line

                  if (tql2(nn, alpha, beta, qqp)) break;

        in phase 3 of bls2.c.  The columns of the two-dimensional array
        qqp which are used to obtain approximate eigenpairs of the matrix
        A'A correspond to the elements of alpha (in ascending order)
        which are approximate eigenvalues of A'A (and hence squares of the
        singular values of the matrix A).  The columns of array qqp would
        have to be reordered to reflect a one-to-one correspondence
        with the newly sorted elements of alpha.

- Information 

        Please address all questions, comments, or corrections to:

        M. W. Berry
        Department of Computer Science
        University of Tennessee
        107 Ayres Hall
        Knoxville, TN  37996-1301
        email: [e-mail]
        phone: [phone]

-File descriptions

        bls2.c requires the include file bls2.h for compilation.
        The local parameters defined in bls2.h are:

         k        remaining # of desired triplets (done when = 0)
         k0       count of triplets found in current iteration
         nb       current block size
         nc       size of current subspace
         ns       number of blocks in current iteration

        The input and output files associated with bls2.c are
        listed below.

             Code           Input         Output
            ------      ------------    ---------
            bls2.c      blp2, matrix    blo2,blv2


       The binary output file blv2 (which contains the 
       approximate right singular vectors followed by the ap-
       proximate left singular vectors) will be created by
       bls2.c if it does not already exist.  If you are
       running on a Unix-based workstation you should uncomment
       the line

                 /*   #define  UNIX_CREAT */

       in the declarations prior to main() in bls2.c.
       UNIX_CREAT specifies the use of the UNIX "creat" system 
       routine with the permissions defined by the PERMS constant

                  #define PERMS 0664

       You may adjust PERMS for the desired permissions on the
       blv2 file (default is Read/Write for user and group,
       and Read for others).  Subsequent runs will be able to
       open and overwrite these files with the default permissions.

       bls2.c obtains its parameters specifying the
       sparse SVD problem to be solved from the input file
       blp2. This parameter file contains the single line

	 <name>   maxit   nc   nb   nums   tol   vtf

       where 

        <name>     is the name of the data set containing the nonzeros of A.
        maxit      is an integer specifying maximum number of (outer)
                   block Lanczos iterations allowed.
        nc         is an integer specifying the upper bound for the
                   Krylov subspace generated via the outer iteration.
        nb         is an integer specifying the initial block size for 
                   the outer iteration.
        nums       is an integer specifying the number of singular triplets
                   desired.
        tol        is a double specifying the residual tolerance for
                   approximated singular triplets.
        vtf        contains the string TRUE or FALSE to indicate when 
                   singular triplets are needed (TRUE) and when only 
                   singular values are needed (FALSE);  If vtf is TRUE,
                   the unformatted output file blv2 will contain the 
                   approximate singular vectors written in the order 

                      u[1], v[1], u[2], v[2], ..., u[ik0], v[ik0]. 

                   Here u[i] and v[i] denote the left and right
                   singular vectors, respectively, corresponding to the 
                   i-th approximate singular value, sing[i].

- Sparse matrix format

        bls2.c is designed to read input matrices that are stored
        in the Harwell-Boeing sparse matrix format.  The nonzeros
        of such matrices are stored in a compressed column-oriented
        format.  The row indices and corresponding nonzero values
        are stored by columns with a column start index array
        whose entries contain pointers to the nonzero starting each
        column.  bls2.c reads the sparse matrix data from the input
        file called "matrix".

        Each input file "matrix" should begin with a four-line header
        record followed by three more records containing, in order, 
        the column-start pointers, the row indices, and the nonzero
        numerical values.

        The first line of the header consists of a 72-character title
        and an 8-character key by which the matrices are referenced.
        The second line can be used for comments or to indicate record
        length for each index or value array.  Although this line is 
        generally ignored, A CHARACTER MUST BE PLACED ON THAT LINE.
        The third line contains a three-character string denoting the
        matrix type and the three integers specifying the number of rows,
        columns, and nonzeros.  The fourth line which usually contains
        input format for Fortran-77 I/O is ignored by our ANSI-C code.
        The exact format is

		"%72c %*s %*s %*s %d %d %d %*d"

	for the first three lines of the header,

		line 1      <title>         <key>
		 	(col.  1 - 72) (col. 73 - 80)

		line 2   <string>

		line 3   <matrix type> nrow ncol nnzero 

	and 

		"%*s %*s %*s %*s"

	for the last line of the header.

		line 4   <string1> <string2> <string3> <string4>

        Even though only the title and the integers specifying the
        number of rows, columns, and nonzero elements are read, other
        strings of input must be present in indicated positions.
        Otherwise, the format of the "fscanf" statements must be 
        changed accordingly.