svdpack.cc

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

   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
            leading m by k part of a must contain the matrix A. Otherwise,
	    the leading k by m part of a must contain  the matrix A.

   b        matrix B as a 2-dimensional array.  When transb = NTRANSP, the
            leading k by n part of a must contain the matrix B. Otherwise,
	    the leading n by k part of a 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 2-dimensional array.  On entry, the leading
	    m by n part of c must contain the matrix C, except when
	    beta = 0.  In that case, c need not be set on entry. 
	    On exit, c is overwritten by the m by n matrix
	    (alpha * op(A) * op(B) + beta * C).

 ***********************************************************************/
void svdpack::dgemm2(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;
   double temp, *atemp;

   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 %1d %s\n",
      "*** ON ENTRY TO DGEMM2, PARAMETER NUMBER",info,"HAD AN ILLEGAL VALUE");
      exit(info);
   }

   if (transa) {
      nrowa = k;
      ncola = m;
   }
   else { 
      nrowa = m;
      ncola = k;
   }
   if (transb) {
      nrowb = n;
      ncolb = k;
   }
   else {
      nrowb = k;
      ncolb = n;
   }
   if (!m || !n || ((alpha == ZERO || !k) && beta == ONE))
      return;

   if (alpha == ZERO) {
      if (beta == ZERO) 
         for (i = 0; i < m; i++)
            for (j = 0; j < n; j++) c[i][j] = ZERO;

      else if (beta != ONE)
         for (i = 0; i < m; i++)
            for (j = 0; j < n; j++) c[i][j] *= beta;
      return;
   }

   if (beta == ZERO)
      for (i = 0; i < m; i++)
         for (j = 0; j < n; j++) c[i][j] = ZERO;

   else if (beta != ONE)
      for (i = 0; i < m; i++)
         for (j = 0; j < n; j++) c[i][j] *= beta;

   if (!transb) { 

      switch(transa) {

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

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

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

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



/***********************************************************************
 *                                                                     *
 *                          dsbmv()                                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   The function performs the matrix-vector operation

	  y := alpha * A * y + beta * y,

   where alpha and beta are scalars, x and y are n-element vectors and
   A is an n by n symmetric band matrix, with k super-diagonals.

   Parameters
   ----------

   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 symmetric matrix,
	     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.  Unchanged on exit.

   y         linear array of dimension of at least n.  Before entry,
	     y must contain the n elements of vector y.  On exit, y is
	     overwritten by the updated vector y.


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

   MISC      imax  

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

void svdpack::dsbmv(long n, long k, double alpha, double **a, 
           double *x, double beta, double *y)

{
   long info, j, i, l;
   double *ptrtemp, temp1, temp2;

   info = 0;
   if ( n < 0 )                                      info = 1;
   else if ( k < 0 )                                 info = 2;

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

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

   ptrtemp = y; 

   /* form y := beta * y */
   if (beta == ZERO) 
      for (i = 0; i < n; i++) *ptrtemp++ = ZERO;
   else if (beta != ONE) 
      for (i = 0; i < n; i++) *ptrtemp++ *= beta;

   if (alpha == ZERO) return;

   for (j = 0; j < n; j++) {
      temp1 = alpha * x[j];
      temp2 = ZERO;
      l = k - j;
      for (i = imax(0, j - k); i < j; i++) {
         y[i] = y[i] + temp1 * a[j][l+i];
         temp2 = temp2 + a[j][l+i] * x[i];
      }
      y[j] = y[j] + temp1 * a[j][k] + alpha * temp2;
   }
}



/***********************************************************************
 *                                                                     *
 *				tql2()   			       *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   tql2() is a translation of a Fortran version of the Algol
   procedure TQL2, Num. Math. 11, 293-306(1968) by Dowdler, Martin, 
   Reinsch and Wilkinson.
   Handbook for Auto. Comp., vol.II-Linear Algebra, 227-240(1971).  

   This function finds the eigenvalues and eigenvectors of a symmetric
   tridiagonal matrix by the QL method.


   Arguments
   ---------

   (input)                                                             
   n      order of the symmetric tridiagonal matrix           
   d      contains the diagonal elements of the input matrix        
   e      contains the subdiagonal elements of the input matrix in its
            first n-1 positions.
   z      contains the identity matrix				    
                                                                   
   (output)                                                       
   d      contains the eigenvalues in ascending order.  if an error
            exit is made, the eigenvalues are correct but unordered for
            for indices 0,1,...,ierr.				   
   e      has been destroyed.					  
   z      contains orthonormal eigenvectors of the symmetric   
            tridiagonal (or full) matrix.  if an error exit is made,
            z contains the eigenvectors associated with the stored 
          eigenvalues.					
   ierr   set to zero for normal return, j if the j-th eigenvalue has
            not been determined after 30 iterations.		    


   Functions used
   --------------
   UTILITY	fsign
   MISC		pythag

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

long svdpack::tql2(long n, double *d, double *e, double **z)

{
   long j, last, l, l1, l2, m, i, k, iteration;
   double tst1, tst2, g, r, s, s2, c, c2, c3, p, f, h, el1, dl1;
   if (n == 1) return(0);
   f = ZERO;
   last = n - 1;
   tst1 = ZERO;
   e[last] = ZERO;

   for (l = 0; l < n; l++) {
      iteration = 0;
      h = fabs(d[l]) + fabs(e[l]);
      if (tst1 < h) tst1 = h;

      /* look for small sub-diagonal element */
      for (m = l; m < n; m++) {
	 tst2 = tst1 + fabs(e[m]);
	 if (tst2 == tst1) break;
      }
      if (m != l) {
	 while (iteration < 30) {
	    iteration += 1;

            /*  form shift */
	    l1 = l + 1;
	    l2 = l1 + 1;
	    g = d[l];
            p = (d[l1] - g) / (2.0 * e[l]);
	    r = pythag(p, ONE);
	    d[l] = e[l] / (p + fsign(r, p));
	    d[l1] = e[l] * (p + fsign(r, p));
	    dl1 = d[l1];
	    h = g - d[l];
	    if (l2 < n) 
	       for (i = l2; i < n; i++) d[i] -= h;
            f += h;

	    /* QL transformation */
	    p = d[m];
	    c = ONE;
	    c2 = c;
	    el1 = e[l1];
	    s = ZERO;
	    i = m - 1;
	    while (i >= l) {
	       c3 = c2;
	       c2 = c;
	       s2 = s;
	       g = c * e[i];
	       h = c * p;
	       r = pythag(p, e[i]);
	       e[i + 1] = s * r;
	       s = e[i] / r;
	       c = p / r;
	       p = c * d[i] - s * g;
	       d[i + 1]= h + s * (c * g + s * d[i]);

	       /*  form vector */
	       for (k = 0; k < n; k ++) {
	          h = z[i + 1][k];
	          z[i + 1][k] = s * z[i][k] + c * h;
	          z[i][k] = c * z[i][k] - s * h;
	       }
	       i--;
	    }
	    p = -s * s2 * c3 *el1 * e[l] / dl1;
	    e[l] = s * p;
	    d[l] = c * p;
	    tst2 = tst1 + fabs(e[l]);
	    if (tst2 <= tst1) break;
	    if (iteration == 30) 
	       return(l);
         }
      }
      d[l] += f;
   }

   /* order the eigenvalues */
   for (l = 1; l < n; l++) {
      i = l - 1;
      k = i;
      p = d[i];
      for (j = l; j < n; j++) {
	 if (d[j] < p) {
	    k = j;
	    p = d[j];
	 }
      }
      /* ...and corresponding eigenvectors */
      if (k != i) {
	 d[k] = d[i];
	 d[i] = p;
	  for (j = 0; j < n; j ++) {
	     p = z[i][j];
	     z[i][j] = z[k][j];
	     z[k][j] = p;
	  }
      }   
   }
   return(0);
}


double svdpack::pythag(double a, double b)

{
   double p, r, s, t, u, temp;

   p = dmax(fabs(a), fabs(b));
   if (p != 0.0) {
      temp = dmin(fabs(a), fabs(b)) / p;
      r = temp * temp; 
      t = 4.0 + r;
      while (t != 4.0) {
	 s = r / t;
	 u = 1.0 + 2.0 * s;
	 p *= u;
	 temp = s / u;
	 r *= temp * temp;
	 t = 4.0 + r;
      }
   }
   return(p);
}


/**************************************************************
 *                                                            *
 * multiplication of A'A by the transpose of an nc by n       *
 * matrix X.  A is nrow by ncol and is stored using the       *
 * Harwell-Boeing compressed column sparse matrix format.     *
 * The transpose of the n by nc product matrix Y is returned  *
 * in the two-dimensional array y.                            *
 *                                                            *
 *              Y = A'A * X'  and  y := Y'                    *
 *                                                            *
 **************************************************************/

void svdpack::opm(long nrow, long ncol, long nc, double **x, double **y)
{
   long i, j, k, end;
   
   mxvcount  += nc;
   mtxvcount += nc;

   for (i = 0; i < nc; i++) 
      for (j = 0; j < nrow; j++) ztempp[i][j] = ZERO;
   for (i = 0; i < nc; i++) 
      for (j = 0; j < ncol; j++) y[i][j] = ZERO;

   /* multiply by sparse matrix A */
   for (i = 0; i < ncol; i++) {
      end = pointr[i+1];
      for (j = pointr[i]; j < end; j++)
         for (k = 0; k < nc; k++)
	    ztempp[k][rowind[j]] += value[j] * x[k][i]; 
   }

   /* multiply by transpose of A (A') */
   for (k = 0; k < nc; k++)
      for (i = 0; i < ncol; i++) {
         end = pointr[i+1];
         for (j = pointr[i]; j < end; j++)