las1.c

求矩阵奇异分解svd算法

void lanso(long n, long lanmax, long maxprs, double endl,
	   double endr, double *ritz, double *bnd, double *wptr[])

{
   double *r, *alf, *eta, *oldeta, *bet, *wrk;
   long ll, first, last, ENOUGH, id1, id2, id3, i, l;

   r = wptr[0];
   alf = wptr[6];
   eta = wptr[7];
   oldeta = wptr[8];
   bet = wptr[9];
   wrk = wptr[5];
   j = 0;

   /* take the first step */
   stpone(n, wptr);
   if (!rnm || ierr) return;
   eta[0] = eps1;
   oldeta[0] = eps1;
   ll = 0;
   first = 1;
   last = imin(maxprs + imax(8,maxprs), lanmax);
   ENOUGH = FALSE;
   id1 = 0;
   while (id1 < maxprs && !ENOUGH) {
      if (rnm <= tol) rnm = 0.0;

      /* the actual lanczos loop */
      lanczos_step(n, first, last, wptr, alf, eta, oldeta, bet, &ll,
		   &ENOUGH);
      if (ENOUGH) j = j - 1;
      else j = last - 1;
      first = j + 1;
      bet[j+1] = rnm;

      /* analyze T */
      l = 0;
      for (id2 = 0; id2 < j; id2++) {
	 if (l > j) break;
         for (i = l; i <= j; i++) if (!bet[i+1]) break;
	 if (i > j) i = j;

	 /* now i is at the end of an unreduced submatrix */
	 dcopy(i-l+1, &alf[l],   1, &ritz[l],  -1);
	 dcopy(i-l,   &bet[l+1], 1, &wrk[l+1], -1);

	 imtqlb(i-l+1, &ritz[l], &wrk[l], &bnd[l]);

	 if (ierr) {
	    printf("IMTQLB FAILED TO CONVERGE (IERR = %d)\n",
		    ierr);
	    printf("L = %d    I = %d\n", l, i);
            for (id3 = l; id3 <= i; id3++) 
	       printf("%d  %lg  %lg  %lg\n",
	       id3, ritz[id3], wrk[id3], bnd[id3]);
	 }
         for (id3 = l; id3 <= i; id3++) 
	    bnd[id3] = rnm * fabs(bnd[id3]);
	 l = i + 1;
      }

      /* sort eigenvalues into increasing order */
      dsort2((j+1) / 2, j + 1, ritz, bnd);

      /* massage error bounds for very close ritz values */
      error_bound(&ENOUGH, endl, endr, ritz, bnd);

      /* should we stop? */
      if (neig < maxprs) {
	 if (!neig) last = first + 9;
	 else last = first + imax(3, 1 + ((j-5) * (maxprs-neig)) / neig);
	 last = imin(last, lanmax);
      }
      else ENOUGH = TRUE;
      ENOUGH = ENOUGH || first >= lanmax;
      id1++;
   }
   store(n, STORQ, j, wptr[1]);
   return;
}
#include <math.h>
#define STORQ 1
#define RETRQ 2
#define STORP 3
#define RETRP 4
#define TRUE  1
#define FALSE 0
#define MAXLL 2

extern double rnm, anorm, tol, eps, eps1, reps, eps34;
extern long ierr, j;
double ddot(long, double *,long, double *, long);
void   dscal(long, double, double *,long);
void   daxpy(long, double, double *,long, double *, long);
void   datx(long, double, double *,long, double *, long);
void   dcopy(long, double *, long, double *, long);
void   purge(long, long, double *, double *, double *, double *,
	     double *, double *, double *);
void   ortbnd(double *, double *, double *, double *);
double startv(long, double *[]);
void   store(long, long, long, double *);
long   imin(long, long);
long   imax(long, long);

/***********************************************************************
 *                                                                     *
 *			lanczos_step()                                 *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   Function embodies a single Lanczos step

   Arguments 
   ---------

   (input)
   n        dimension of the eigenproblem for matrix B
   first    start of index through loop				      
   last     end of index through loop				     
   wptr	    array of pointers pointing to work space		    
   alf	    array to hold diagonal of the tridiagonal matrix T
   eta      orthogonality estimate of Lanczos vectors at step j   
   oldeta   orthogonality estimate of Lanczos vectors at step j-1
   bet      array to hold off-diagonal of T                     
   ll       number of intitial Lanczos vectors in local orthog. 
              (has value of 0, 1 or 2)			
   enough   stop flag			

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

   BLAS		ddot, dscal, daxpy, datx, dcopy
   USER		store
   LAS		purge, ortbnd, startv
   UTILITY	imin, imax

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

void lanczos_step(long n, long first, long last, double *wptr[],
		  double *alf, double *eta, double *oldeta,
		  double *bet, long *ll, long *enough)

{
   double t, *mid;
   long i;

   for (j=first; j<last; j++) {
      mid     = wptr[2];
      wptr[2] = wptr[1];
      wptr[1] = mid;
      mid     = wptr[3];
      wptr[3] = wptr[4];
      wptr[4] = mid;

      store(n, STORQ, j-1, wptr[2]);
      if (j-1 < MAXLL) store(n, STORP, j-1, wptr[4]);
      bet[j] = rnm;

      /* restart if invariant subspace is found */
      if (!bet[j]) {
	 rnm = startv(n, wptr);
	 if (ierr) return;
	 if (!rnm) *enough = TRUE;
      }
      if (*enough) break;

      /* take a lanczos step */
      t = 1.0 / rnm;
      datx(n, t, wptr[0], 1, wptr[1], 1);
      dscal(n, t, wptr[3], 1);
      opb(n, wptr[3], wptr[0]);
      daxpy(n, -rnm, wptr[2], 1, wptr[0], 1);
      alf[j] = ddot(n, wptr[0], 1, wptr[3], 1);
      daxpy(n, -alf[j], wptr[1], 1, wptr[0], 1);

      /* orthogonalize against initial lanczos vectors */
      if (j <= MAXLL && (fabs(alf[j-1]) > 4.0 * fabs(alf[j])))
	 *ll = j;  
      for (i=0; i < imin(*ll, j-1); i++) {
	 store(n, RETRP, i, wptr[5]);
	 t = ddot(n, wptr[5], 1, wptr[0], 1);
	 store(n, RETRQ, i, wptr[5]);
         daxpy(n, -t, wptr[5], 1, wptr[0], 1);
	 eta[i] = eps1;
	 oldeta[i] = eps1;
      }

      /* extended local reorthogonalization */
      t = ddot(n, wptr[0], 1, wptr[4], 1);
      daxpy(n, -t, wptr[2], 1, wptr[0], 1);
      if (bet[j] > 0.0) bet[j] = bet[j] + t;
      t = ddot(n, wptr[0], 1, wptr[3], 1);
      daxpy(n, -t, wptr[1], 1, wptr[0], 1);
      alf[j] = alf[j] + t;
      dcopy(n, wptr[0], 1, wptr[4], 1);
      rnm = sqrt(ddot(n, wptr[0], 1, wptr[4], 1));
      anorm = bet[j] + fabs(alf[j]) + rnm;
      tol = reps * anorm;

      /* update the orthogonality bounds */
      ortbnd(alf, eta, oldeta, bet);

      /* restore the orthogonality state when needed */
      purge(n,*ll,wptr[0],wptr[1],wptr[4],wptr[3],wptr[5],eta,oldeta);
      if (rnm <= tol) rnm = 0.0;
   }
   return;
}
extern double rnm, eps, eps1, reps, eps34;
extern long j;
void   dswap(long, double *, long, double *, long);

/***********************************************************************
 *                                                                     *
 *                          ortbnd()                                   *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   Funtion updates the eta recurrence

   Arguments 
   ---------

   (input)
   alf      array to hold diagonal of the tridiagonal matrix T         
   eta      orthogonality estimate of Lanczos vectors at step j        
   oldeta   orthogonality estimate of Lanczos vectors at step j-1     
   bet      array to hold off-diagonal of T                          
   n        dimension of the eigenproblem for matrix B		    
   j        dimension of T					  
   rnm	    norm of the next residual vector			 
   eps1	    roundoff estimate for dot product of two unit vectors

   (output)
   eta      orthogonality estimate of Lanczos vectors at step j+1     
   oldeta   orthogonality estimate of Lanczos vectors at step j        


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

   BLAS		dswap

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

void ortbnd(double *alf, double *eta, double *oldeta, double *bet)

{
   long i;
   if (j < 1) return;
   if (rnm) {
      if (j > 1) {
	 oldeta[0] = (bet[1] * eta[1] + (alf[0]-alf[j]) * eta[0] -
		      bet[j] * oldeta[0]) / rnm + eps1;
      }
      for (i=1; i<=j-2; i++) 
	 oldeta[i] = (bet[i+1] * eta[i+1] + (alf[i]-alf[j]) * eta[i] +
		      bet[i] * eta[i-1] - bet[j] * oldeta[i])/rnm + eps1;
   }
   oldeta[j-1] = eps1;
   dswap(j, oldeta, 1, eta, 1);  
   eta[j] = eps1;
   return;
}
#include <math.h>
#define STORQ 1
#define RETRQ 2
#define STORP 3
#define RETRP 4
#define TRUE  1
#define FALSE 0

extern double tol, rnm, eps, eps1, reps, eps34;
extern long j;
void   store(long, long, long, double *);
void   daxpy(long, double, double *, long, double *, long);
void   dcopy(long, double *, long, double *, long);
long   idamax(long, double *, long);
double ddot(long, double *, long, double *, long);

/***********************************************************************
 *                                                                     *
 *				purge()                                *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   Function examines the state of orthogonality between the new Lanczos
   vector and the previous ones to decide whether re-orthogonalization 
   should be performed


   Arguments 
   ---------

   (input)
   n        dimension of the eigenproblem for matrix B		       
   ll       number of intitial Lanczos vectors in local orthog.       
   r        residual vector to become next Lanczos vector            
   q        current Lanczos vector			           
   ra       previous Lanczos vector
   qa       previous Lanczos vector
   wrk      temporary vector to hold the previous Lanczos vector
   eta      state of orthogonality between r and prev. Lanczos vectors 
   oldeta   state of orthogonality between q and prev. Lanczos vectors
   j        current Lanczos step				     

   (output)
   r	    residual vector orthogonalized against previous Lanczos 
	      vectors
   q        current Lanczos vector orthogonalized against previous ones


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

   BLAS		daxpy,  dcopy,  idamax,  ddot
   USER		store

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

void purge(long n, long ll, double *r, double *q, double *ra,  
	   double *qa, double *wrk, double *eta, double *oldeta)

{
   double t, tq, tr, reps1;
   long k, iteration, flag, i;

   if (j < ll+2) return; 

   k = idamax(j - (ll+1), &eta[ll], 1) + ll;
   if (fabs(eta[k]) > reps) {
      reps1 = eps1 / reps;
      iteration = 0;
      flag = TRUE;
      while (iteration < 2 && flag) {
	 if (rnm > tol) {

	    /* bring in a lanczos vector t and orthogonalize both 
	     * r and q against it */
	    tq = 0.0;
	    tr = 0.0;
            for (i = ll; i < j; i++) {
	       store(n,  RETRQ,  i,  wrk);
	       t   = -ddot(n, qa, 1, wrk, 1);
	       tq += fabs(t);
	       daxpy(n,  t,  wrk,  1,  q,  1);
	       t   = -ddot(n, ra, 1, wrk, 1);
	       tr += fabs(t);
	       daxpy(n, t, wrk, 1, r, 1);
	    }
	    dcopy(n, q, 1, qa, 1);
	    t   = -ddot(n, r, 1, qa, 1);
	    tr += fabs(t);
	    daxpy(n, t, q, 1, r, 1);
	    dcopy(n, r, 1, ra, 1);
	    rnm = sqrt(ddot(n, ra, 1, r, 1));
	    if (tq <= reps1 && tr <= reps1 * rnm) flag = FALSE;
	 }
	 iteration++;
      }
      for (i = ll; i <= j; i++) { 
	 eta[i] = eps1;
	 oldeta[i] = eps1;
      }
   }
   return;
}
#include <math.h>

extern double rnm, anorm, tol, eps, reps;
extern long j, ierr;
void   daxpy(long, double, double *,long, double *, long);
void   datx(long, double, double *,long, double *, long);
void   dcopy(long, double *, long, double *, long);
double ddot(long, double *,long, double *, long);
void   dscal(long, double, double *,long);
double startv(long, double *[]);
void   opb(long, double *, double *);
void   store(long, long, long, double *);

/***********************************************************************
 *                                                                     *
 *                         stpone()                                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   Function performs the first step of the Lanczos algorithm.  It also
   does a step of extended local re-orthogonalization.

   Arguments 
   ---------

   (input)
   n      dimension of the eigenproblem for matrix B

   (output)
   ierr   error flag
   wptr   array of pointers that point to work space that contains
	    wptr[0]             r[j]
	    wptr[1]             q[j]
	    wptr[2]             q[j-1]
	    wptr[3]             p
	    wptr[4]             p[j-1]
	    wptr[6]             diagonal elements of matrix T