s_las2.cc

矩阵奇异分解(svd)最新c++版本

	    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;
		  r = (d[i] - g) * s + 2.0 * c * b;
		  p = s * r;
		  d[i+1] = g + p;
		  g = c * r - b;
		  f = bnd[i+1];
		  bnd[i+1] = s * bnd[i] + c * f;
		  bnd[i] = c * bnd[i] - s * f;
		  i--;
	       }
	    }       /* end while (underflow != FALSE && i >= l) */
	    /*........ recover from underflow .........*/
	    if (underflow) {
	       d[i+1] -= p;
	       e[m] = 0.0;
	    }
	    else {
	       d[l] -= p;
	       e[l] = g;
	       e[m] = 0.0;
	    }
	 } 		       		   /* end if (m != l) */
	 else {

            /* order the eigenvalues */
	    exchange = TRUE;
	    if (l != 0) {
	       i = l;
	       while (i >= 1 && exchange == TRUE) {
	          if (p < d[i-1]) {
		     d[i] = d[i-1];
		     bnd[i] = bnd[i-1];
	             i--;
	          }
	          else exchange = FALSE;
	       }
	    }
	    if (exchange) i = 0;
	    d[i] = p;
	    bnd[i] = f; 
	    iteration = 31;
	 }
      }			       /* end while (iteration <= 30) */
   }				   /* end for (l=0; l<n; l++) */
   return;
}


/***********************************************************************
 *                                                                     *
 *			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 s_las2::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);
      mopb(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;
}



/***********************************************************************
 *                                                                     *
 *                          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 s_las2::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;
}


/***********************************************************************
 *                                                                     *
 *				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 s_las2::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;
}


/***************************************************************** 
 * Function finds the index of element having max. absolute value*
 * based on FORTRAN 77 routine from Linpack by J. Dongarra       *
 *****************************************************************/ 

long s_las2::idamax(long n,double *dx,long incx)

{
   long ix,i,imax;
   double dtemp, dmax;

   if (n < 1) return(-1);
   if (n == 1) return(0);
   if (incx == 0) return(-1);

   if (incx < 0) ix = (-n+1) * incx;
   else ix = 0;
   imax = ix;
   dx += ix;
   dmax = fabs(*dx);
   for (i=1; i < n; i++) {
      ix += incx;
      dx += incx;
      dtemp = fabs(*dx);
      if (dtemp > dmax) {
	 dmax = dtemp;
	 imax = ix;
      }
   }
   return(imax);
}


/***********************************************************************
 *                                                                     *
 *			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 s_las2::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]))