las1.c

求矩阵奇异分解svd算法

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

   BLAS		daxpy, datx, dcopy, ddot, dscal
   USER		store, opb
   LAS		startv

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

void stpone(long n, double *wrkptr[])

{
   double 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 long j, ierr;
extern double eps;
double ddot(long, double *,long, double *, long);
void   daxpy(long, double, double *,long, double *, long);
void   dcopy(long, double *, long, double *, long);
static double random(long *);
void   store(long, long, long, double *);
void   opb(long, double *, double *);

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

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

double startv(long n, double *wptr[])

{
   double rnm2, *r, t;
   long irand;
   long id, i;

   /* get initial vector; default is random */
   rnm2 = ddot(n, wptr[0], 1, wptr[0], 1);
   irand = 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()                               *
 *                        (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)

 ***********************************************************************/
static 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 <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>
extern double tol, eps34, eps;
extern long j, neig;
long   idamax(long, double *, long);
double dmin(double, double);

/***********************************************************************
 *                                                                     *
 *			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(long *enough, double endl, double endr, 
		 double *ritz, double *bnd)
{
   long mid, i;
   double 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 long ierr;
double pythag(double, double);
double fsign(double, double);

/***********************************************************************
 *                                                                     *
 *				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(long n, double d[], double e[], double bnd[])

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

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

   double 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;