Documentation updates

Improved consistency of descriptions and better descriptions of stop_fn.

GKD.m

@@ -1,5 +1,5 @@
 function [U,S,V,HIST,UD] = GKD(A,numVals,varargin)
-%% GKD   Find a few singular values and vectors for large sparse matrices
+%% GKD   Find a few singular triplets for large sparse matrices
 %
 % [U,S,V,H,UD] = GKD(A,numVals) returns numVals singular triplets
 %       closest to target with convergence history
@@ -8,46 +8,95 @@
 % with orthonormal columns, S is K-by-K diagonal, and V is N-by-K with
 % orthonormal columns.
 %
+% HIST returns the total time, number of iterations, matrix vector products,
+% number of restarts, number of converged values (residual < normA*tol), and
+% the smallest unconverged residual.
+%
+% UD returns any userdata needed or generated by custom stopping functions.
+% The standard stopping criteria for GKD do not utilize this structure. 
+%
+% If A is a function, A(x,'notransp') must return the matrix vector
+% product A*x and A(x,'transp') must return A'*x. Additionally, the 
+% optional name value pairs 'm' and 'n' must be configured with the 
+% number of rows and columns of A respectively.
+%
 % [U,S,V,H,UD] = GKD(A,numVals,...) configures additional options
 %
 % Options for GKD (Name, Value pairs)
 % tol           residual norm tolerance
 % SIGMA         Sets which values are desired ('L' or 'S')
-% target_fn  Function used for expanding the basis
-%    [index,userdata] = target_fn(solverdata,userdata)
-% stop_fn     Function used for stopping the solver
-%    [done,numVals,userdata] = stop_fn(numVals,solverdata,userdata)
+% target_fn     Function used for expanding the basis ('resid' or 'prog_tol')
+% stop_fn       Function used for stopping the solver
+%                   [done,numVals,userdata] = stop_fn(numVals,solverdata,userdata)
 % maxMV         maximum number of matrix-vector multiplications
+% maxTime       maximum solver time
 % normA         norm(A) estimate
-% display      Prints partial history to console if set
-% v0            Initial vector for V
-% b             Block size
+% display       Prints partial history to console if set
+% v0            initial vectors with min(size(A)) rows
+% b             block size
 % minRestart    number of vectors to maintain after restart
 % maxBasis      max number of basis vectors in V,U
-% numOld        number of +k vectors to keep during restart
-% maxQMR         max number of inner solver iterations
+% numOld        number of +k vectors to keep during restart (must be <= b)
+% maxQMR        max number of QMR iterations
 % seed          random seed
 % m             number of rows in A (if A is a function_handle)
 % n             number of cols in A (if A is a function_handle)
 % P             preconditioner for AtA
+% userdata      external data needed for custom stopping functions
 %
 % Default Options Settings
 %
 % tol          1e-6
 % SIGMA        'L' (Largest)
-% target_fn    Targeting based on first values with residuals above 'tol'
-% stop_fn      Stopping based on residual tolerance ('tol')
-% maxMV        -1 (No stopping based on matvecs)
+% target_fn    'resid'
+% stop_fn      'resid' (Stops when the first numVals residuals are < tol)
+% maxMV        inf (No stopping based on matvecs)
+% maxTime      inf (No time based stopping)
 % normA        Largest value seen (Accurate when SIGMA = 'L')
 % display      0 (Off)
 % v0           Gaussian random vectors
 % b            1
-% minRestart        numVals+max(b,15)
-% maxBasis     max(minRestart+2*b,floor(1.3*minRestart))
-% numOld       1
+% minRestart   max(7,p.numVals+5)
+% maxBasis     max([15,p.minRestart+4*p.b,floor(1.3*p.minRestart)])
+% numOld       -1 (Uses current block size for +k)
 % maxQMR       0
 % seed         'shuffle' (sets rng based on current time)
-% P            1 (Identity matrix)
+% P            []
+% userdata     []
+%
+% Using the stop_fn option
+% 
+% Requires a funtion_handle with the following form:
+%       [done,numVals,userdata] = stop_fn(numVals,solverdata,userdata)
+% Inputs:
+% numVals      The number of singular values desired. This is given as an 
+%              input so the user function does not need to assign the
+%              numVals output variable unless desired.
+% solverdata   Includes the following information: 
+%                   mvp: number of matrix vector products,
+%                   outerits: number of outer iterations,
+%                   time: total time elapsed,
+%                   normA: ||A||,
+%                   k: current basis size,
+%                   s: a (k x 1) vector of ritz values,
+%                   b: current block size,
+%                   tol: minimum residual tolerance,
+%                   rn: a (maxBasis x 1) vector of the most recent residual norms
+% userdata     Inputs the userdata structure given as an optional parameter. 
+%              This can be updated inside of the stopping function as long
+%              as the updates are returned in the third return argument
+%
+% Outputs:
+% done         Convergence flag. If set, the solver will stop and return
+%              the most recent numVals Ritz vectors and values. 
+% numVals      Determines the number of values to return when the 'done'
+%              flag is set. Can be set separately, but should always be 
+%              less than or equal to the numVals input. Doing this without 
+%              setting the 'done' flag may change targeting to improve
+%              convergence.
+% userdata     Allows for updates to the userdata structure to be stored
+%              for later iterations.
+
 
 
 %% Input Parsing %%
@@ -72,7 +121,7 @@
 addParameter(p,'m',-1);
 addParameter(p,'n',-1);
 addParameter(p,'P',[]);
-addParameter(p,'user_data',[]);
+addParameter(p,'userdata',[]);
 
 parse(p,A,numVals,varargin{:});
 
@@ -181,7 +230,7 @@
         'numVals',p.numVals,'normA',normA,'k',k,'s',sr, ...
         'b',p.b,'tol',p.tol,'rn',allrun);
 
-    [index,p.user_data] = p.target_fn(solver_data,p.user_data);
+    [index,p.userdata] = p.target_fn(solver_data,p.userdata);
 
     u = U(:,1:k)*ur(:,index);
     v = V(:,1:k)*vr(:,index);
@@ -200,16 +249,15 @@
         length(find(allrun < normA*p.tol)), min(allrun(allrun > normA*p.tol))/normA];
 
     solver_data.rn = allrun;
-    [done,p.numVals,p.user_data] = p.stop_fn(p.numVals,solver_data,p.user_data);
+    [done,p.numVals,p.userdata] = p.stop_fn(p.numVals,solver_data,p.userdata);
 
     if isempty(index) || done
         U = U(:,1:k)*ur(:,1:min(p.numVals,k));
         V = V(:,1:k)*vr(:,1:min(p.numVals,k));
-        S = sr(1:min(p.numVals,k));
+        S = diag(sr(1:min(p.numVals,k)));
         break;
     end
 
-
     cb_size = size(ru,2);
     if p.maxQMR > 0 || ~isempty(p.P)
         for j = 1:cb_size
@@ -276,7 +324,7 @@
     V = Temp;
 end
 
-UD = p.user_data;
+UD = p.userdata;
 
 end