bls2.c

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

extern double **vvp, **uvtmpp, **yp, **ppp, **qqp;
extern double *p, *q, *t, *z, *alpha, *beta;
extern long iconv, nn, iter;

void   dsbmv(long, long, double, double **, double *, double, double *);
void   dgemm2(long, long, long, long, long, double, double **, double **,
            double, double **);
void   orthg(long, long, long, double **, double **, double *);
void   formbigs(long, long, double **, double **, double **);
double ddot(long, double *,long, double *, long);
void   daxpy(long, double, double *,long, double *, long);
double random(long *);
double enorm(long, double *);
void   opm(long, long, long, double **, double **);

/***********************************************************************
 *                                                                     *
 *                        block2()                                     *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   This function implements the first two phases of the hybrid block
   Lanczos procedure.  In the first phase, which is also known as the
   block Lanczos outer iteration, a symmetric block tridiagonal matrix S
   is formed.  The eigenvalues of the matrix S approximate those 
   of matrix B, where B = A'A and A is the original sparse matrix.  
   Total (or complete) re-orthogonalization is used.

   In the second phase, single-vector Lanczos tridiagonalization is used
   to reduce (preserving eigenvalues of S) the block matrix S to a 
   symmetric tridiagonal matrix T. 


   Arguments
   ---------

   (input)
   w, wp     work space
   sp        diagonal blocks (symmetric submatrices) of the
		symmetric block tridiagonal matrix S 
   rp        super-diagonal blocks (upper-triangular submatrices)
		of the symmetric block tridiagonal matrix S 
   bigsp     symmetric block tridiagonal matrix S
   m         row dimension of sparse matrix A
   n         column dimension of sparse matrix A
   nb        current block size
   ns        number of blocks in current iteration
   irand     seed for random number generator
   

   (output - globally defined)
   alpha     diagonal elements of symmetric tridiagonal matrix T (reduced 
                from matrix S)
   beta      off-diagonal elements of symmetric tridiagonal matrix T 
		(reduced from matrix S)
   tres      residuals of approximate eigenvalues determined from
	     a previous set of block Lanczos outer iterations.
   ppp       array of eigenvectors of S


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

   Defined and documented in bls2.h


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

   BLAS         ddot, daxpy, enorm, dgemm2, orthg, dsbmv
   USER         opm
   MISC         random
   BLS2         formbigs

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

void block2(double **sp, double **rp, double **bigsp, long m, long n, 
            long nb, long ns, long *irand)

{
   long jinc, nk, nj, i, j, k, blkptr;
   double *ptr, dum, pnm;

   for (i = 0; i < nn; i++)
      for (j = 0; j < nb; j++) sp[i][j] = ZERO;
   /* ns (number of blocks) is assumed to be at least 2 */
   nk = nn - nb;
   for (i = 0; i < nk; i++)
      for (j = 0; j < nb; j++) rp[i][j] = ZERO;

   opm(m, n, nb, vvp, yp);

   dgemm2(NTRANSP, TRANSP, nb, nb, n, ONE, vvp, yp, ZERO, sp);

   blkptr = 0;
   for (j = 1; j < ns; j++) {

      dgemm2(TRANSP, NTRANSP, nb, n, nb, -ONE, &sp[blkptr],
            &vvp[blkptr], ONE, yp);
      if (j > 1)
         dgemm2(NTRANSP, NTRANSP, nb, n, nb, -ONE, &rp[blkptr - nb],
               &vvp[blkptr - nb], ONE, yp);

      if (j == 1 && iter > 1)
         for (i = 0; i < nb; i++) tres[i] = enorm(n, yp[i]);

      if (iconv) {
         nk = nb * j;
         nj = nk + iconv;

         ptr = vv;
         for (k = 0; k < nk; k++) 
            for (i = 0; i < n; i++) uvtmpp[k][i] = *ptr++; 
         ptr = v0;
         for (k = nk; k < nj; k++) 
            for (i = 0; i < n; i++) uvtmpp[k][i] = *ptr++; 
      }
      else {
         nj = nb * j;
         ptr = vv;
         for (k = 0; k < nj; k++) 
            for (i = 0; i < n; i++) uvtmpp[k][i] = *ptr++;
      }
      ptr = y;
      for (k = 0; k < nb; k++) 
         for (i = 0; i < n; i++) {
	    uvtmpp[nj + k][i] = *ptr;
	    *ptr++ = ZERO;
         }
      orthg(nb, nj, n, yp, uvtmpp, ztemp);

      for (k = 0; k < nb; k++) 
         for (i = k; i < nb; i++) rp[blkptr + k][i] = yp[i][k];

      jinc = blkptr + nb;
      ptr = vvp[jinc];
      for (k = nj; k < nj + nb; k++) 
         for (i = 0; i < n; i++) *ptr++ = uvtmpp[k][i]; 
  
      opm(m, n, nb, &vvp[jinc], yp);
      dgemm2(NTRANSP, TRANSP, nb, nb, n, ONE, &vvp[jinc], yp, ZERO, &sp[jinc]);
      blkptr += nb;
   }
   formbigs(nn, nb, rp, sp, bigsp);

   for (i = 0; i < nn; i++) p[i] = random(irand);
   pnm = enorm(nn,p);

   ptr = pp;
   for (i = 0; i < nn; i++) {
      p[i] /= pnm;
      *ptr++ = p[i];
      z[i] = ZERO;
   }
   for (j = 0; j < nn; j++) {
      dsbmv(nn, nb, ONE, bigsp, p, ONE, z);
      alpha[j] = ddot(nn, p, 1, z, 1);
      if (j == nn - 1) break;

      /* compute Z[j] := Z[j] - alpha[j] * P */
      daxpy(nn, -alpha[j], p, 1, z, 1);

      /* orthogonalize Z w.r.t. previous PP's */
      for (k = 0; k <= j; k++)
         for (i = 0; i < nn; i++) uvtmpp[k][i] = ppp[k][i];
      if (j) {
         ptr = uvtmpp[j + 1];
         for (i = 0; i < nn; i++) *ptr++ = z[i];
	 orthg(1, j+1, nn, yp, uvtmpp, ztemp);
	 ptr = uvtmpp[j + 1];
         for (i = 0; i < nn; i++) z[i] = *ptr++;
      }
      else {
         dum = -ddot(nn, uvtmp, 1, z, 1);
	 daxpy(nn, dum, uvtmp, 1, z, 1);
      }
      beta[j] = enorm(nn,z);
      if (beta[j] != ZERO) {

	 /* compute P[j+1] := Z[j] / beta[j] */
         ptr = ppp[j + 1];
         for (i = 0; i < nn; i++) {
	    t[i] = p[i];
	    p[i] = z[i] /beta[j];
	    *ptr++ = p[i];
	    z[i] = -beta[j] * t[i];
	 }
      }
      else return;
   }
   return;
}
#define       ZERO       0.0
#define       ONE        1.0
#define       TRANSP     1
#define       NTRANSP    0
#define       CONTINUE   1
#define       DONE       0

extern double *v0, *ztemp, *uvtmp, *vv, *alpha, *beta;
extern double *p, *q, *t, *z;
extern double **uvtmpp, **vvp, **yp;
extern long iconv, nn, iter;

double ddot(long, double *,long, double *, long);
double enorm(long, double *);
void   daxpy(long, double, double *,long, double *, long);
void   opb(long, long, double *, double *);
void   orthg(long, long, long, double **, double **, double *);

/***********************************************************************
 *                                                                     *
 *                        polong2()                                     *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   This function is a single-vector Lanczos tridiagonalization procedure
   (degenerate case of block size = 1 for block2.c) which is used when
   normal deflation reduces the current block size to 1.

   The function returns DONE when the number of remaining triplets to be 
   approximated is less than or equal to zero.  Otherwise, it returns
   status CONTINUE.
   

   Arguments
   ---------

   (input)
   k        current # of desired triplets
   m        row dimension of the sparse matrix A whose SVD is sought
   n        column dimension of the sparse matrix A whose SVD is sought
   tol     user-specified tolerance for approximate singular triplets
   wp      work space

   (output)
   sing    linear array containing the iko approximate singular values
   res     linear array containing the iko residuals of the approximate 
	   singular triplets
   alpha   diagonal elements of symmetric tridiagonal matrix from the 
	   inner recursion
   beta    off-diagonal elements of symmetric tridiagonal matrix from the 
	   inner recursion

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

   Defined and documented in bls2.h


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

   BLAS         ddot, enorm, daxpy, orthg
   USER         opb

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

long polong2(long *k, long m, long n, 
            double *v, double *res, double *eig, double tol)

{
   double *ptr, dum, znm;
   long i, j, jj, convpj;

   for (j = 0; j < n; j++) p[j] = vv[j];
   for (j = 0; j < n; j++) z[j] = ZERO;

   for (j = 0; j < nn; j++) {
      opb(m, n, p, q);
      daxpy(n, ONE, q, 1, z, 1);
      alpha[j] = ddot(n, p, 1, z, 1);
      if (j == nn - 1) return(CONTINUE);

      /* compute Z[j] := Z[j] - alpha[j] * P */
      if (alpha[j] != ZERO) {
	 daxpy(n, -alpha[j], p, 1, z, 1);
	 if (!j) {
	    if ((znm = enorm(n, z)) <= tol) {
	       ptr = &v[iconv * n];
               for (i = 0; i < n; i++) *ptr++ = p[i];
	       eig[iconv] = alpha[0];
	       res[iconv] = znm;
               *k -= 1;
	       iter -= 1;
	       return(DONE);
            }
	 }

	 /* orthogonalize Z w.r.t. converged right S-vectors and previous 
	    VV's */
         convpj = iconv + j;
	 ptr = v0;
	 for (jj = 0; jj < iconv; jj++)
	    for (i = 0; i < n; i++)
	       uvtmpp[jj][i] = *ptr++;
	 ptr = vv;
	 for (jj = iconv; jj <= convpj; jj++)
	    for (i = 0; i < n; i++)
	       uvtmpp[jj][i] = *ptr++;
         if (convpj) {
	    ptr = uvtmpp[convpj + 1];
	    for (i = 0; i < n; i++) *ptr++ = z[i];
	    orthg(1, convpj + 1, n, yp, uvtmpp, ztemp);
	    ptr = uvtmpp[convpj + 1];
	    for (i = 0; i < n; i++) z[i] = *ptr++;
	 }
	 else {
	    dum = -ddot(n, uvtmp, 1, z, 1);
	    daxpy(n, dum, uvtmp, 1, z, 1);
	 }

	 /* compute beta[j] */
	 beta[j] = enorm(n,z);
	 if (beta[j] != ZERO) {

            /* compute P[j+1] := Z[j] / beta[j] */
	    ptr = vvp[j + 1];
	    for (i = 0; i < n; i++) {
	       t[i] = p[i];
	       p[i] = z[i] / beta[j];
	       *ptr++ = p[i];
	       z[i] = -beta[j] * t[i];
	    }
         }
	 else return(CONTINUE);
      }
      else return(CONTINUE);
   }
}
#include <stdio.h>
#include <fcntl.h>

#define  NCMAX  5200
#define  NZMAX  800000
#define  ZERO   0.0

/***********************************************************************
 *								       *
 *		              validate()			       *
 *								       *
 ***********************************************************************/
/***********************************************************************

   Description
   -----------
   Function validates input parameters and returns error code (long)  

   Arguments 
   ---------
  (input)
  fp_out1   polonger to output file
  nnzero    number of nonzero elements in input sparse matrix A
  nrow      row dimension of A
  ncol      column dimension of A
  maxsubsp  maximum dimension of Krylov subspace allowed
  blksize   initial block size
  nums      number of singular values desired
  tol       user-specified tolerance for approximate singular triplets

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

long validate(FILE *fp_out1, long nnzero, long nrow, long ncol,
	      long maxsubsp, long blksize, long nums, double tol)

{
   long error_index;
   char *error[8];

   error_index = 0;
   error[1] = " ***** SORRY, YOUR MATRIX IS TOO BIG *****";
   error[2] = " ***** NCOL MUST NOT BE GREATER THAN NROW *****";
   error[3] = " ***** TOLERANCE IS INVALID *****";
   error[4] = " ***** MAXIMUM SUBSPACE DIMENSION IS INVALID *****";
   error[5] = " ***** INITIAL BLOCK SIZE MUST BE GREATER THAN 1 *****";
   error[6] = " ***** NUMBER OF SINGULAR VALUES DESIRED IS INVALID *****";
   error[7] = " ***** INIT BLK SIZE MUST BE LESS THAN NO. OF S-VALUES DESIRED *****";

   if (ncol > NCMAX || nnzero > NZMAX)           error_index = 1;
   else if (ncol > nrow)                         error_index = 2;
   else if (tol < ZERO)                          error_index = 3;
   else if (maxsubsp > ncol || maxsubsp <= 0)    error_index = 4;
   else if (blksize <= 1 || blksize > maxsubsp)  error_index = 5;
   else if (nums > maxsubsp || nums <= 0)        error_index = 6;
   else if (blksize > nums)                      error_index = 7;

   if (error_index) fprintf(fp_out1, "%s\n", error[error_index]);
   return(error_index);
}
extern long mxvcount, mtxvcount;
extern long *pointr, *rowind;
extern double *value,*ztemp;

#define         ZERO         0.0

/**************************************************************
 *                                                            *
 * multiplication of matrix B by vector x, where B = A'A,     *
 * and A is nrow by ncol (nrow >> ncol) and is stored using   *
 * the Harwell-Boeing compressed column sparse matrix format. *
 * Hence, B is of order n = ncol.  y stores product vector.   *
 *                                                            *