tms1.c

求矩阵奇异分解svd算法

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

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

   UTILITY	dmax, dmin

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

float pythag(float a, float b)

{
   float 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()                               *
 *                        (float 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 float precision random number between (0,1)

 ***********************************************************************/
float random(short *iy)

{
   static short m2 = 0;
   static short ia, ic, mic;
   static float halfm, s;

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

      /* compute multiplier and increment for linear congruential 
       * method */
      ia = 8 * (int)(halfm * atan(1.0) / 8.0) + 5;
      ic = 2 * (int)(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((float)(*iy) * s);
}
#include <stdio.h>
#include "tmsc.h"

extern int ncol,nrow;
extern int *pointr,*rowind;
extern float *value; 

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

  int i,j,upper;

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

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

  return;
}

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


/***********************************************************************
 *                                                                     *
 *                         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 mxn 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 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 the mxn part is stored.

***********************************************************************/
void dgemm3(int transa, int transb, int m, int n, int k, 
            float alpha, float **a, int ira, int ica, float **b,
            int irb, int icb, float beta, float **c, int irc,
            int icc)
{
   int info;
   int i, j, l, nrowa, ncola, nrowb, ncolb;
   float 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 %1d %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;
      }
   }
}
#include <stdio.h>
#include <math.h>
#include "tmsc.h"

float dasum(int n, float *dx, int incx)
{
  /**************************************************************
   *  Function forms the sum of the absolute values.            *
   *  Uses unrolled loops for increment equal to one.           *
   *  Based on Fortran-77 routine from Linpack by J.Dongarra.
   **************************************************************/

  float dtemp,dsum;
  int   i,ix,m,mp1;

  dsum = ZERO;
  dtemp = ZERO;
  if(n <= 0) return;
  if(incx != 1) {

  /* code for increment not equal to 1 */

    ix = 0;
    if(incx < 0) ix = (-n+1)*incx + 1;
    for(i=0; i<n; i++) {
       dtemp += fabs(dx[ix]);
       ix    += incx;
    }
    dsum = dtemp;
    return(dsum);
  }  

  /* code for increment equal to 1 */
  
  /* clean-up loop */

  m = n % 6;
  if(m) {
    for(i=0; i<m; i++)
       dtemp += fabs(dx[i]);
  }
  if(n>=6) {
    for(i=m; i<n; i+=6)
       dtemp += fabs(dx[i]) + fabs(dx[i+1]) + fabs(dx[i+2]) +
                fabs(dx[i+3]) + fabs(dx[i+4]) + fabs(dx[i+5]);
  }
  dsum = dtemp;
  return(dsum);
}
#include <math.h>

#define ZERO 0.0e+0
#define ONE 1.0e+0
float fsign(float, float);


/***********************************************************************
 *                                                                     *
 *                              tred2()                                *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

  Description
  -----------
  
  tred2() is a translation of the algol procedure TRED2, Num. Math. 11, 
  181-195 (1968) by Martin, Reinsch, and Wikinson.  Handbook for Auto.
  Comp., Vol. II- Linear Algebra, 212-226 (1971)

  This subroutine reduces a real symmetric matrix to a symmetric
  tridiagonal matrix using and accumulating orthogonal similarity
  transformations.

  Arguments
  ---------

  (input)
  offset index of the leading element of the matrix to be
         tridiagonalized. The matrix tridiagonalized should be 
         stored in a[offset:n-1, offset:n-1]

  n	 order of the matrix

  a	 contains the real symmetric input matrix. Only the upper
	 triangle of the matrix need be supplied

  (output)
  d	 contains the diagonal elements of the tridiagonal matrix.
  
  e	 contains the subdiagonal elements of the tridiagonal matrix
	 in its first n-1 positions.

  z	 contains the orthogonal transformation matrix produced in the
	 reduction.

  a and z may coincide. If distinct, a is unaltered.

  Functions used:
  UTILITY: fsign

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

void tred2(int offset, int n, float **a, float *d, float *e, float **z)
{
 int jj,ii,i,j,k,l, jp1;
 float *zptr,*aptr,h, scale, f, g,  hh, tmp;

 int i1;

 for (i=offset;i<n;i++) 
  { 
   for (j=i;j<n;j++)
     {
      z[j][i]=a[i][j];   /*fix this later?.. the rest of the routine 
                           assumes that z has the lower triangular part
                           of the symmetric matrix */
     }
   d[i]=a[i][n-1];
  }


  if (n==1) 
   {
    for (i=offset;i<n;i++)
     {
       d[i]=z[n-1][i];
       z[n-1][i]=ZERO;
     }
    z[n-1][n-1]=ONE;
    e[1]=ZERO;
    return;
   }

  /*for i = n step -1 until 2 do*/

  for (ii=3;ii<n+2-offset;ii++)
   {
     i=n+2-ii;
     l=i-1;
     h=ZERO; 
     scale=ZERO;

    /*scale row (algol tol then not needed)*/
     if (l>=1)