svdpack.cc

矩阵奇异分解(svd)最新c++版本

   version, vectors are accumulated and B is assumed to be a square matrix.


   Parameters
   ---------

   (input)
   n           order of B
   s           linear array containing the diagonal elements of B
   e           linear array containing the off-diagonal elements of B

   (output)
   s           contains the singular values of the original matrix B
   up          2-dimensional array containing left singular vectors of B
   vp          2-dimensional array containing right singular vectors of B


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

   BLAS         dscal, dswap, drot, drotg 
   MISC         dmax

 ***********************************************************************/
void svdpack::qriter2(long n, double *s, double *e, double **up, double **vp)

{
   long negligible, iter, m, mm1, k, l, qrcase;
   double ztest, test, sl, g, t, smm1;
   double f, t1, cs, sn, scale, sm, el, emm1, b, c, shift;

   m     = n - 1;
   iter  = 0;
   while (m >= 0) {
      if (iter >= MAXIT) return;
      negligible = FALSE;

      /* this portion of the code inspects for negligible elements
       * in the s and e arrays.  On completion the variable qrcase
       * is set as follows:
       * qrcase = CASE1 if s[m] and e[l-1] are negligible and l < m
       * qrcase = CASE2 if s[l] is negligible and l < m
       * qrcase = CASE3 if s[l], ..., s[m] are not negligible (QR step),
       *                e[l-1] is negligible and l < m
       * qrcase = CONVERGENCE if e[m-1] is negligible */

      for (l = m - 1; l >= 0; l--) {
	 test = fabs(s[l]) + fabs(s[l + 1]);
	 ztest = test + fabs(e[l]);
	 if (ztest == test) {
	    e[l] = ZERO;
	    negligible = TRUE;
	 }
	 if (negligible) break;
      }

      if (l == m - 1) qrcase = CONVERGE;
      else {
         negligible = FALSE;
         for (k = m; k > l; k--) {
	    test = ZERO;
	    if (k != m) test += fabs(e[k]);
	    if (k != l + 1) test += fabs(e[k-1]);
	    ztest = test + fabs(s[k]);
	    if (ztest == test) {
	       s[k] = ZERO;
	       negligible = TRUE;
	    }
	    if (negligible) break;
	 }
	 if (k == l) qrcase = CASE3;
	 else if (k == m) qrcase = CASE1;
	 else {
	    qrcase = CASE2;
	    l = k;
	 }
      }
      l += 1;

      switch(qrcase) {

         /* deflate negligible s[m] */
	 case CASE1:    mm1 = m - 1;
			f = e[mm1];
			e[mm1] = ZERO;
			for (k = mm1; k >= l; k--) {
			   t1 = s[k];
			   drotg(&t1, &f, &cs, &sn);
			   s[k] = t1;
			   if (k != l) {
			      f = -sn * e[k - 1];
			      e[k - 1] *= cs;
			   }
			   drot(n, vp[k], vp[m], cs, sn);
                        }
			break;

         /* split at negligible s[l] */
	 case CASE2:    f = e[l - 1];
			e[l - 1] = ZERO;
			for (k = l; k <= m; k++) {
			   t1 = s[k];
			   drotg(&t1, &f, &cs, &sn);
			   s[k] = t1;
			   f = -sn * e[k];
			   e[k] *= cs;
			   drot(n, up[k], up[l - 1], cs, sn);
                        }
			break;

         /* perform one QR step */
	 case CASE3:    f = e[l - 1];
			/* calculate the shift */
			scale = dmax(fabs(s[m]), fabs(s[m - 1]));
			if (scale < fabs(e[m - 1])) scale = fabs(e[m - 1]);
			if (scale < fabs(s[l])) scale = fabs(s[l]);
			if (scale < fabs(e[l])) scale = fabs(e[l]);
			sm = s[m] / scale;
			smm1 = s[m - 1] / scale;
			emm1 = e[m - 1] / scale;
			sl = s[l] / scale;
			el = e[l] / scale;
			b = ((smm1 + sm) * (smm1 - sm) + emm1 * emm1) / 2.0;
			c = (sm * emm1);
			c *= c;
			shift = ZERO;
			if (b != ZERO || c !=ZERO) {
			   shift = sqrt(b * b + c);
			   if (b < ZERO) shift = -shift;
			   shift = c / (b + shift);
			}
			f = (sl + sm) * (sl - sm) + shift;
			g = sl * el;

			/* chase zeros */
			mm1 = m - 1;
			for (k = l; k <= mm1; k++) {
			   drotg(&f, &g, &cs, &sn);
			   if (k != l) e[k - 1] = f;
			   f = cs * s[k] + sn * e[k];
			   e[k] = cs * e[k] - sn * s[k];
			   g = sn * s[k + 1];
			   s[k + 1] = cs * s[k + 1];
			   drot(n, vp[k], vp[k + 1], cs, sn);
			   drotg(&f, &g, &cs, &sn);
			   s[k] = f;
			   f = cs * e[k] + sn * s[k + 1];
			   s[k + 1] = -sn * e[k] + cs * s[k + 1];
			   g = sn * e[k + 1];
			   e[k + 1] = cs * e[k + 1];
			   if (k < n - 1)
			      drot(n, up[k], up[k + 1], cs, sn);
			}
			e[mm1] = f;
			iter += 1;
			break;

         /* convergence */
	 case CONVERGE: if (s[l] < ZERO) {
			   /* make singular value positive */
	                   s[l] = -s[l];
			   dscal (n, -ONE, vp[l], 1); 
                        }
			/* order the singular value */
			while (l < n - 1) {
			   if (s[l] < s[l + 1]) {
			      t = s[l];
			      s[l] = s[l + 1];
			      s[l + 1] = t;
			      if (l < n - 1) 
				 dswap(n, vp[l], 1, vp[l + 1], 1);
			      if (l < n - 1) 
				 dswap(n, up[l], 1, up[l + 1], 1);
                              l += 1;
			   }
			   else break;
			}
			iter = 0;
			m -= 1;
			break;
      }
   }
}
/***************************************************************** 
 * applies a plane rotation;  assume a stride of 1 for dx and dy *
 * based on FORTRAN 77 routine from Linpack by J. Dongarra       *
 *****************************************************************/ 

void svdpack::drot(long n, double *dx, double *dy, double c, double s)

{
   long i;
   double temp;

   if (n <= 0) return;

   for (i = 0; i < n; i++) {
      temp = c * (*dx) + s * (*dy);
      *dy = c * (*dy) - s * (*dx);
      dy++;
      *dx++ = temp;
   }
   return;
}



/***************************************************************** 
 * constructs Givens plane rotation                              *
 * based on FORTRAN 77 routine from Linpack by J. Dongarra       *
 *****************************************************************/ 

void svdpack::drotg(double *da, double *db, double *c, double *s)

{
   double r, roe, scale, z, temp1, temp2;

   roe = *db;
   temp1 = fabs(*da);
   temp2 = fabs(*db);
   if (temp1 > temp2) roe = *da;
   scale = temp1 + temp2;

   if (scale != ZERO) {
      temp1 = *da / scale;
      temp2 = *db / scale;
      r = scale * sqrt(temp1 * temp1 + temp2 * temp2);
      r *= fsign(ONE, roe);
      *c = *da / r;
      *s = *db / r;
   }
   else {
      *c = ONE;
      *s = ZERO;
      r = ZERO;
   }
   z = *s;

   temp1 = fabs(*c);
   if (temp1 > ZERO && temp1 <= *s) z = ONE / *c;

   *da = r;
   *db = z;
   return;
}

/************************************************************** 
 * Function forms the dot product of two vectors.      	      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

double svdpack::ddot(long n,double *dx,long incx,double *dy,long incy)

{
   long i;
   double dot_product;

   if (n <= 0 || incx == 0 || incy == 0) return(0.0);
   dot_product = 0.0;
   if (incx == 1 && incy == 1) 
      for (i=0; i < n; i++) dot_product += (*dx++) * (*dy++);
   else {
      if (incx < 0) dx += (-n+1) * incx;
      if (incy < 0) dy += (-n+1) * incy;
      for (i=0; i < n; i++) {
         dot_product += (*dx) * (*dy);
         dx += incx;
         dy += incy;
      }
   }
   return(dot_product);
}
/************************************************************** 
 * Function scales a vector by a constant.     		      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

void svdpack::dscal(long n,double da,double *dx,long incx)

{
   long i;

   if (n <= 0 || incx == 0) return;
   if (incx < 0) dx += (-n+1) * incx;
   for (i=0; i < n; i++) {
      *dx *= da;
      dx += incx;
   }
   return;
}
/************************************************************** 
 * Constant times a vector plus a vector     		      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

void svdpack::daxpy (long n,double da,double *dx,long incx,double *dy,long incy)

{
   long i;

   if (n <= 0 || incx == 0 || incy == 0 || da == 0.0) return;
   if (incx == 1 && incy == 1) 
      for (i=0; i < n; i++) {
	 *dy += da * (*dx++);
	 dy++;
      }
   else {
      if (incx < 0) dx += (-n+1) * incx;
      if (incy < 0) dy += (-n+1) * incy;
      for (i=0; i < n; i++) {
         *dy += da * (*dx);
         dx += incx;
         dy += incy;
      }
   }
   return;
}

/***********************************************************************
 *                                                                     *
 *				mrandom()                               *
 *                        (double precision)                           *
 ***********************************************************************/
/***********************************************************************

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

   This is a translation of a Fortran-77 uniform random number
   generator.  The code is based  on  theory and suggestions  given in
   D. E. Knuth (1969),  vol  2.  The argument to the function should 
   be initialized to an arbitrary integer prior to the first call to 
   random.  The calling program should  not  alter  the value of the
   argument between subsequent calls to random.  Random returns values
   within the the interval (0,1).


   Arguments 
   ---------

   (input)
   iy	   an integer seed whose value must not be altered by the caller
	   between subsequent calls

   (output)
   random  a double precision random number between (0,1)

 ***********************************************************************/
double svdpack::mrandom(long *iy)

{
   static long m2 = 0;
   static long ia, ic, mic;
   static double halfm, s;

   /* If first entry, compute (max int) / 2 */
   if (!m2) {
      m2 = 1 << (8 * (int)sizeof(int) - 2); 
      halfm = m2;

      /* compute multiplier and increment for linear congruential 
       * method */
      ia = 8 * (long)(halfm * atan(1.0) / 8.0) + 5;
      ic = 2 * (long)(halfm * (0.5 - sqrt(3.0)/6.0)) + 1;
      mic = (m2-ic) + m2;

      /* s is the scale factor for converting to floating point */
      s = 0.5 / halfm;
   }

   /* compute next random number */
   *iy = *iy * ia;

   /* for computers which do not allow integer overflow on addition */
   if (*iy > mic) *iy = (*iy - m2) - m2;

   *iy = *iy + ic;

   /* for computers whose word length for addition is greater than
    * for multiplication */
   if (*iy / 2 > m2) *iy = (*iy - m2) - m2;
  
   /* for computers whose integer overflow affects the sign bit */
   if (*iy < 0) *iy = (*iy + m2) + m2;

   return((double)(*iy) * s);
}

/************************************************************** 
 * Function copies a vector x to a vector y	     	      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

void svdpack::dcopy(long n,double *dx,long incx,double *dy,long incy)

{
   long i;

   if (n <= 0 || incx == 0 || incy == 0) return;
   if (incx == 1 && incy == 1) 
      for (i=0; i < n; i++) *dy++ = *dx++;

   else {
      if (incx < 0) dx += (-n+1) * incx;
      if (incy < 0) dy += (-n+1) * incy;
      for (i=0; i < n; i++) {
         *dy = *dx;
         dx += incx;
         dy += incy;
      }
   }
   return;
}

/***********************************************************************
 *                                                                     *
 *                         dgemm2()                                    *
 *                                                                     *
 * A C-translation of the level 3 BLAS routine DGEMM by Dongarra,      *
 * Duff, du Croz, and Hammarling (see LAPACK Users' Guide).            *
 * In this version, the arrays which store the matrices used in this   *
 * matrix-matrix multiplication are accessed as two-dimensional arrays.*
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   dgemm2() performs one of the matrix-matrix operations

	      C := alpha * op(A) * op(B) + beta * C,

   where op(X) = X or op(X) = X', alpha and beta are scalars, and A, B
   and C are matrices, with op(A) an m by k matrix, op(B) a k by n
   matrix and C an m by n matrix.


   Parameters
   ----------

   (input)
   transa   TRANSP indicates op(A) = A' is to be used in the multiplication
	    NTRANSP indicates op(A) = A is to be used in the multiplication

   transb   TRANSP indicates op(B) = B' is to be used in the multiplication
	    NTRANSP indicates op(B) = B is to be used in the multiplication

   m        on entry, m specifies the number of rows of the matrix op(A)
	    and of the matrix C.  m must be at least zero.  Unchanged
	    upon exit.