s_las1.cc

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

   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_las1::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;
}




/***********************************************************************
 *                                                                     *
 *                        ritvec()                                     *
 * 	    Function computes the singular vectors of matrix A	       *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   This function is invoked by landr() only if eigenvectors of the cyclic
   eigenproblem are desired.  When called, ritvec() computes the 
   singular vectors of A and writes the result to an unformatted file.


   Parameters
   ----------

   (input)
   n	      dimension of the eigenproblem for matrix B
   kappa      relative accuracy of ritz values acceptable as 
		eigenvalues of B (singular values of A)
   ritz       array of ritz values
   bnd        array of error bounds
   alf        array of diagonal elements of the tridiagonal matrix T
   bet        array of off-diagonal elements of T
   w1, w2     work space
   fp_out2    pointer to unformatted output file
   j	      number of Lanczos iterations performed

   (output)
   xv1        array of eigenvectors of B (singular vectors of A)
   ierr	      error code
	      0 for normal return from imtql2()
	      k if convergence did not occur for k-th eigenvalue in 
		imtql2()
   nsig       number of accepted ritz values based on kappa

   (local)
   s	      work array which is initialized to the identity matrix of
	      order (j + 1) upon calling imtql2().  After the call, 
	      s contains the orthonormal eigenvectors of the symmetric 
	      tridiagonal matrix T

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

   BLAS		dscal, dcopy, daxpy
   USER		store
   		imtql2

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

void s_las1::ritvec(long fpo2, FILE* fp, bool ascii, long n, 
		    double kappa, double *ritz, double *bnd,
		    double *alf, double *bet, double *w1, double *w2)

{
   long js, jsq, i, k, size, id, id2, tmp;
   double *s;

   js = j + 1;
   jsq = js * js;
   size = sizeof(double) * n;

   if(!(s = (double *) malloc (jsq * sizeof(double)))) {
      perror("MALLOC FAILED in RITVEC()");
      exit(errno);
   }

   /* initialize s as an identity matrix */
   for (i = 0; i < jsq; i++) s[i] = 0.0;
   for (i = 0; i < jsq; i+= (js + 1)) s[i] = 1.0;
   dcopy(js, alf, 1, w1, -1);
   dcopy(j, &bet[1], 1, &w2[1], -1);

   /* compute eigenvalues of T */
   imtql2(js, js, w1, w2, s);
   if (ierr) return;

   /* on return w1 contains eigenvalues in ascending order and s 
    * contains the corresponding eigenvectors */

   /* SUVRIT: Commented out the following header info and replaced by
      svdpack::write_header */
   //write(fpo2, (char *)&n, sizeof(n));
   //write(fpo2, (char *)&js, sizeof(js));
   //write(fpo2, (char *)&kappa, sizeof(kappa));
   nsig = 0;   
   for (k = 0; k < js; k++) {
      if (bnd[k] <=  kappa*fabs(ritz[k]) )
       ++nsig;
   }
   write_header(fpo2, fp, ascii, nrow, ncol, nsig, "las1");

   nsig = 0;
   id = 0;
   id2 = jsq-js;
   for (k = 0; k < js; k++) {
      tmp = id2;
      if (bnd[k] <= kappa * fabs(ritz[k]) ) {
         for (i = 0; i < n; i++) w1[i] = 0.0;
         for (i = 0; i < js; i++) {
	    store(n, RETRQ, i, w2);
	    daxpy(n, s[tmp], w2, 1, w1, 1);
	    tmp -= js;
         }
	 if (ascii) {
	   for (int zz = 0; zz < size/sizeof(double); zz++) {
	     fprintf(fp, "%f ", w1[zz]);
	   }
	   fprintf(fp, "\n");
	 } 
	 else 
	   write(fpo2, (char *)w1, size);

         /* store the w1 vector sequentially in array xv1 */
         for (i = 0; i < n; i++) xv1[id++] = w1[i];
         nsig += 1; 
      }
      id2++;
   }
   free(s);
   return;
}


/***********************************************************************
 *                                                                     *
 *                     store()                                         *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   store() is a user-supplied function which, based on the input
   operation flag, stores to or retrieves from memory a vector.


   Arguments 
   ---------

   (input)
   n       length of vector to be stored or retrieved
   isw     operation flag:
	     isw = 1 request to store j-th Lanczos vector q(j)
	     isw = 2 request to retrieve j-th Lanczos vector q(j)
	     isw = 3 request to store q(j) for j = 0 or 1
	     isw = 4 request to retrieve q(j) for j = 0 or 1
   s	   contains the vector to be stored for a "store" request 

   (output)
   s	   contains the vector retrieved for a "retrieve" request 

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

   BLAS		dcopy

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

void s_las1::store(long n, long isw, long j, double *s)

{
   switch(isw) {
   case STORQ:	dcopy(n, s, 1, &a[(j+MAXLL) * n], 1);
		break;
   case RETRQ:	dcopy(n, &a[(j+MAXLL) * n], 1, s, 1);
		break;
   case STORP:	if (j >= MAXLL) {
		   fprintf(stderr,"store: (STORP) j >= MAXLL \n");
		   break;
		}
		dcopy(n, s, 1, &a[j*n], 1);
		break;
   case RETRP:	if (j >= MAXLL) {
		   fprintf(stderr,"store: (RETRP) j >= MAXLL \n");
		   break;
		}
   		dcopy(n, &a[j*n], 1, s, 1);
		break;
   }
   return;
}



/***********************************************************************
 *                                                                     *
 *                          lanso()                                    *
 *                                                                     *
 ***********************************************************************/
/***********************************************************************

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

   Function determines when the restart of the Lanczos algorithm should 
   occur and when it should terminate.

   Arguments 
   ---------

   (input)
   n         dimension of the eigenproblem for matrix B
   lanmax    upper limit of desired number of lanczos steps           
   maxprs    upper limit of desired number of eigenpairs             
   endl      left end of interval containing unwanted eigenvalues
   endr      right end of interval containing unwanted eigenvalues
   ritz      array to hold the ritz values                       
   bnd       array to hold the error bounds                          
   wptr      array of pointers that point to work space:            
  	       wptr[0]-wptr[5]  six vectors of length n		
  	       wptr[6] array to hold diagonal of the tridiagonal matrix T
  	       wptr[9] array to hold off-diagonal of T	
  	       wptr[7] orthogonality estimate of Lanczos vectors at 
		 step j
 	       wptr[8] orthogonality estimate of Lanczos vectors at 
		 step j-1

   (output)
   j         number of Lanczos steps actually taken
   neig      number of ritz values stabilized
   ritz      array to hold the ritz values
   bnd       array to hold the error bounds
   ierr      (globally declared) error flag
	     ierr = 8192 if stpone() fails to find a starting vector
	     ierr = k if convergence did not occur for k-th eigenvalue
		    in imtqlb()
	     ierr = 0 otherwise


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

   LAS		stpone, error_bound, lanczos_step
   MISC		dsort2
   UTILITY	imin, imax

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

void s_las1::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;
}