rtresgev

 RTRESGEV   Find a few eigenpairs of a generalized Hermitian definite eigenvalue problem 
   using the Riemannian Trust-Region method.

   Find a few smallest generalized eigenvalues and eigenvectors of the matrices A and B.
   Uses an implentation of the Riemannian Trust-Region algorithm.

   D = RTRESGEV(A) returns a vector of A's 6 smallest signed eigenvalues.
   A must be Hermitian and should be large and sparse.

   [V,D] = RTRESGEV(A) returns a diagonal matrix D of A's 6 smallest signed
   eigenvalues and a matrix V whose columns are the corresponding eigenvectors.

   [V,D,OPS] = RTRESGEV(A) also returns a structure containing information
   on the performance of the algorithm.

   RTRESGEV(A,B) solves the generalized eigenvalue problem A*V == B*V*D.  
   B must be Hermitian positive definite and the same size as A.
   The positive definiteness of B is taken for granted, and will not
   be tested.
   RTRESGEV(A,[],...) indicates the standard eigenvalue problem A*V == V*D.

   RTRESGEV(A,K) and RTRESGEV(A,B,K) return the K smallest signed eigenvalues.

   RTRESGEV(A,K,OPTS) and RTRESGEV(A,B,K,OPTS) specify options:
     OPTS.verbosity: verbose printing at each outer step [0 | {1} | 2]
     OPTS.x0: basis for the initial iterate [{randomly generated}]
     OPTS.tol: outer convergence tolerance [scalar | {1e-6}]
     OPTS.outerstop: mechanism for outer convergence 
          'gradabs' absolute tolerance: quit when |grad_i| < tol
          'gradrel' relative tolerance: quit when |grad_i|/|grad_0| < tol  
          'resrel'  relative tolerance: quit when |res_i|/|lambda_i| < tol (default)

   RTRESGEV(AFUN,N,K,OPTS) and RTRESGEV(AFUN,N,BFUN,K,OPTS) accepts the
   functions AFUN and BFUN instead of the matrix A, where N is the size of A.
   The functions are called like:
      y = feval(AFUN,x);
   and
      y = feval(BFUN,x);

     ---------------------- ADVANCED OPTIONS -----------------------

     Note, any field in OPTS that is not listed below will be passed  
     out of RTRESGEV in OPS. For example,                             
       >> options.desc = 'my favorite test';                            
       >> [V,L,O] = RTRESGEV(...,options);                              
       >> O.desc
       ans = my favorite test

     OPTS.fout: file identifier for printing {1}
          All print statements are directed to fout.
          Typical choices are:
             1, standard output,
             2, standard error,
             or a file handle returned by fopen.
          RTRESGEV will not close the file handle.
     OPTS.debug: debugging checks and output [{0} | 1]
     OPTS.useSA: use 2-D subspace acceleration [false]
                 This is only in effect early in the algorithm.
     OPTS.p: block size for problem [integer | {k}]
          Must satisfy k <= OPTS.p <= n
     OPTS.paramtest: parameter test          [{0} | 1]
          Do not actually run the algorithm, but instead
          test the given parameters, set the defaults, and return them
          in a struct.
          Example: 
             defaultopts = rtresgev(A,B,k,opts);
     OPTS.Vtest: space against which to test at each step
          Canonical sines are computed w.r.t. the euclidean inner product and 
          are stored in OPS.sines. (See note in code header for methodology.)
     OPTS.Ltest: values against which to test at each step 
          Errors between the values in Ltest and the current eigenvalue 
          estimates are stored in OPS.verrors. OPTS.Ltest must contain
          p values, stored in a vector or a diagonal matrix.
     OPTS.maxinner: max number of inner iterations
          Default is the dimension of the tangent space:
             d = n*p - (p^2+p)/2
     OPTS.innerstop: available stopping criteria for the model minimization
          A bitmask of three bits, the default is 3=1+2=maxinner or RTR model grad-based
             bitand(innerstop,1): maxinner iterations reached
             bitand(innerstop,2): RTR model gradient-based stop is enabled (kappa/theta)
             bitand(innerstop,4): stop inner when outer stopping criteria is satisfied
          At least one of these bits must be set.
     OPTS.maxouter: max number of outer iterations  [integer | {100}]
     OPTS.minouter: min number of outer iterations  [integer]
          When OPTS.randomize == 1, default is 2
          Otherwise, default is 0.
          Useful when OPTS.randomize == 1, to allow the algorithm
          opportunity to escape a non-optimal critical point.
          Note, the algorithm always performs one outer step before
          checking the outer stopping criteria, as long as the gradient
          is non-zero.
     OPTS.randomize: randomize model minimization   [{0} | 1]
          This randomizes the initial solution to the model minimization.
          Otherwise, initial solution is the origin. See OPTS.min_outer.
          Randomization must be disabled if a preconditioner is used.
     OPTS.kappa: convergence: see paper   [scalar | {0.1}]
          kappa effects a linear reduction in the 
          residual before truncating the inner solver.
          Must satisfy 0 < kappa <= 1
     OPTS.theta: convergence: see paper   [scalar | {1.0}]
          theta effects a superlinear convergence rate
          by introducing a more strict stopping criteria on 
          the inner problem, as the algorithm nears solution.
          The algorithm seeks order theta+1 convergence.
          Must satisfy 0 <= theta <= 2 for a maximum cubic 
          rate of convergence
     OPTS.rho_prime: TR: acceptance ratio [scalar | {0.1}]
          The ratio between the reduction in the Rayleigh
          quotient and reduction in the quadratic model
          is used as an acceptance criteria for a new proposal.
          Should satisfy 0 < rho_prime < .25
     OPTS.Delta_bar: TR: max trust region size [scalar | {inf}]
     OPTS.Delta0: TR: initial trust region size [scalar | {p*sqrt(3)}]
          Must satisfy 0 <= Delta0 <= Delta_bar
     OPTS.Prec: preconditioner for the inner iteration.

     Regarding Preconditioning: A preconditioner for the inner iteration may be 
        specified via OPTS.Prec. The preconditioner will be invoked as follows:
           Peta = feval(OPTS.Prec,X,BX,eta);
        The input to the preconditioner will always be B-orthogonal.
        The preconditioner must return a matrix which is B-orthogonal to X.
        See, for example, 'precDLP'

    Examples:
       A = delsq(numgrid('C',15));
       B = speye(15);
       [V,D] = rtresgev(A,B,5);

 See also IRTRESGEV