octave_packages/linear-algebra-2.2.0/lobpcg.m

   1 ## Copyright (C) 2000-2011 A.V. Knyazev <Andrew.Knyazev@ucdenver.edu>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU Lesser General Public License as published by the Free
   5 ## Software Foundation; either version 3 of the License, or (at your option) any
   6 ## later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
  11 ## for more details.
  12 ##
  13 ## You should have received a copy of the GNU Lesser General Public License
  14 ## along with this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {[@var{blockVectorX}, @var{lambda}] =} lobpcg (@var{blockVectorX}, @var{operatorA})
  18 ## @deftypefnx {Function File} {[@var{blockVectorX}, @var{lambda}, @var{failureFlag}] =} lobpcg (@var{blockVectorX}, @var{operatorA})
  19 ## @deftypefnx {Function File} {[@var{blockVectorX}, @var{lambda}, @var{failureFlag}, @var{lambdaHistory}, @var{residualNormsHistory}] =} lobpcg (@var{blockVectorX}, @var{operatorA}, @var{operatorB}, @var{operatorT}, @var{blockVectorY}, @var{residualTolerance}, @var{maxIterations}, @var{verbosityLevel})
  20 ## Solves Hermitian partial eigenproblems using preconditioning.
  21 ##
  22 ## The first form outputs the array of algebraic smallest eigenvalues @var{lambda} and
  23 ## corresponding matrix of orthonormalized eigenvectors @var{blockVectorX} of the
  24 ## Hermitian (full or sparse) operator @var{operatorA} using input matrix
  25 ## @var{blockVectorX} as an initial guess, without preconditioning, somewhat
  26 ## similar to:
  27 ##
  28 ## @example
  29 ## # for real symmetric operator operatorA
  30 ## opts.issym  = 1; opts.isreal = 1; K = size (blockVectorX, 2);
  31 ## [blockVectorX, lambda] = eigs (operatorA, K, 'SR', opts);
  32 ##
  33 ## # for Hermitian operator operatorA
  34 ## K = size (blockVectorX, 2);
  35 ## [blockVectorX, lambda] = eigs (operatorA, K, 'SR');
  36 ## @end example
  37 ##
  38 ## The second form returns a convergence flag. If @var{failureFlag} is 0 then
  39 ## all the eigenvalues converged; otherwise not all converged.
  40 ##
  41 ## The third form computes smallest eigenvalues @var{lambda} and corresponding eigenvectors
  42 ## @var{blockVectorX} of the generalized eigenproblem Ax=lambda Bx, where
  43 ## Hermitian operators @var{operatorA} and @var{operatorB} are given as functions, as
  44 ## well as a preconditioner, @var{operatorT}. The operators @var{operatorB} and
  45 ## @var{operatorT} must be in addition @emph{positive definite}. To compute the largest
  46 ## eigenpairs of @var{operatorA}, simply apply the code to @var{operatorA} multiplied by
  47 ## -1. The code does not involve @emph{any} matrix factorizations of @var{operatorA} and
  48 ## @var{operatorB}, thus, e.g., it preserves the sparsity and the structure of
  49 ## @var{operatorA} and @var{operatorB}.
  50 ##
  51 ## @var{residualTolerance} and @var{maxIterations} control tolerance and max number of
  52 ## steps, and @var{verbosityLevel} = 0, 1, or 2 controls the amount of printed
  53 ## info. @var{lambdaHistory} is a matrix with all iterative lambdas, and
  54 ## @var{residualNormsHistory} are matrices of the history of 2-norms of residuals
  55 ##
  56 ## Required input:
  57 ## @itemize @bullet
  58 ## @item
  59 ## @var{blockVectorX} (class numeric) - initial approximation to eigenvectors,
  60 ## full or sparse matrix n-by-blockSize. @var{blockVectorX} must be full rank.
  61 ## @item
  62 ## @var{operatorA} (class numeric, char, or function_handle) - the main operator
  63 ## of the eigenproblem, can be a matrix, a function name, or handle
  64 ## @end itemize
  65 ##
  66 ## Optional function input:
  67 ## @itemize @bullet
  68 ## @item
  69 ## @var{operatorB} (class numeric, char, or function_handle) - the second operator,
  70 ## if solving a generalized eigenproblem, can be a matrix, a function name, or
  71 ## handle; by default if empty, @code{operatorB = I}.
  72 ## @item
  73 ## @var{operatorT} (class char or function_handle) - the preconditioner, by
  74 ## default @code{operatorT(blockVectorX) = blockVectorX}.
  75 ## @end itemize
  76 ##
  77 ## Optional constraints input:
  78 ## @itemize @bullet
  79 ## @item
  80 ## @var{blockVectorY} (class numeric) - a full or sparse n-by-sizeY matrix of
  81 ## constraints, where sizeY < n. @var{blockVectorY} must be full rank. The
  82 ## iterations will be performed in the (operatorB-) orthogonal complement of the
  83 ## column-space of @var{blockVectorY}.
  84 ## @end itemize
  85 ##
  86 ## Optional scalar input parameters:
  87 ## @itemize @bullet
  88 ## @item
  89 ## @var{residualTolerance} (class numeric) - tolerance, by default, @code{residualTolerance = n * sqrt (eps)}
  90 ## @item
  91 ## @var{maxIterations} - max number of iterations, by default, @code{maxIterations = min (n, 20)}
  92 ## @item
  93 ## @var{verbosityLevel} - either 0 (no info), 1, or 2 (with pictures); by
  94 ## default, @code{verbosityLevel = 0}.
  95 ## @end itemize
  96 ##
  97 ## Required output:
  98 ## @itemize @bullet
  99 ## @item
 100 ## @var{blockVectorX} and @var{lambda} (class numeric) both are computed
 101 ## blockSize eigenpairs, where @code{blockSize = size (blockVectorX, 2)}
 102 ## for the initial guess @var{blockVectorX} if it is full rank.
 103 ## @end itemize
 104 ##
 105 ## Optional output:
 106 ## @itemize @bullet
 107 ## @item
 108 ## @var{failureFlag} (class integer) as described above.
 109 ## @item
 110 ## @var{lambdaHistory} (class numeric) as described above.
 111 ## @item
 112 ## @var{residualNormsHistory} (class numeric) as described above.
 113 ## @end itemize
 114 ##
 115 ## Functions @code{operatorA(blockVectorX)}, @code{operatorB(blockVectorX)} and
 116 ## @code{operatorT(blockVectorX)} must support @var{blockVectorX} being a matrix, not
 117 ## just a column vector.
 118 ##
 119 ## Every iteration involves one application of @var{operatorA} and @var{operatorB}, and
 120 ## one of @var{operatorT}.
 121 ##
 122 ## Main memory requirements: 6 (9 if @code{isempty(operatorB)=0}) matrices of the
 123 ## same size as @var{blockVectorX}, 2 matrices of the same size as @var{blockVectorY}
 124 ## (if present), and two square matrices of the size 3*blockSize.
 125 ##
 126 ## In all examples below, we use the Laplacian operator in a 20x20 square
 127 ## with the mesh size 1 which can be generated in MATLAB by running:
 128 ## @example
 129 ## A = delsq (numgrid ('S', 21));
 130 ## n = size (A, 1);
 131 ## @end example
 132 ##
 133 ## or in MATLAB and Octave by:
 134 ## @example
 135 ## [~,~,A] = laplacian ([19, 19]);
 136 ## n = size (A, 1);
 137 ## @end example
 138 ##
 139 ## Note that @code{laplacian} is a function of the specfun octave-forge package.
 140 ##
 141 ## The following Example:
 142 ## @example
 143 ## [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, 1e-5, 50, 2);
 144 ## @end example
 145 ##
 146 ## attempts to compute 8 first eigenpairs without preconditioning, but not all
 147 ## eigenpairs converge after 50 steps, so failureFlag=1.
 148 ##
 149 ## The next Example:
 150 ## @example
 151 ## blockVectorY = [];
 152 ## lambda_all = [];
 153 ## for j = 1:4
 154 ##   [blockVectorX, lambda] = lobpcg (randn (n, 2), A, blockVectorY, 1e-5, 200, 2);
 155 ##   blockVectorY           = [blockVectorY, blockVectorX];
 156 ##   lambda_all             = [lambda_all' lambda']';
 157 ##   pause;
 158 ## end
 159 ## @end example
 160 ##
 161 ## attemps to compute the same 8 eigenpairs by calling the code 4 times with
 162 ## blockSize=2 using orthogonalization to the previously founded eigenvectors.
 163 ##
 164 ## The following Example:
 165 ## @example
 166 ## R       = ichol (A, struct('michol', 'on'));
 167 ## precfun = @@(x)R\(R'\x);
 168 ## [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, [], @@(x)precfun(x), 1e-5, 60, 2);
 169 ## @end example
 170 ##
 171 ## computes the same eigenpairs in less then 25 steps, so that failureFlag=0
 172 ## using the preconditioner function @code{precfun}, defined inline. If @code{precfun}
 173 ## is defined as an octave function in a file, the function handle
 174 ## @code{@@(x)precfun(x)} can be equivalently replaced by the function name @code{precfun}. Running:
 175 ##
 176 ## @example
 177 ## [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, speye (n), @@(x)precfun(x), 1e-5, 50, 2);
 178 ## @end example
 179 ##
 180 ## produces similar answers, but is somewhat slower and needs more memory as
 181 ## technically a generalized eigenproblem with B=I is solved here.
 182 ##
 183 ## The following example for a mostly diagonally dominant sparse matrix A
 184 ## demonstrates different types of preconditioning, compared to the standard
 185 ## use of the main diagonal of A:
 186 ##
 187 ## @example
 188 ## clear all; close all;
 189 ## n       = 1000;
 190 ## M       = spdiags ([1:n]', 0, n, n);
 191 ## precfun = @@(x)M\x;
 192 ## A       = M + sprandsym (n, .1);
 193 ## Xini    = randn (n, 5);
 194 ## maxiter = 15;
 195 ## tol     = 1e-5;
 196 ## [~,~,~,~,rnp] = lobpcg (Xini, A, tol, maxiter, 1);
 197 ## [~,~,~,~,r]   = lobpcg (Xini, A, [], @@(x)precfun(x), tol, maxiter, 1);
 198 ## subplot (2,2,1), semilogy (r'); hold on;
 199 ## semilogy (rnp', ':>');
 200 ## title ('No preconditioning (top)'); axis tight;
 201 ## M(1,2)  = 2;
 202 ## precfun = @@(x)M\x; % M is no longer symmetric
 203 ## [~,~,~,~,rns] = lobpcg (Xini, A, [], @@(x)precfun(x), tol, maxiter, 1);
 204 ## subplot (2,2,2), semilogy (r'); hold on;
 205 ## semilogy (rns', '--s');
 206 ## title ('Nonsymmetric preconditioning (square)'); axis tight;
 207 ## M(1,2)  = 0;
 208 ## precfun = @@(x)M\(x+10*sin(x)); % nonlinear preconditioning
 209 ## [~,~,~,~,rnl] = lobpcg (Xini, A, [], @@(x)precfun(x), tol, maxiter, 1);
 210 ## subplot (2,2,3),  semilogy (r'); hold on;
 211 ## semilogy (rnl', '-.*');
 212 ## title ('Nonlinear preconditioning (star)'); axis tight;
 213 ## M       = abs (M - 3.5 * speye (n, n));
 214 ## precfun = @@(x)M\x;
 215 ## [~,~,~,~,rs] = lobpcg (Xini, A, [], @@(x)precfun(x), tol, maxiter, 1);
 216 ## subplot (2,2,4),  semilogy (r'); hold on;
 217 ## semilogy (rs', '-d');
 218 ## title ('Selective preconditioning (diamond)'); axis tight;
 219 ## @end example
 220 ##
 221 ## @heading References
 222 ## This main function @code{lobpcg} is a version of the preconditioned conjugate
 223 ## gradient method (Algorithm 5.1) described in A. V. Knyazev, Toward the Optimal
 224 ## Preconditioned Eigensolver:
 225 ## Locally Optimal Block Preconditioned Conjugate Gradient Method,
 226 ## SIAM Journal on Scientific Computing 23 (2001), no. 2, pp. 517-541.
 227 ## @uref{http://dx.doi.org/10.1137/S1064827500366124}
 228 ##
 229 ## @heading Known bugs/features
 230 ## @itemize @bullet
 231 ## @item
 232 ## an excessively small requested tolerance may result in often restarts and
 233 ## instability. The code is not written to produce an eps-level accuracy! Use
 234 ## common sense.
 235 ## @item
 236 ## the code may be very sensitive to the number of eigenpairs computed,
 237 ## if there is a cluster of eigenvalues not completely included, cf.
 238 ## @example
 239 ## operatorA = diag ([1 1.99 2:99]);
 240 ## [blockVectorX, lambda] = lobpcg (randn (100, 1),operatorA, 1e-10, 80, 2);
 241 ## [blockVectorX, lambda] = lobpcg (randn (100, 2),operatorA, 1e-10, 80, 2);
 242 ## [blockVectorX, lambda] = lobpcg (randn (100, 3),operatorA, 1e-10, 80, 2);
 243 ## @end example
 244 ## @end itemize
 245 ##
 246 ## @heading Distribution
 247 ## The main distribution site: @uref{http://math.ucdenver.edu/~aknyazev/}
 248 ##
 249 ## A C-version of this code is a part of the @uref{http://code.google.com/p/blopex/}
 250 ## package and is directly available, e.g., in PETSc and HYPRE.
 251 ## @end deftypefn
 252
 253 function [blockVectorX,lambda,varargout] = lobpcg(blockVectorX,operatorA,varargin)
 254   %Begin
 255
 256   % constants
 257
 258   CONVENTIONAL_CONSTRAINTS = 1;
 259   SYMMETRIC_CONSTRAINTS = 2;
 260
 261   %Initial settings
 262
 263   failureFlag = 1;
 264   if nargin < 2
 265       error('BLOPEX:lobpcg:NotEnoughInputs',...
 266       strcat('There must be at least 2 input agruments: ',...
 267           'blockVectorX and operatorA'));
 268   end
 269   if nargin > 8
 270       warning('BLOPEX:lobpcg:TooManyInputs',...
 271           strcat('There must be at most 8 input agruments ',...
 272           'unless arguments are passed to a function'));
 273   end
 274
 275   if ~isnumeric(blockVectorX)
 276       error('BLOPEX:lobpcg:FirstInputNotNumeric',...
 277           'The first input argument blockVectorX must be numeric');
 278   end
 279   [n,blockSize]=size(blockVectorX);
 280   if blockSize > n
 281       error('BLOPEX:lobpcg:FirstInputFat',...
 282       'The first input argument blockVectorX must be tall, not fat');
 283   end
 284   if n < 6
 285       error('BLOPEX:lobpcg:MatrixTooSmall',...
 286           'The code does not work for matrices of small sizes');
 287   end
 288
 289   if isa(operatorA,'numeric')
 290       nA = size(operatorA,1);
 291       if any(size(operatorA) ~= nA)
 292           error('BLOPEX:lobpcg:MatrixNotSquare',...
 293               'operatorA must be a square matrix or a string');
 294       end
 295       if size(operatorA) ~= n
 296           error('BLOPEX:lobpcg:MatrixWrongSize',...
 297           ['The size ' int2str(size(operatorA))...
 298               ' of operatorA is not the same as ' int2str(n)...
 299               ' - the number of rows of blockVectorX']);
 300       end
 301   end
 302
 303   count_string = 0;
 304
 305   operatorT = [];
 306   operatorB = [];
 307   residualTolerance = [];
 308   maxIterations = [];
 309   verbosityLevel = [];
 310   blockVectorY = []; sizeY = 0;
 311   for j = 1:nargin-2
 312       if isequal(size(varargin{j}),[n,n])
 313           if isempty(operatorB)
 314               operatorB = varargin{j};
 315           else
 316               error('BLOPEX:lobpcg:TooManyMatrixInputs',...
 317           strcat('Too many matrix input arguments. ',...
 318           'Preconditioner operatorT must be an M-function'));
 319           end
 320       elseif isequal(size(varargin{j},1),n) && size(varargin{j},2) < n
 321           if isempty(blockVectorY)
 322               blockVectorY = varargin{j};
 323               sizeY=size(blockVectorY,2);
 324           else
 325               error('BLOPEX:lobpcg:WrongConstraintsFormat',...
 326               'Something wrong with blockVectorY input argument');
 327           end
 328       elseif ischar(varargin{j}) || isa(varargin{j},'function_handle')
 329           if count_string == 0
 330               if isempty(operatorB)
 331                   operatorB = varargin{j};
 332                   count_string = count_string + 1;
 333               else
 334                   operatorT = varargin{j};
 335               end
 336           elseif count_string == 1
 337               operatorT = varargin{j};
 338           else
 339               warning('BLOPEX:lobpcg:TooManyStringFunctionHandleInputs',...
 340                   'Too many string or FunctionHandle input arguments');
 341           end
 342       elseif isequal(size(varargin{j}),[n,n])
 343           error('BLOPEX:lobpcg:WrongPreconditionerFormat',...
 344           'Preconditioner operatorT must be an M-function');
 345       elseif max(size(varargin{j})) == 1
 346           if isempty(residualTolerance)
 347               residualTolerance = varargin{j};
 348           elseif isempty(maxIterations)
 349               maxIterations = varargin{j};
 350           elseif isempty(verbosityLevel)
 351               verbosityLevel = varargin{j};
 352           else
 353               warning('BLOPEX:lobpcg:TooManyScalarInputs',...
 354                   'Too many scalar parameters, need only three');
 355           end
 356       elseif isempty(varargin{j})
 357           if isempty(operatorB)
 358               count_string = count_string + 1;
 359           elseif ~isempty(operatorT)
 360               count_string = count_string + 1;
 361           elseif ~isempty(blockVectorY)
 362               error('BLOPEX:lobpcg:UnrecognizedEmptyInput',...
 363                  ['Unrecognized empty input argument number ' int2str(j+2)]);
 364           end
 365       else
 366           error('BLOPEX:lobpcg:UnrecognizedInput',...
 367               ['Input argument number ' int2str(j+2) ' not recognized.']);
 368       end
 369   end
 370
 371   if verbosityLevel
 372       if issparse(blockVectorX)
 373           fprintf(['The sparse initial guess with %i colunms '...
 374           'and %i raws is detected  \n'],n,blockSize);
 375       else
 376           fprintf(['The full initial guess with %i colunms '...
 377               'and %i raws is detected  \n'],n,blockSize);
 378       end
 379       if ischar(operatorA)
 380           fprintf('The main operator is detected as an M-function %s \n',...
 381               operatorA);
 382       elseif isa(operatorA,'function_handle')
 383           fprintf('The main operator is detected as an M-function %s \n',...
 384               func2str(operatorA));
 385       elseif issparse(operatorA)
 386           fprintf('The main operator is detected as a sparse matrix \n');
 387       else
 388           fprintf('The main operator is detected as a full matrix \n');
 389       end
 390       if isempty(operatorB)
 391           fprintf('Solving standard eigenvalue problem, not generalized \n');
 392       elseif ischar(operatorB)
 393           fprintf(['The second operator of the generalized eigenproblem \n'...
 394           'is detected as an M-function %s \n'],operatorB);
 395       elseif isa(operatorB,'function_handle')
 396           fprintf(['The second operator of the generalized eigenproblem \n'...
 397           'is detected as an M-function %s \n'],func2str(operatorB));
 398       elseif issparse(operatorB)
 399           fprintf(strcat('The second operator of the generalized',...
 400               'eigenproblem \n is detected as a sparse matrix \n'));
 401       else
 402           fprintf(strcat('The second operator of the generalized',...
 403               'eigenproblem \n is detected as a full matrix \n'));
 404       end
 405       if isempty(operatorT)
 406           fprintf('No preconditioner is detected \n');
 407       elseif ischar(operatorT)
 408           fprintf('The preconditioner is detected as an M-function %s \n',...
 409               operatorT);
 410       elseif isa(operatorT,'function_handle')
 411           fprintf('The preconditioner is detected as an M-function %s \n',...
 412               func2str(operatorT));
 413       end
 414       if isempty(blockVectorY)
 415           fprintf('No matrix of constraints is detected \n')
 416       elseif issparse(blockVectorY)
 417           fprintf('The sparse matrix of %i constraints is detected \n',sizeY);
 418       else
 419           fprintf('The full matrix of %i constraints is detected \n',sizeY);
 420       end
 421       if issparse(blockVectorY) ~= issparse(blockVectorX)
 422           warning('BLOPEX:lobpcg:SparsityInconsistent',...
 423               strcat('The sparsity formats of the initial guess and ',...
 424               'the constraints are inconsistent'));
 425       end
 426   end
 427
 428   % Set defaults
 429
 430   if isempty(residualTolerance)
 431       residualTolerance = sqrt(eps)*n;
 432   end
 433   if isempty(maxIterations)
 434       maxIterations = min(n,20);
 435   end
 436   if isempty(verbosityLevel)
 437       verbosityLevel = 0;
 438   end
 439
 440   if verbosityLevel
 441       fprintf('Tolerance %e and maximum number of iterations %i \n',...
 442           residualTolerance,maxIterations)
 443   end
 444
 445   %constraints preprocessing
 446   if isempty(blockVectorY)
 447       constraintStyle = 0;
 448   else
 449       %    constraintStyle = SYMMETRIC_CONSTRAINTS; % more accurate?
 450       constraintStyle = CONVENTIONAL_CONSTRAINTS;
 451   end
 452
 453   if constraintStyle == CONVENTIONAL_CONSTRAINTS
 454
 455       if isempty(operatorB)
 456           gramY = blockVectorY'*blockVectorY;
 457       else
 458           if isnumeric(operatorB)
 459               blockVectorBY = operatorB*blockVectorY;
 460           else
 461               blockVectorBY = feval(operatorB,blockVectorY);
 462           end
 463           gramY=blockVectorY'*blockVectorBY;
 464       end
 465       gramY=(gramY'+gramY)*0.5;
 466       if isempty(operatorB)
 467           blockVectorX = blockVectorX - ...
 468               blockVectorY*(gramY\(blockVectorY'*blockVectorX));
 469       else
 470           blockVectorX =blockVectorX - ...
 471               blockVectorY*(gramY\(blockVectorBY'*blockVectorX));
 472       end
 473
 474   elseif constraintStyle == SYMMETRIC_CONSTRAINTS
 475
 476       if ~isempty(operatorB)
 477           if isnumeric(operatorB)
 478               blockVectorY = operatorB*blockVectorY;
 479           else
 480               blockVectorY = feval(operatorB,blockVectorY);
 481           end
 482       end
 483       if isempty(operatorT)
 484           gramY = blockVectorY'*blockVectorY;
 485       else
 486           blockVectorTY = feval(operatorT,blockVectorY);
 487           gramY = blockVectorY'*blockVectorTY;
 488       end
 489       gramY=(gramY'+gramY)*0.5;
 490       if isempty(operatorT)
 491           blockVectorX = blockVectorX - ...
 492               blockVectorY*(gramY\(blockVectorY'*blockVectorX));
 493       else
 494           blockVectorX = blockVectorX - ...
 495               blockVectorTY*(gramY\(blockVectorY'*blockVectorX));
 496       end
 497
 498   end
 499
 500   %Making the initial vectors (operatorB-) orthonormal
 501   if isempty(operatorB)
 502       %[blockVectorX,gramXBX] = qr(blockVectorX,0);
 503       gramXBX=blockVectorX'*blockVectorX;
 504       if ~isreal(gramXBX)
 505           gramXBX=(gramXBX+gramXBX')*0.5;
 506       end
 507       [gramXBX,cholFlag]=chol(gramXBX);
 508       if  cholFlag ~= 0
 509           error('BLOPEX:lobpcg:ConstraintsTooTight',...
 510              'The initial approximation after constraints is not full rank');
 511       end
 512       blockVectorX = blockVectorX/gramXBX;
 513   else
 514       %[blockVectorX,blockVectorBX] = orth(operatorB,blockVectorX);
 515       if isnumeric(operatorB)
 516           blockVectorBX = operatorB*blockVectorX;
 517       else
 518           blockVectorBX = feval(operatorB,blockVectorX);
 519       end
 520       gramXBX=blockVectorX'*blockVectorBX;
 521       if ~isreal(gramXBX)
 522           gramXBX=(gramXBX+gramXBX')*0.5;
 523       end
 524       [gramXBX,cholFlag]=chol(gramXBX);
 525       if  cholFlag ~= 0
 526           error('BLOPEX:lobpcg:InitialNotFullRank',...
 527               sprintf('%s\n%s', ...
 528               'The initial approximation after constraints is not full rank',...
 529               'or/and operatorB is not positive definite'));
 530       end
 531       blockVectorX = blockVectorX/gramXBX;
 532       blockVectorBX = blockVectorBX/gramXBX;
 533   end
 534
 535   % Checking if the problem is big enough for the algorithm,
 536   % i.e. n-sizeY > 5*blockSize
 537   % Theoretically, the algorithm should be able to run if
 538   % n-sizeY > 3*blockSize,
 539   % but the extreme cases might be unstable, so we use 5 instead of 3 here.
 540   if n-sizeY < 5*blockSize
 541       error('BLOPEX:lobpcg:MatrixTooSmall','%s\n%s', ...
 542       'The problem size is too small, relative to the block size.',...
 543       'Try using eig() or eigs() instead.');
 544   end
 545
 546   % Preallocation
 547   residualNormsHistory=zeros(blockSize,maxIterations);
 548   lambdaHistory=zeros(blockSize,maxIterations+1);
 549   condestGhistory=zeros(1,maxIterations+1);
 550
 551   blockVectorBR=zeros(n,blockSize);
 552   blockVectorAR=zeros(n,blockSize);
 553   blockVectorP=zeros(n,blockSize);
 554   blockVectorAP=zeros(n,blockSize);
 555   blockVectorBP=zeros(n,blockSize);
 556
 557   %Initial settings for the loop
 558   if isnumeric(operatorA)
 559       blockVectorAX = operatorA*blockVectorX;
 560   else
 561       blockVectorAX = feval(operatorA,blockVectorX);
 562   end
 563
 564   gramXAX = full(blockVectorX'*blockVectorAX);
 565   gramXAX = (gramXAX + gramXAX')*0.5;
 566   % eig(...,'chol') uses only the diagonal and upper triangle -
 567   % not true in MATLAB
 568   % Octave v3.2.3-4, eig() does not support inputting 'chol'
 569   [coordX,gramXAX]=eig(gramXAX,eye(blockSize));
 570
 571   lambda=diag(gramXAX); %eig returns non-ordered eigenvalues on the diagonal
 572
 573   if issparse(blockVectorX)
 574       coordX=sparse(coordX);
 575   end
 576
 577   blockVectorX  =  blockVectorX*coordX;
 578   blockVectorAX = blockVectorAX*coordX;
 579   if ~isempty(operatorB)
 580       blockVectorBX = blockVectorBX*coordX;
 581   end
 582   clear coordX
 583
 584   condestGhistory(1)=-log10(eps)/2;  %if too small cause unnecessary restarts
 585
 586   lambdaHistory(1:blockSize,1)=lambda;
 587
 588   activeMask = true(blockSize,1);
 589   % currentBlockSize = blockSize; %iterate all
 590   %
 591   % restart=1;%steepest descent
 592
 593   %The main part of the method is the loop of the CG method: begin
 594   for iterationNumber=1:maxIterations
 595
 596       %     %Computing the active residuals
 597       %     if isempty(operatorB)
 598       %         if currentBlockSize > 1
 599       %             blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
 600       %                 blockVectorX(:,activeMask)*spdiags(lambda(activeMask),0,currentBlockSize,currentBlockSize);
 601       %         else
 602       %             blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
 603       %                 blockVectorX(:,activeMask)*lambda(activeMask);
 604       %                   %to make blockVectorR full when lambda is just a scalar
 605       %         end
 606       %     else
 607       %         if currentBlockSize > 1
 608       %             blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
 609       %                 blockVectorBX(:,activeMask)*spdiags(lambda(activeMask),0,currentBlockSize,currentBlockSize);
 610       %         else
 611       %             blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
 612       %                 blockVectorBX(:,activeMask)*lambda(activeMask);
 613       %                       %to make blockVectorR full when lambda is just a scalar
 614       %         end
 615       %     end
 616
 617       %Computing all residuals
 618       if isempty(operatorB)
 619           if blockSize > 1
 620               blockVectorR = blockVectorAX - ...
 621                   blockVectorX*spdiags(lambda,0,blockSize,blockSize);
 622           else
 623               blockVectorR = blockVectorAX - blockVectorX*lambda;
 624               %to make blockVectorR full when lambda is just a scalar
 625           end
 626       else
 627           if blockSize > 1
 628               blockVectorR=blockVectorAX - ...
 629                   blockVectorBX*spdiags(lambda,0,blockSize,blockSize);
 630           else
 631               blockVectorR = blockVectorAX - blockVectorBX*lambda;
 632               %to make blockVectorR full when lambda is just a scalar
 633           end
 634       end
 635
 636       %Satisfying the constraints for the active residulas
 637       if constraintStyle == SYMMETRIC_CONSTRAINTS
 638           if isempty(operatorT)
 639               blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
 640                   blockVectorY*(gramY\(blockVectorY'*...
 641                   blockVectorR(:,activeMask)));
 642           else
 643               blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
 644                   blockVectorY*(gramY\(blockVectorTY'*...
 645                   blockVectorR(:,activeMask)));
 646           end
 647       end
 648
 649       residualNorms=full(sqrt(sum(conj(blockVectorR).*blockVectorR)'));
 650       residualNormsHistory(1:blockSize,iterationNumber)=residualNorms;
 651
 652       %index antifreeze
 653       activeMask = full(residualNorms > residualTolerance) & activeMask;
 654       %activeMask = full(residualNorms > residualTolerance);
 655       %above allows vectors back into active, which causes problems with frosen Ps
 656       %activeMask = full(residualNorms > 0);      %iterate all, ignore freeze
 657
 658       currentBlockSize = sum(activeMask);
 659       if  currentBlockSize == 0
 660           failureFlag=0; %all eigenpairs converged
 661           break
 662       end
 663
 664       %Applying the preconditioner operatorT to the active residulas
 665       if ~isempty(operatorT)
 666           blockVectorR(:,activeMask) = ...
 667               feval(operatorT,blockVectorR(:,activeMask));
 668       end
 669
 670       if constraintStyle == CONVENTIONAL_CONSTRAINTS
 671           if isempty(operatorB)
 672               blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
 673                   blockVectorY*(gramY\(blockVectorY'*...
 674                   blockVectorR(:,activeMask)));
 675           else
 676               blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
 677                   blockVectorY*(gramY\(blockVectorBY'*...
 678                   blockVectorR(:,activeMask)));
 679           end
 680       end
 681
 682       %Making active (preconditioned) residuals orthogonal to blockVectorX
 683       if isempty(operatorB)
 684           blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
 685               blockVectorX*(blockVectorX'*blockVectorR(:,activeMask));
 686       else
 687           blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
 688               blockVectorX*(blockVectorBX'*blockVectorR(:,activeMask));
 689       end
 690
 691       %Making active residuals orthonormal
 692       if isempty(operatorB)
 693           %[blockVectorR(:,activeMask),gramRBR]=...
 694           %qr(blockVectorR(:,activeMask),0); %to increase stability
 695           gramRBR=blockVectorR(:,activeMask)'*blockVectorR(:,activeMask);
 696           if ~isreal(gramRBR)
 697               gramRBR=(gramRBR+gramRBR')*0.5;
 698           end
 699           [gramRBR,cholFlag]=chol(gramRBR);
 700           if  cholFlag == 0
 701               blockVectorR(:,activeMask) = blockVectorR(:,activeMask)/gramRBR;
 702           else
 703               warning('BLOPEX:lobpcg:ResidualNotFullRank',...
 704                   'The residual is not full rank.');
 705               break
 706           end
 707       else
 708           if isnumeric(operatorB)
 709               blockVectorBR(:,activeMask) = ...
 710                   operatorB*blockVectorR(:,activeMask);
 711           else
 712               blockVectorBR(:,activeMask) = ...
 713                   feval(operatorB,blockVectorR(:,activeMask));
 714           end
 715           gramRBR=blockVectorR(:,activeMask)'*blockVectorBR(:,activeMask);
 716           if ~isreal(gramRBR)
 717               gramRBR=(gramRBR+gramRBR')*0.5;
 718           end
 719           [gramRBR,cholFlag]=chol(gramRBR);
 720           if  cholFlag == 0
 721               blockVectorR(:,activeMask) = ...
 722                   blockVectorR(:,activeMask)/gramRBR;
 723               blockVectorBR(:,activeMask) = ...
 724                   blockVectorBR(:,activeMask)/gramRBR;
 725           else
 726               warning('BLOPEX:lobpcg:ResidualNotFullRankOrElse',...
 727               strcat('The residual is not full rank or/and operatorB ',...
 728               'is not positive definite.'));
 729               break
 730           end
 731
 732       end
 733       clear gramRBR;
 734
 735       if isnumeric(operatorA)
 736           blockVectorAR(:,activeMask) = operatorA*blockVectorR(:,activeMask);
 737       else
 738           blockVectorAR(:,activeMask) = ...
 739               feval(operatorA,blockVectorR(:,activeMask));
 740       end
 741
 742       if iterationNumber > 1
 743
 744           %Making active conjugate directions orthonormal
 745           if isempty(operatorB)
 746               %[blockVectorP(:,activeMask),gramPBP] = qr(blockVectorP(:,activeMask),0);
 747               gramPBP=blockVectorP(:,activeMask)'*blockVectorP(:,activeMask);
 748               if ~isreal(gramPBP)
 749                   gramPBP=(gramPBP+gramPBP')*0.5;
 750               end
 751               [gramPBP,cholFlag]=chol(gramPBP);
 752               if  cholFlag == 0
 753                   blockVectorP(:,activeMask) = ...
 754                       blockVectorP(:,activeMask)/gramPBP;
 755                   blockVectorAP(:,activeMask) = ...
 756                       blockVectorAP(:,activeMask)/gramPBP;
 757               else
 758                   warning('BLOPEX:lobpcg:DirectionNotFullRank',...
 759                       'The direction matrix is not full rank.');
 760                   break
 761               end
 762           else
 763               gramPBP=blockVectorP(:,activeMask)'*blockVectorBP(:,activeMask);
 764               if ~isreal(gramPBP)
 765                   gramPBP=(gramPBP+gramPBP')*0.5;
 766               end
 767               [gramPBP,cholFlag]=chol(gramPBP);
 768               if  cholFlag == 0
 769                   blockVectorP(:,activeMask) = ...
 770                       blockVectorP(:,activeMask)/gramPBP;
 771                   blockVectorAP(:,activeMask) = ...
 772                       blockVectorAP(:,activeMask)/gramPBP;
 773                   blockVectorBP(:,activeMask) = ...
 774                       blockVectorBP(:,activeMask)/gramPBP;
 775               else
 776                   warning('BLOPEX:lobpcg:DirectionNotFullRank',...
 777                  strcat('The direction matrix is not full rank ',...
 778               'or/and operatorB is not positive definite.'));
 779                   break
 780               end
 781           end
 782           clear gramPBP
 783       end
 784
 785       condestGmean = mean(condestGhistory(max(1,iterationNumber-10-...
 786           round(log(currentBlockSize))):iterationNumber));
 787
 788       %  restart=1;
 789
 790       % The Raileight-Ritz method for [blockVectorX blockVectorR blockVectorP]
 791
 792       if  residualNorms > eps^0.6
 793           explicitGramFlag = 0;
 794       else
 795           explicitGramFlag = 1;  %suggested by Garrett Moran, private
 796       end
 797
 798       activeRSize=size(blockVectorR(:,activeMask),2);
 799       if iterationNumber == 1
 800           activePSize=0;
 801           restart=1;
 802       else
 803           activePSize=size(blockVectorP(:,activeMask),2);
 804           restart=0;
 805       end
 806
 807       gramXAR=full(blockVectorAX'*blockVectorR(:,activeMask));
 808       gramRAR=full(blockVectorAR(:,activeMask)'*blockVectorR(:,activeMask));
 809       gramRAR=(gramRAR'+gramRAR)*0.5;
 810
 811       if explicitGramFlag
 812           gramXAX=full(blockVectorAX'*blockVectorX);
 813           gramXAX=(gramXAX'+gramXAX)*0.5;
 814           if isempty(operatorB)
 815               gramXBX=full(blockVectorX'*blockVectorX);
 816               gramRBR=full(blockVectorR(:,activeMask)'*...
 817                   blockVectorR(:,activeMask));
 818               gramXBR=full(blockVectorX'*blockVectorR(:,activeMask));
 819           else
 820               gramXBX=full(blockVectorBX'*blockVectorX);
 821               gramRBR=full(blockVectorBR(:,activeMask)'*...
 822                   blockVectorR(:,activeMask));
 823               gramXBR=full(blockVectorBX'*blockVectorR(:,activeMask));
 824           end
 825           gramXBX=(gramXBX'+gramXBX)*0.5;
 826           gramRBR=(gramRBR'+gramRBR)*0.5;
 827
 828       end
 829
 830       for cond_try=1:2,           %cond_try == 2 when restart
 831
 832           if ~restart
 833               gramXAP=full(blockVectorAX'*blockVectorP(:,activeMask));
 834               gramRAP=full(blockVectorAR(:,activeMask)'*...
 835                   blockVectorP(:,activeMask));
 836               gramPAP=full(blockVectorAP(:,activeMask)'*...
 837                   blockVectorP(:,activeMask));
 838               gramPAP=(gramPAP'+gramPAP)*0.5;
 839
 840               if explicitGramFlag
 841                   gramA = [ gramXAX     gramXAR     gramXAP
 842                       gramXAR'    gramRAR     gramRAP
 843                       gramXAP'     gramRAP'    gramPAP ];
 844               else
 845                   gramA = [ diag(lambda)  gramXAR  gramXAP
 846                       gramXAR'      gramRAR  gramRAP
 847                       gramXAP'      gramRAP'  gramPAP ];
 848               end
 849
 850               clear gramXAP  gramRAP gramPAP
 851
 852               if isempty(operatorB)
 853                   gramXBP=full(blockVectorX'*blockVectorP(:,activeMask));
 854                   gramRBP=full(blockVectorR(:,activeMask)'*...
 855                       blockVectorP(:,activeMask));
 856               else
 857                   gramXBP=full(blockVectorBX'*blockVectorP(:,activeMask));
 858                   gramRBP=full(blockVectorBR(:,activeMask)'*...
 859                       blockVectorP(:,activeMask));
 860                   %or blockVectorR(:,activeMask)'*blockVectorBP(:,activeMask);
 861               end
 862
 863               if explicitGramFlag
 864                   if isempty(operatorB)
 865                       gramPBP=full(blockVectorP(:,activeMask)'*...
 866                           blockVectorP(:,activeMask));
 867                   else
 868                       gramPBP=full(blockVectorBP(:,activeMask)'*...
 869                           blockVectorP(:,activeMask));
 870                   end
 871                   gramPBP=(gramPBP'+gramPBP)*0.5;
 872                   gramB = [ gramXBX  gramXBR  gramXBP
 873                       gramXBR' gramRBR  gramRBP
 874                       gramXBP' gramRBP' gramPBP ];
 875                   clear   gramPBP
 876               else
 877                   gramB=[eye(blockSize) zeros(blockSize,activeRSize) gramXBP
 878                       zeros(blockSize,activeRSize)' eye(activeRSize) gramRBP
 879                       gramXBP' gramRBP' eye(activePSize) ];
 880               end
 881
 882               clear gramXBP  gramRBP;
 883
 884           else
 885
 886               if explicitGramFlag
 887                   gramA = [ gramXAX   gramXAR
 888                       gramXAR'    gramRAR  ];
 889                   gramB = [ gramXBX  gramXBR
 890                       gramXBR' eye(activeRSize)  ];
 891                   clear gramXAX gramXBX gramXBR
 892               else
 893                   gramA = [ diag(lambda)  gramXAR
 894                       gramXAR'        gramRAR  ];
 895                   gramB = eye(blockSize+activeRSize);
 896               end
 897
 898               clear gramXAR gramRAR;
 899
 900           end
 901
 902           condestG = log10(cond(gramB))+1;
 903           if (condestG/condestGmean > 2 && condestG > 2 )|| condestG > 8
 904               %black magic - need to guess the restart
 905               if verbosityLevel
 906                   fprintf('Restart on step %i as condestG %5.4e \n',...
 907                       iterationNumber,condestG);
 908               end
 909               if cond_try == 1 && ~restart
 910                   restart=1; %steepest descent restart for stability
 911               else
 912                   warning('BLOPEX:lobpcg:IllConditioning',...
 913                       'Gramm matrix ill-conditioned: results unpredictable');
 914               end
 915           else
 916               break
 917           end
 918
 919       end
 920
 921       [gramA,gramB]=eig(gramA,gramB);
 922       lambda=diag(gramB(1:blockSize,1:blockSize));
 923       coordX=gramA(:,1:blockSize);
 924
 925       clear gramA gramB
 926
 927       if issparse(blockVectorX)
 928           coordX=sparse(coordX);
 929       end
 930
 931       if ~restart
 932           blockVectorP =  blockVectorR(:,activeMask)*...
 933               coordX(blockSize+1:blockSize+activeRSize,:) + ...
 934               blockVectorP(:,activeMask)*...
 935               coordX(blockSize+activeRSize+1:blockSize + ...
 936               activeRSize+activePSize,:);
 937           blockVectorAP = blockVectorAR(:,activeMask)*...
 938               coordX(blockSize+1:blockSize+activeRSize,:) + ...
 939               blockVectorAP(:,activeMask)*...
 940               coordX(blockSize+activeRSize+1:blockSize + ...
 941               activeRSize+activePSize,:);
 942           if ~isempty(operatorB)
 943               blockVectorBP = blockVectorBR(:,activeMask)*...
 944                   coordX(blockSize+1:blockSize+activeRSize,:) + ...
 945                   blockVectorBP(:,activeMask)*...
 946                   coordX(blockSize+activeRSize+1:blockSize+activeRSize+activePSize,:);
 947           end
 948       else %use block steepest descent
 949           blockVectorP =   blockVectorR(:,activeMask)*...
 950               coordX(blockSize+1:blockSize+activeRSize,:);
 951           blockVectorAP = blockVectorAR(:,activeMask)*...
 952               coordX(blockSize+1:blockSize+activeRSize,:);
 953           if ~isempty(operatorB)
 954               blockVectorBP = blockVectorBR(:,activeMask)*...
 955                   coordX(blockSize+1:blockSize+activeRSize,:);
 956           end
 957       end
 958
 959       blockVectorX = blockVectorX*coordX(1:blockSize,:) + blockVectorP;
 960       blockVectorAX=blockVectorAX*coordX(1:blockSize,:) + blockVectorAP;
 961       if ~isempty(operatorB)
 962           blockVectorBX=blockVectorBX*coordX(1:blockSize,:) + blockVectorBP;
 963       end
 964       clear coordX
 965       %%end RR
 966
 967       lambdaHistory(1:blockSize,iterationNumber+1)=lambda;
 968       condestGhistory(iterationNumber+1)=condestG;
 969
 970       if verbosityLevel
 971           fprintf('Iteration %i current block size %i \n',...
 972               iterationNumber,currentBlockSize);
 973           fprintf('Eigenvalues lambda %17.16e \n',lambda);
 974           fprintf('Residual Norms %e \n',residualNorms');
 975       end
 976   end
 977   %The main step of the method was the CG cycle: end
 978
 979   %Postprocessing
 980
 981   %Making sure blockVectorX's "exactly" satisfy the blockVectorY constrains??
 982
 983   %Making sure blockVectorX's are "exactly" othonormalized by final "exact" RR
 984   if isempty(operatorB)
 985       gramXBX=full(blockVectorX'*blockVectorX);
 986   else
 987       if isnumeric(operatorB)
 988           blockVectorBX = operatorB*blockVectorX;
 989       else
 990           blockVectorBX = feval(operatorB,blockVectorX);
 991       end
 992       gramXBX=full(blockVectorX'*blockVectorBX);
 993   end
 994   gramXBX=(gramXBX'+gramXBX)*0.5;
 995
 996   if isnumeric(operatorA)
 997       blockVectorAX = operatorA*blockVectorX;
 998   else
 999       blockVectorAX = feval(operatorA,blockVectorX);
1000   end
1001   gramXAX = full(blockVectorX'*blockVectorAX);
1002   gramXAX = (gramXAX + gramXAX')*0.5;
1003
1004   %Raileigh-Ritz for blockVectorX, which is already operatorB-orthonormal
1005   [coordX,gramXBX] = eig(gramXAX,gramXBX);
1006   lambda=diag(gramXBX);
1007
1008   if issparse(blockVectorX)
1009       coordX=sparse(coordX);
1010   end
1011
1012   blockVectorX  =   blockVectorX*coordX;
1013   blockVectorAX  =  blockVectorAX*coordX;
1014   if ~isempty(operatorB)
1015       blockVectorBX  =  blockVectorBX*coordX;
1016   end
1017
1018   %Computing all residuals
1019   if isempty(operatorB)
1020       if blockSize > 1
1021           blockVectorR = blockVectorAX - ...
1022               blockVectorX*spdiags(lambda,0,blockSize,blockSize);
1023       else
1024           blockVectorR = blockVectorAX - blockVectorX*lambda;
1025           %to make blockVectorR full when lambda is just a scalar
1026       end
1027   else
1028       if blockSize > 1
1029           blockVectorR=blockVectorAX - ...
1030               blockVectorBX*spdiags(lambda,0,blockSize,blockSize);
1031       else
1032           blockVectorR = blockVectorAX - blockVectorBX*lambda;
1033           %to make blockVectorR full when lambda is just a scalar
1034       end
1035   end
1036
1037   residualNorms=full(sqrt(sum(conj(blockVectorR).*blockVectorR)'));
1038   residualNormsHistory(1:blockSize,iterationNumber)=residualNorms;
1039
1040   if verbosityLevel
1041       fprintf('Final Eigenvalues lambda %17.16e \n',lambda);
1042       fprintf('Final Residual Norms %e \n',residualNorms');
1043   end
1044
1045   if verbosityLevel == 2
1046       whos
1047       figure(491)
1048       semilogy((abs(residualNormsHistory(1:blockSize,1:iterationNumber-1)))');
1049       title('Residuals for Different Eigenpairs','fontsize',16);
1050       ylabel('Eucledian norm of residuals','fontsize',16);
1051       xlabel('Iteration number','fontsize',16);
1052       %axis tight;
1053       %axis([0 maxIterations+1 1e-15 1e3])
1054       set(gca,'FontSize',14);
1055       figure(492);
1056       semilogy(abs((lambdaHistory(1:blockSize,1:iterationNumber)-...
1057           repmat(lambda,1,iterationNumber)))');
1058       title('Eigenvalue errors for Different Eigenpairs','fontsize',16);
1059       ylabel('Estimated eigenvalue errors','fontsize',16);
1060       xlabel('Iteration number','fontsize',16);
1061       %axis tight;
1062       %axis([0 maxIterations+1 1e-15 1e3])
1063       set(gca,'FontSize',14);
1064       drawnow;
1065   end
1066
1067   varargout(1)={failureFlag};
1068   varargout(2)={lambdaHistory(1:blockSize,1:iterationNumber)};
1069   varargout(3)={residualNormsHistory(1:blockSize,1:iterationNumber-1)};
1070 end