las1.c

求矩阵奇异分解svd算法

void stpone(int n, float *wrkptr[])

{
   float t, *alf;
   alf = wrkptr[6];

   /* get initial vector; default is random */
   rnm = startv(n, wrkptr);
   if (rnm == 0.0 || ierr != 0) return;

   /* normalize starting vector */
   t = 1.0 / rnm;
   datx(n, t, wrkptr[0], 1, wrkptr[1], 1);
   dscal(n, t, wrkptr[3], 1);

   /* take the first step */
   opb(n, wrkptr[3], wrkptr[0]);
   alf[0] = ddot(n, wrkptr[0], 1, wrkptr[3], 1);
   daxpy(n, -alf[0], wrkptr[1], 1, wrkptr[0], 1);
   t = ddot(n, wrkptr[0], 1, wrkptr[3], 1);
   daxpy(n, -t, wrkptr[1], 1, wrkptr[0], 1);
   alf[0] += t;
   dcopy(n, wrkptr[0], 1, wrkptr[4], 1);
   rnm = sqrt(ddot(n, wrkptr[0], 1, wrkptr[4], 1));
   anorm = rnm + fabs(alf[0]);
   tol = reps * anorm;
   return;
}
#include <stdio.h>
#include <math.h>
#define STORQ 1
#define RETRQ 2
#define STORP 3
#define RETRP 4
extern int j, ierr;
extern float eps;
float ddot(int, float *,int, float *, int);
void   daxpy(int, float, float *,int, float *, int);
void   dcopy(int, float *, int, float *, int);
static float random(short *);
void   store(int, int, int, float *);
void   opb(int, float *, float *);

/***********************************************************************
 *                                                                     *
 *                         startv()                                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   Function delivers a starting vector in r and returns |r|; it returns 
   zero if the range is spanned, and ierr is non-zero if no starting 
   vector within range of operator can be found.

   Parameters 
   ---------

   (input)
   n      dimension of the eigenproblem matrix B
   wptr   array of pointers that point to work space
   j      starting index for a Lanczos run
   eps    machine epsilon (relative precision)

   (output)
   wptr   array of pointers that point to work space that contains
	  r[j], q[j], q[j-1], p[j], p[j-1]
   ierr   error flag (nonzero if no starting vector can be found)

   Functions used
   --------------

   BLAS		ddot, dcopy, daxpy
   USER		opb, store
   MISC		random

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

float startv(int n, float *wptr[])

{
   float rnm2, *r, t;
   short irand;
   int id, i;

   /* get initial vector; default is random */
   rnm2 = ddot(n, wptr[0], 1, wptr[0], 1);
   irand = (short) 918273 + j;
   r = wptr[0];
   for (id = 0; id < 3; id++) {
      if (id > 0 || j > 0 || rnm2 == 0) 
	 for (i = 0; i < n; i++) r[i] = random(&irand);
      dcopy(n, wptr[0], 1, wptr[3], 1);

      /* apply operator to put r in range (essential if m singular) */
      opb(n, wptr[3], wptr[0]);
      dcopy(n, wptr[0], 1, wptr[3], 1);
      rnm2 = ddot(n, wptr[0], 1, wptr[3], 1);
      if (rnm2 > 0.0) break;
   }

   /* fatal error */
   if (rnm2 <= 0.0) {
      ierr = 8192;
      return(-1);
   }
   if (j > 0) {
      for (i = 0; i < j; i++) {
         store(n, RETRQ, i, wptr[5]);
	 t = -ddot(n, wptr[3], 1, wptr[5], 1);
	 daxpy(n, t, wptr[5], 1, wptr[0], 1);
      }

      /* make sure q[j] is orthogonal to q[j-1] */
      t = ddot(n, wptr[4], 1, wptr[0], 1);
      daxpy(n, -t, wptr[2], 1, wptr[0], 1);
      dcopy(n, wptr[0], 1, wptr[3], 1);
      t = ddot(n, wptr[3], 1, wptr[0], 1);
      if (t <= eps * rnm2) t = 0.0;
      rnm2 = t;
   }
   return(sqrt(rnm2));
}

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

 ***********************************************************************/
static 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 <math.h>
float dmax(float, float);
float dmin(float, float);

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

   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>
extern float tol, eps34, eps;
extern int j, neig;
int   idamax(int, float *, int);
float dmin(float, float);

/***********************************************************************
 *                                                                     *
 *			error_bound()                                  *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   Function massages error bounds for very close ritz values by placing 
   a gap between them.  The error bounds are then refined to reflect 
   this.


   Arguments 
   ---------

   (input)
   endl     left end of interval containing unwanted eigenvalues
   endr     right end of interval containing unwanted eigenvalues
   ritz     array to store the ritz values
   bnd      array to store the error bounds
   enough   stop flag


   Functions used
   --------------

   BLAS		idamax
   UTILITY	dmin

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

void error_bound(int *enough, float endl, float endr, 
		 float *ritz, float *bnd)
{
   int mid, i;
   float gapl, gap;

   /* massage error bounds for very close ritz values */
   mid = idamax(j + 1, bnd, 1);

   for (i=((j+1) + (j-1)) / 2; i >= mid + 1; i -= 1)
      if (fabs(ritz[i-1] - ritz[i]) < eps34 * fabs(ritz[i])) 
         if (bnd[i] > tol && bnd[i-1] > tol) {
	    bnd[i-1] = sqrt(bnd[i] * bnd[i] + bnd[i-1] * bnd[i-1]);
	    bnd[i] = 0.0;
         }
	 

   for (i=((j+1) - (j-1)) / 2; i <= mid - 1; i +=1 ) 
      if (fabs(ritz[i+1] - ritz[i]) < eps34 * fabs(ritz[i])) 
	 if (bnd[i] > tol && bnd[i+1] > tol) {
	    bnd[i+1] = sqrt(bnd[i] * bnd[i] + bnd[i+1] * bnd[i+1]);
	    bnd[i] = 0.0;
         }

   /* refine the error bounds */
   neig = 0;
   gapl = ritz[j] - ritz[0];
   for (i = 0; i <= j; i++) {
      gap = gapl;
      if (i < j) gapl = ritz[i+1] - ritz[i];
      gap = dmin(gap, gapl);
      if (gap > bnd[i]) bnd[i] = bnd[i] * (bnd[i] / gap);
      if (bnd[i] <= 16.0 * eps * fabs(ritz[i])) {
	 neig += 1;
	 if (!*enough) *enough = endl < ritz[i] && ritz[i] < endr;
      }
   }   
   return;
}
#include <math.h>
#define		TRUE	1
#define		FALSE	0
extern int ierr;
float pythag(float, float);
float fsign(float, float);

/***********************************************************************
 *                                                                     *
 *				imtqlb()			       *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   imtqlb() is a translation of a Fortran version of the Algol
   procedure IMTQL1, Num. Math. 12, 377-383(1968) by Martin and 
   Wilkinson, as modified in Num. Math. 15, 450(1970) by Dubrulle.  
   Handbook for Auto. Comp., vol.II-Linear Algebra, 241-248(1971).  
   See also B. T. Smith et al, Eispack Guide, Lecture Notes in 
   Computer Science, Springer-Verlag, (1976).

   The function finds the eigenvalues of a symmetric tridiagonal
   matrix by the implicit 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
          last n-1 positions.  e[0] is arbitrary	             

   (output)
   d      contains the eigenvalues in ascending order.  if an error
            exit is made, the eigenvalues are correct and ordered for
            indices 0,1,...ierr, but may not be the smallest eigenvalues.
   e      has been destroyed.					    
   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

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

void imtqlb(int n, float d[], float e[], float bnd[])

{
   int last, l, m, i, iteration;

   /* various flags */
   int exchange, convergence, underflow;	

   float b, test, g, r, s, c, p, f;

   if (n == 1) return;
   ierr = 0;
   bnd[0] = 1.0;
   last = n - 1;
   for (i = 1; i < n; i++) {
      bnd[i] = 0.0;
      e[i-1] = e[i];
   }
   e[last] = 0.0;
   for (l = 0; l < n; l++) {
      iteration = 0;
      while (iteration <= 30) {
	 for (m = l; m < n; m++) {
	    convergence = FALSE;
	    if (m == last) break;
	    else {
	       test = fabs(d[m]) + fabs(d[m+1]);
	       if (test + fabs(e[m]) == test) convergence = TRUE;
	    }
	    if (convergence) break;
	 }
	    p = d[l]; 
	    f = bnd[l]; 
	 if (m != l) {
	    if (iteration == 30) {
	       ierr = l;
	       return;
	    }
	    iteration += 1;
	    /*........ form shift ........*/
	    g = (d[l+1] - p) / (2.0 * e[l]);
	    r = pythag(g, 1.0);
	    g = d[m] - p + e[l] / (g + fsign(r, g));
	    s = 1.0;
	    c = 1.0;
	    p = 0.0;
	    underflow = FALSE;
	    i = m - 1;
	    while (underflow == FALSE && i >= l) {
	       f = s * e[i];
	       b = c * e[i];
	       r = pythag(f, g);
	       e[i+1] = r;
	       if (r == 0.0) underflow = TRUE;
	       else {
		  s = f / r;
		  c = g / r;
		  g = d[i+1] - p;