svdpack.cc

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

      /* form Y := alpha * A * X + Y */
      case NTRANSP:  for(i = 0; i < m; i++) {
                        ptrtemp = *a++;
		        temp = ZERO;
		        for(j = 0; j < n; j++) 
			   temp += *ptrtemp++ * x[j];
			y[i] += alpha * temp;
		     }
		     break;
		     
      /* form Y := alpha * A' * X + Y */
      case TRANSP:   for(i = 0; i < m; i++) { 
                        ptrtemp = *a++;
			if (x[i] != ZERO) {
			   temp = alpha * x[i];
			   for(j = 0; j < n; j++)
			      y[j] += temp * (*ptrtemp++);
			}
		     }
		     break;
   }
}



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

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

   dgemm() 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.

   Note that the arrays storing matrices B and C are linear arrays while
   the array of A is two-dimensional.


   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.

   n        on entry, n specifies the number of columns of the matrix op(B)
	    and of the matrix C.  n must be at least zero.  Unchanged
	    upon exit.

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

   alpha    a scalar multiplier

   a        matrix A as a 2-dimensional array.  When transa = NTRANSP, the
	    first k columns of the first m rows must contain the matrix A.
	    Otherwise, the first m columns of the first k rows must contain
	    the matrix A.

   b        matrix B as a linear array.  The leading (k * n) elements of
	    b must contain the matrix B.

   beta     a scalar multiplier.  When beta is supplied as zero then C
	    need not be set on input.

   c        matrix C as a linear array.  Before entry, the leading (m * n)
	    elements of c must contain the matrix C except when beta = 0.
	    In this case, c need not be set on entry.
	    On exit, c is overwritten by the (m * n) elements of matrix
	    (alpha * op(A) * op(B) + beta * C).

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

void svdpack::dgemm(long transa, long transb, long m, long n, long k, 
           double alpha, double **a, double *b, double beta, double *c)

{
   long info;
   long i, j, l, nrowa, ncola, nrowb, ncolb, nc;
   double temp, *atemp, *btemp1, *ptrtemp, *ctemp;

   info = 0;
   if      ( transa != TRANSP && transa != NTRANSP ) info = 1;
   else if ( transb != TRANSP && transb != NTRANSP ) info = 2;
   else if ( m < 0 ) 				     info = 3;
   else if ( n < 0 )				     info = 4;
   else if ( k < 0 )        			     info = 5;

   if (info) {
      fprintf(stderr, "%s %1ld %s\n",
      "*** ON ENTRY TO DGEMM, PARAMETER NUMBER",info,"HAD AN ILLEGAL VALUE");
      return; //exit(info);
   }

   if (transa) {
      nrowa = k;
      ncola = m;
   }
   else { 
      nrowa = m;
      ncola = k;
   }
   if (transb) {
      nrowb = n;
      ncolb = k;
   }
   else {
      nrowb = k;
      ncolb = n;
   }
   nc = m * n;

   if (!m || !n || ((alpha == ZERO || !k) && beta == ONE))
      return;

   ctemp = c; 
   if (alpha == ZERO) {
      if (beta == ZERO)
         for (i = 0; i < nc; i++) *ctemp++ = ZERO;
      else if (beta != ONE)
         for (i = 0; i < nc; i++) *ctemp++ *= beta;
      return;
   }

   if (beta == ZERO) 
      for (i = 0; i < nc; i++) *ctemp++ = ZERO;
   else if (beta != ONE) 
      for (i = 0; i < nc; i++) *ctemp++ *= beta;

   if (!transb) { 

      switch(transa) {

	 /* form C := alpha * A * B + beta * C */
	 case NTRANSP:  ptrtemp = c;
		        for(l = 0; l < nrowa; l++) {
                           atemp = *a++;
      	       		   btemp1 = b;
      		  	   for(j = 0; j < ncola; j++) {
	 	     	      temp = *atemp * alpha;
	             	      ctemp = ptrtemp;
         	     	      for(i = 0; i < ncolb; i++) 
	    			 (*ctemp++) += temp * (*btemp1++);
	 	     	      atemp++;
      		  	   }
   	       		   ptrtemp = ctemp;
   	    	        }
			break;

	 /* form C := alpha * A' * B + beta * C */
	 case TRANSP:   ptrtemp = b;
	                for(l = 0; l < nrowa; l++) {
                           atemp = *a++;
      	       		   ctemp = c;
      		  	   for(j = 0; j < ncola; j++) {
	 	     	      temp = *atemp * alpha;
	             	      btemp1 = ptrtemp;
         	     	      for(i = 0; i < ncolb; i++) 
	    			 (*ctemp++) += temp * (*btemp1++);
	 	     	      atemp++;
      		  	   }
   	       		   ptrtemp = btemp1;
   	    		}
			break;
      }
   }
   else { 
      ctemp = c;

      switch(transa) {

	 /* form C := alpha * A * B' + beta * C */
	 case NTRANSP: for(l = 0; l < nrowa; l++) {
      	       		   btemp1 = b;
      		  	   for(j = 0; j < nrowb; j++) {
	 	     	      atemp = *a;
         	     	      for(i = 0; i < ncolb; i++) 
	    			 *ctemp += (*atemp++) * alpha * (*btemp1++);
	 	     	      ctemp++;
      		  	   }
			   a++;
   	    		}
			break;

	 /* form C := alpha * A' * B' + beta * C */
	 case TRANSP:   for(i = 0; i < ncola; i++) {
			   btemp1 = b;
			   for (l = 0; l < nrowb; l++) {
      	       		      temp = ZERO;
	 		      for(j = 0; j < nrowa; j++) 
			         temp += a[j][i] * (*btemp1++);
	    		      *ctemp++ += alpha * temp;
			   }
   	    		}
			break;
      }
   }
}



/***********************************************************************
 *                                                                     *
 *                        enorm()                                      *
 *  a C translation of the Fortran-77 version by Burton, Garbow,       *
 *  Hillstrom and More of Argonne National Laboratory.                 *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   given an n-vector x, this function calculates the Euclidean norm of x.
   The Euclidean norm is computed by accumulating the sum of squares in 
   three different sums.  The sums of squares for the small and large 
   components are scaled so that no overflows occur.  Non-destructive 
   underflows are permitted.  Underflows and overflows do not occur in the 
   computation of the unscaled sum of squares for the intermediate components. 
   The definitions of small, intermediate and large components depend on two
   constants, rdwarf and rgiant.  The restrictions on these constants are 
   that rdwarf**2 not underflow and rgiant**2 not overflow.  The constants
   given here are suitable for every known computer.
   The function returns the Euclidean norm of vector x in double precision.


   Parameters
   ----------

   n         number of elements in vector x
   x         linear array of vector x whose Euclidean norm is to be 
		calculated

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

double svdpack::enorm(long n, double *x)

{
   double norm2, agiant, floatn, s1, s2, s3, xabs, x1max, x3max;
   long i;

   s1     = ZERO;
   s2     = ZERO;
   s3     = ZERO;
   x1max  = ZERO;
   x3max  = ZERO;
   floatn = (double)n;
   agiant = RGIANT / floatn;

   for (i = 0; i < n; i++) {
      xabs = fabs(x[i]);
      /* summing components of vector that need no scaling */
      if (xabs > RDWARF && xabs < agiant)
	 s2 += xabs * xabs;
      else {
	 /* underflow... */
	 if (xabs <= RDWARF) {
            if (xabs > x3max) {
	       s3 = ONE + s3 * (x3max/xabs) * (x3max/xabs);
	       x3max = xabs;
	    }
	    else if (xabs != 0)
	       s3 += (xabs/x3max) * (xabs/x3max);
	 }
	 /* overflow... */
	 else {
	    /* summing large components of vector */
	    if (xabs <= x1max)
	       s1 += (xabs/x1max) * (xabs/x1max);
            else {
	       s1 = ONE + s1 * (x1max/xabs) * (x1max/xabs);
	       x1max = xabs;
	    }
	 }
      }
   }
   if (s1 != ZERO)
      norm2 = x1max * sqrt(s1 + (s2/x1max) / x1max);
   else if (s2 != ZERO) {
      if (s2 >= x3max)
	 norm2 = sqrt(s2 * (ONE + (x3max/s2) * (x3max*s3)));
      else 
	 norm2 = sqrt(x3max * ((s2/x3max) + (x3max*s3)));
   }
   else
      norm2 = x3max * sqrt(s3);
   return(norm2);
}


/***********************************************************************
 *                                                                     *
 *                          dtbmv()                                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   The function performs one of the matrix-vector operations

	x := A * x,  or  x := A' * x,

   where A is an upper-triangular matrix.

   Parameters
   ----------

   trans     if trans = TRANSP, A' is to be used in the multiplication
             if trans = NTRANSP, A is to be used in the multiplication

   n         number of rows of matrix A; n must be at least 0.  Unchanged
	     upon exit.

   k         number of super-diagonals of matrix A

   a         2-dimensional array whose leading n by (k + 1) part must 
	     contain the upper triangular band part of the matrix of 
	     coefficients, supplied row by row, with the leading diagonal
	     of the matrix in column (k + 1) of the array, the first super-
	     diagonal starting at position 2 in column k, and so on.
	     The top left k by k triangle of the array A is not referenced.

   x         linear array of dimension of at least n.  Before entry,
	     x must contain the n elements of vector x.  On exit, x is
	     overwritten with the transformed vector x.


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

   MISC      imax  

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

void svdpack::dtbmv(long trans, long n, long k, double **a, double *x)

{
   long info, j, i, l, end;
   double temp;

   info = 0;
   if      ( trans != TRANSP && trans != NTRANSP )   info = 1;
   else if ( n < 0 )                                 info = 2;
   else if ( k < 0 )                                 info = 3;

   if (info) {
      fprintf(stderr, "%s %1ld %s\n",
      "*** ON ENTRY TO DTBMV, PARAMETER NUMBER",info,"HAD AN ILLEGAL VALUE");
      return; //exit(info);
   }

   switch(trans) {
      case NTRANSP:  for (j = 0; j < n; j++) {
                        temp = x[j];
                        l = k - j;
                        for (i = imax(0, j - k); i < j; i++) 
                           x[i] += temp * a[j][l+i];
                        x[j] *= a[j][k];
                     }
		     break;

      case TRANSP:   for (j = n - 1; j >= 0; j--) {
			temp = x[j] * a[j][k];
			l = k - j;
			end = imax(0, j - k);
			for (i = j - 1; i >= end; i--)
			   temp += x[i] * a[j][l+i];
                        x[j] = temp;
		     }
		     break;
   }
}


/************************************************************** 
 * Function interchanges two vectors		     	      *
 * Based on Fortran-77 routine from Linpack by J. Dongarra    *
 **************************************************************/ 

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

{
   long i;
   double dtemp;

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

/* am keeping these temporarily here will do away with later on */

#define               MAXIT     30
#define               CASE1     1
#define               CASE2     2
#define               CASE3     3
#define               CONVERGE  4

/***********************************************************************
 *                                                                     *
 *                        qriter2()                                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   This function reduces an upper bidiagonal matrix B to diagonal form.
   It is a C translation of a portion of DSVDC from Linpack.  In this