bls1.c

求矩阵奇异分解svd算法

   }
   /* compute alpha[0] */
   alpha[0] = enorm(m, t);
   for (j = 0; j < nn; j++) {
      if (alpha[j] != ZERO) {
	 
	 /* compute Q[j] := T[j] / alpha[j] */
	 ptr = uup[j];
         for (i = 0; i < m; i++) {
	    q[i] = t[i] / alpha[j];
	    *ptr++ = q[i];
	 }
	 if (j == nn - 1) return(CONTINUE);

	 /* compute Z[j] := A^T * Q[J] - alpha[j] * P[j] */
	 opat(n, q, z);
	 daxpy(n, -alpha[j], p, 1, z, 1);
	 if (!j) {
	    if ((znm = enorm(n, z)) <= tol) {
	       ptr = &u0[iconv * m];
               for (i = 0; i < m; i++) *ptr++ = q[i];
	       ptr = &v0[iconv * n];
               for (i = 0; i < n; i++) *ptr++ = p[i];
	       sing[iconv] = alpha[0];
	       res[iconv] = znm;
	       iconv += 1;
               *k -= 1;
	       if (!*k) {
		  iter -= 1;
		  return(DONE);
	       }
            }
	 }
         convpj = iconv + j;
	 ptr = v0;

	 /* orthogonalize Z w.r.t. converged right S-vectors and 
	  * previous VV's */
	 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, wp, uvtmpp, temp);
	    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);
	 }
	 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++) {
	       p[i] = z[i] /beta[j];
	       *ptr++ = p[i];
	    }
	    /* compute T[j+1] := A * P[j+1] - beta[j] * Q[j] */
	    opa(m, n, p, t);
	    daxpy(m, -beta[j], q, 1, t, 1);

	    /* orthogonalize T w.r.t. converged left S-vectors and 
	     * previous UU's */

            convpj = iconv + j;
	    ptr = u0;
	    for (jj = 0; jj < iconv; jj++)
	       for (i = 0; i < m; i++)
	          uvtmpp[jj][i] = *ptr++;
	    ptr = uu;
	    for (jj = iconv; jj <= convpj; jj++)
	       for (i = 0; i < m; i++)
	          uvtmpp[jj][i] = *ptr++;
            if (convpj) {
	       ptr = uvtmpp[convpj + 1];
	       for (i = 0; i < m; i++) *ptr++ = t[i];
	       orthg(1, convpj + 1, m, wp, uvtmpp, temp);
	       ptr = uvtmpp[convpj + 1];
	       for (i = 0; i < m; i++) t[i] = *ptr++;
	    }
	    else {
	       dum = -ddot(m, uvtmp, 1, t, 1);
	       daxpy(m, dum, uvtmp, 1, t, 1);
	    }
	    alpha[j + 1] = enorm(m, t);
         }
	 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   pointer 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);
}
#define       ZERO    0.0
extern long mxvcount;
extern long *pointr, *rowind;
extern double *value;

/**************************************************************
 *                                                            *
 * multiplication of matrix A by vector x, where A is m by n  *
 * (m >> n) and is stored using the Harwell-Boeing compressed *
 * column sparse matrix format.  y stores product vector.     *
 *                                                            *
 **************************************************************/

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

{
   long end,i,j;
   
   mxvcount += 1;
   for (i = 0; i < m; i++) y[i] = ZERO;

   for (i = 0; i < n; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)
	 y[rowind[j]] += value[j] * x[i]; 
   }
   return;
}
#define       ZERO    0.0
extern long mtxvcount;
extern long *pointr, *rowind;
extern double *value;

/**************************************************************
 *                                                            *
 * multiplication of an n by m matrix A' by a vector X, store *
 * in Y.                                                      *
 *                                                            *
 **************************************************************/

void opat(long n, double *x, double *y)
{
   long end,i,j;
   
   mtxvcount += 1;
   for (i = 0; i < n; i++) y[i] = ZERO;

   for (i = 0; i < n; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)
	 y[i] += value[j] * x[rowind[j]]; 
   }
   return;
}
double fsign(double a,double b)
/************************************************************** 
 * returns |a| if b is positive; else fsign returns -|a|      *
 **************************************************************/ 
{

   if ((a>=0.0 && b>=0.0) || (a<0.0 && b<0.0))return(a);
   if  ((a<0.0 && b>=0.0) || (a>=0.0 && b<0.0))return(-a);
}

double dmax(double a, double b)
/************************************************************** 
 * returns the larger of two double precision numbers         *
 **************************************************************/ 
{

   if (a > b) return(a);
   else return(b);
}

double dmin(double a, double b)
/************************************************************** 
 * returns the smaller of two double precision numbers        *
 **************************************************************/ 
{

   if (a < b) return(a);
   else return(b);
}

long imin(long a, long b)
/************************************************************** 
 * returns the smaller of two integers                        *
 **************************************************************/ 
{

   if (a < b) return(a);
   else return(b);
}

long imax(long a,long b)
/************************************************************** 
 * returns the larger of two integers                         *
 **************************************************************/ 
{

   if (a > b) return(a);
   else return(b);
}
#include <math.h>
#include <stdio.h>

#define		TRUE	 1
#define		FALSE  	 0
#define		TRANSP   1
#define		NTRANSP  0
#define		ZERO	 0.0
#define		ONE	 1.0
#define		CONST    100.0

void   dgemv(long, long, long, double, double **, double *, double, double *);
double ddot(long, double *,long, double *, long);
void   dscal(long, double, double *,long);
void   daxpy(long, double, double *,long, double *, long);
void   dcopy(long, double *, long, double *, long);

/***********************************************************************
 *                                                                     *
 *                        orthg()                                      *
 *         Gram-Schmidt orthogonalization procedure                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   The p by n matrix Z stored row-wise in rows f to (f+p-1) of
   array X is reorthogonalized w.r.t. the first f rows of array X.
   The resulting matrix Z is then factored into the product of a
   p by n orthonormal matrix (stored over matrix Z) and a p by p
   upper-triangular matrix (stored in the first p rows and columns 
   of array B).  (Based on orthog from Rutishauser) 


   Parameters
   ----------

   (input)
   p           number of consecutive vectors of array x (stored row-wise)
	       to be orthogonalized
   f           number of rows against which the next p rows are to be
	       orthogonalized
   n           column dimension of x
   x           2-dimensional array whose p rows are to be orthogonalized
	       against its first f rows
   temp        work array


   (output)
   x           output matrix whose f+p rows are orthonormalized
   b           p by p upper-triangular matrix


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

   BLAS         dgemv, ddot, dscal, daxpy, dcopy

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

void orthg(long p, long f, long n, double **b, double **x, double *temp)

{
   long fp, k, km1;
   long orig, small;
   double t, s;

   if (!p) return;
   if (f == 0 && p > n) {
      fprintf(stderr,"%s\n",
         "*** ON ENTRY TO ORTHG, MATRIX TO BE ORTHONORMALIZED IS SINGULAR");
      exit(-1);
   }
   fp = f + p;

   for (k = f; k < fp; k++) {
      km1 = k - 1;
      orig = TRUE;

      while(TRUE) {
         t = ZERO;

	 if (km1 >= 0) {
	    if (km1 > 0) {
	       dgemv(NTRANSP, k, n, ONE, x, x[k], ZERO, temp);
	       t += ddot(k, temp, 1, temp, 1);
	    }

	    else {
	       temp[0] = ddot(n, x[0], 1, x[k], 1);
	       t += temp[0] * temp[0];
	    }

	    if (orig && km1 >= f) 
               dcopy(k - f, &temp[f], 1, &b[k - f][0], 1); 

            if (km1 > 0) 
	       dgemv(TRANSP, k, n, -ONE, x, temp, ONE, &x[k][0]);
            else
	       daxpy(n, -temp[0], x[0], 1, x[k], 1);
         }

	 if (km1 < 0 || p != 1) {
	    s = ddot(n, x[k], 1, x[k], 1);
	    t += s;
	    if (s > t/CONST) {
	       small = FALSE;
	       s = sqrt(s);
               b[k - f][k - f] = s;
	       if (s != ZERO) s = ONE/s;
	       dscal(n, s, x[k], 1);
	    }
	    else {
	       small = TRUE;
	       orig  = FALSE;
	    }
	 }
	 if (small == FALSE || p == 1) break;
      }
   }
}
#include <stdio.h>

#define		TRANSP   1
#define		NTRANSP  0
#define		ZERO	 0.0
#define		ONE	 1.0

/***********************************************************************
 *                                                                     *
 *                         dgemv()                                     *
 * A C-translation of the level 2 BLAS routine DGEMV by Dongarra,      *
 * du Croz, and Hammarling, and Hanson (see LAPACK Users' Guide).      *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   dgemv() performs one of the matrix-vector operations

   y := alpha * A * x + beta * y  or  y := alpha * A' * x + beta * y

   where alpha and beta are scalars, X, Y are vectors and A is an
   m by n matrix.

void dgemv(long transa, long m, long n, 
           double alpha, double **a, double *x, double beta, double *y)

   Parameters