tms2.c

求矩阵奇异分解svd算法

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

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

long tql2(long offset, 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 = offset; 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 = offset; 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+offset; 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 = offset; j < n; j ++) {
	     p = z[i][j];
	     z[i][j] = z[k][j];
	     z[k][j] = p;
	  }
      }   
   }
   return(0);
}		/*...... end main ............................*/
#include <math.h>
double dmax(double, double);
double dmin(double, double);

/************************************************************** 
 *							      *
 * Function finds sqrt(a^2 + b^2) without overflow or         *
 * destructive underflow.				      *
 *							      *
 **************************************************************/ 
/************************************************************** 

   Funtions used
   -------------

   UTILITY	dmax, dmin

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

double 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);
}
#include <math.h>
/***********************************************************************
 *                                                                     *
 *				random()                               *
 *                        (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 random(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);
}
#include <stdio.h>
#include "tmsc.h"

extern long ncol,nrow;
extern long *pointr,*rowind;
extern double *value; 

void opat(double *x, double *y)
{
/*  multiplication of transpose (nrow by ncol sparse matrix A)
    by the vector x, store in y  */

  long i,j,last;

  for(i=0; i<ncol; i++)
     y[i]=ZERO;

  for(i=0; i<ncol; i++) {
     last = pointr[i+1];
     for(j=pointr[i]; j<last; j++)
        y[i] += value[j]*x[rowind[j]];
  }

return;
}

#include <stdio.h>
#include "tmsc.h"


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

/***********************************************************************
 *                                                                     *
 *                         dgemm3()                                    *
 *                                                                     *
 * 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
   -----------

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

   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.
   ira,ica  row and column indices of matrix a, where mxk part starts.

   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.
   irb,icb  row and column indices of matrix b, where kxn part starts.

   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).
   irc,icc  row and column indices of matrix c, where mxn part starts.

***********************************************************************/
void dgemm3(long transa, long transb, long m, long n, long k, 
            double alpha, double **a, long ira, long ica, double **b,
            long irb, long icb, double beta, double **c, long irc,
            long icc)
{
   long info;
   long i, j, l, nrowa, ncola, nrowb, ncolb;
   double temp;

   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 DGEMM3, PARAMETER NUMBER",info,"HAD AN ILLEGAL VALUE");
      exit(info);
   }

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

   if (alpha == ZERO) {
      if (beta == ZERO) 
         for (j = 0; j < n; j++)
            for (i = 0; i < m; i++) c[icc+j][irc+i] = ZERO;
      else 
         for (j = 0; j < n; j++)
             for (i = 0; i < m; i++) c[icc+j][irc+i] *= beta;
      return;
   }
   if (!transb) { 

      switch(transa) {

	 /* form C := alpha * A * B + beta * C */
	 case NTRANSP: for(j = 0; j < n; j++) {
                          for(i=0; i<m; i++) c[icc+j][irc+i]=0.0;
                          for(l=0; l<k; l++) 
                             if(b[icb+j][irb+l]!=0.0) {
                               temp = alpha*b[icb+j][irb+l];
                               for(i=0; i<m; i++) 
                                c[icc+j][irc+i] += temp*a[ica+l][ira+i];
      		  	      }
                        }
			break;

	 /* form C := alpha * A' * B + beta * C */
	 case TRANSP:   for(j = 0; j < n; j++) {
                           for(i=0; i<m; i++) {
                              temp = 0.0;
      		  	      for(l = 0; l< k;l++) 
                                 temp += a[ica+i][ira+l]*b[icb+j][irb+l];
                              if(beta==0.0) c[icc+j][irc+i]=alpha*temp; 
                              else 
                       c[icc+j][irc+i] = alpha*temp + beta*c[icc+j][irc+i];
      		  	   }
   	    		}
			break;
      }
   }
   else { 
      switch(transa) {

	 /* form C := alpha * A * B' + beta * C */
	 case NTRANSP: for(j=0; j<n; j++) {
      		  	   for(i=0; i<m; i++) c[icc+j][irc+i]=ZERO;
      		  	   for(l=0; l<k; l++) {
	 	     	      temp = alpha*b[l+icb][j+irb];
         	     	      for(i=0; i<m; i++) 
	    			 c[j+icc][i+irc] += temp*a[l+ica][i+ira];
      		  	   }
   	    		}
			break;

	 /* form C := alpha * A' * B' + beta * C */
	 case TRANSP:   for(j=0; j<n; j++) {
			   for (i=0; i<m; i++) {
      	       		      temp = ZERO;
	 		      for(l=0; l<k; l++) 
				 temp += a[i+ica][l+ira]*b[l+icb][j+irb];
                              if(!beta) c[j+icc][i+irc] += alpha * temp;
                              else 
                                 c[j+icc][i+irc]= alpha*temp+
                                                  beta*c[j+icc][i+irc];
			   }
   	    		}
			break;
      }
   }