X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Flinear-algebra-2.2.0%2Flobpcg.m;fp=octave_packages%2Flinear-algebra-2.2.0%2Flobpcg.m;h=30585911b5076a9bd31d7581685d4bb10a5e0455;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/linear-algebra-2.2.0/lobpcg.m b/octave_packages/linear-algebra-2.2.0/lobpcg.m new file mode 100644 index 0000000..3058591 --- /dev/null +++ b/octave_packages/linear-algebra-2.2.0/lobpcg.m @@ -0,0 +1,1070 @@ +## Copyright (C) 2000-2011 A.V. Knyazev +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU Lesser General Public License as published by the Free +## Software Foundation; either version 3 of the License, or (at your option) any +## later version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License +## for more details. +## +## You should have received a copy of the GNU Lesser General Public License +## along with this program; if not, see . + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{blockVectorX}, @var{lambda}] =} lobpcg (@var{blockVectorX}, @var{operatorA}) +## @deftypefnx {Function File} {[@var{blockVectorX}, @var{lambda}, @var{failureFlag}] =} lobpcg (@var{blockVectorX}, @var{operatorA}) +## @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}) +## Solves Hermitian partial eigenproblems using preconditioning. +## +## The first form outputs the array of algebraic smallest eigenvalues @var{lambda} and +## corresponding matrix of orthonormalized eigenvectors @var{blockVectorX} of the +## Hermitian (full or sparse) operator @var{operatorA} using input matrix +## @var{blockVectorX} as an initial guess, without preconditioning, somewhat +## similar to: +## +## @example +## # for real symmetric operator operatorA +## opts.issym = 1; opts.isreal = 1; K = size (blockVectorX, 2); +## [blockVectorX, lambda] = eigs (operatorA, K, 'SR', opts); +## +## # for Hermitian operator operatorA +## K = size (blockVectorX, 2); +## [blockVectorX, lambda] = eigs (operatorA, K, 'SR'); +## @end example +## +## The second form returns a convergence flag. If @var{failureFlag} is 0 then +## all the eigenvalues converged; otherwise not all converged. +## +## The third form computes smallest eigenvalues @var{lambda} and corresponding eigenvectors +## @var{blockVectorX} of the generalized eigenproblem Ax=lambda Bx, where +## Hermitian operators @var{operatorA} and @var{operatorB} are given as functions, as +## well as a preconditioner, @var{operatorT}. The operators @var{operatorB} and +## @var{operatorT} must be in addition @emph{positive definite}. To compute the largest +## eigenpairs of @var{operatorA}, simply apply the code to @var{operatorA} multiplied by +## -1. The code does not involve @emph{any} matrix factorizations of @var{operatorA} and +## @var{operatorB}, thus, e.g., it preserves the sparsity and the structure of +## @var{operatorA} and @var{operatorB}. +## +## @var{residualTolerance} and @var{maxIterations} control tolerance and max number of +## steps, and @var{verbosityLevel} = 0, 1, or 2 controls the amount of printed +## info. @var{lambdaHistory} is a matrix with all iterative lambdas, and +## @var{residualNormsHistory} are matrices of the history of 2-norms of residuals +## +## Required input: +## @itemize @bullet +## @item +## @var{blockVectorX} (class numeric) - initial approximation to eigenvectors, +## full or sparse matrix n-by-blockSize. @var{blockVectorX} must be full rank. +## @item +## @var{operatorA} (class numeric, char, or function_handle) - the main operator +## of the eigenproblem, can be a matrix, a function name, or handle +## @end itemize +## +## Optional function input: +## @itemize @bullet +## @item +## @var{operatorB} (class numeric, char, or function_handle) - the second operator, +## if solving a generalized eigenproblem, can be a matrix, a function name, or +## handle; by default if empty, @code{operatorB = I}. +## @item +## @var{operatorT} (class char or function_handle) - the preconditioner, by +## default @code{operatorT(blockVectorX) = blockVectorX}. +## @end itemize +## +## Optional constraints input: +## @itemize @bullet +## @item +## @var{blockVectorY} (class numeric) - a full or sparse n-by-sizeY matrix of +## constraints, where sizeY < n. @var{blockVectorY} must be full rank. The +## iterations will be performed in the (operatorB-) orthogonal complement of the +## column-space of @var{blockVectorY}. +## @end itemize +## +## Optional scalar input parameters: +## @itemize @bullet +## @item +## @var{residualTolerance} (class numeric) - tolerance, by default, @code{residualTolerance = n * sqrt (eps)} +## @item +## @var{maxIterations} - max number of iterations, by default, @code{maxIterations = min (n, 20)} +## @item +## @var{verbosityLevel} - either 0 (no info), 1, or 2 (with pictures); by +## default, @code{verbosityLevel = 0}. +## @end itemize +## +## Required output: +## @itemize @bullet +## @item +## @var{blockVectorX} and @var{lambda} (class numeric) both are computed +## blockSize eigenpairs, where @code{blockSize = size (blockVectorX, 2)} +## for the initial guess @var{blockVectorX} if it is full rank. +## @end itemize +## +## Optional output: +## @itemize @bullet +## @item +## @var{failureFlag} (class integer) as described above. +## @item +## @var{lambdaHistory} (class numeric) as described above. +## @item +## @var{residualNormsHistory} (class numeric) as described above. +## @end itemize +## +## Functions @code{operatorA(blockVectorX)}, @code{operatorB(blockVectorX)} and +## @code{operatorT(blockVectorX)} must support @var{blockVectorX} being a matrix, not +## just a column vector. +## +## Every iteration involves one application of @var{operatorA} and @var{operatorB}, and +## one of @var{operatorT}. +## +## Main memory requirements: 6 (9 if @code{isempty(operatorB)=0}) matrices of the +## same size as @var{blockVectorX}, 2 matrices of the same size as @var{blockVectorY} +## (if present), and two square matrices of the size 3*blockSize. +## +## In all examples below, we use the Laplacian operator in a 20x20 square +## with the mesh size 1 which can be generated in MATLAB by running: +## @example +## A = delsq (numgrid ('S', 21)); +## n = size (A, 1); +## @end example +## +## or in MATLAB and Octave by: +## @example +## [~,~,A] = laplacian ([19, 19]); +## n = size (A, 1); +## @end example +## +## Note that @code{laplacian} is a function of the specfun octave-forge package. +## +## The following Example: +## @example +## [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, 1e-5, 50, 2); +## @end example +## +## attempts to compute 8 first eigenpairs without preconditioning, but not all +## eigenpairs converge after 50 steps, so failureFlag=1. +## +## The next Example: +## @example +## blockVectorY = []; +## lambda_all = []; +## for j = 1:4 +## [blockVectorX, lambda] = lobpcg (randn (n, 2), A, blockVectorY, 1e-5, 200, 2); +## blockVectorY = [blockVectorY, blockVectorX]; +## lambda_all = [lambda_all' lambda']'; +## pause; +## end +## @end example +## +## attemps to compute the same 8 eigenpairs by calling the code 4 times with +## blockSize=2 using orthogonalization to the previously founded eigenvectors. +## +## The following Example: +## @example +## R = ichol (A, struct('michol', 'on')); +## precfun = @@(x)R\(R'\x); +## [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, [], @@(x)precfun(x), 1e-5, 60, 2); +## @end example +## +## computes the same eigenpairs in less then 25 steps, so that failureFlag=0 +## using the preconditioner function @code{precfun}, defined inline. If @code{precfun} +## is defined as an octave function in a file, the function handle +## @code{@@(x)precfun(x)} can be equivalently replaced by the function name @code{precfun}. Running: +## +## @example +## [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, speye (n), @@(x)precfun(x), 1e-5, 50, 2); +## @end example +## +## produces similar answers, but is somewhat slower and needs more memory as +## technically a generalized eigenproblem with B=I is solved here. +## +## The following example for a mostly diagonally dominant sparse matrix A +## demonstrates different types of preconditioning, compared to the standard +## use of the main diagonal of A: +## +## @example +## clear all; close all; +## n = 1000; +## M = spdiags ([1:n]', 0, n, n); +## precfun = @@(x)M\x; +## A = M + sprandsym (n, .1); +## Xini = randn (n, 5); +## maxiter = 15; +## tol = 1e-5; +## [~,~,~,~,rnp] = lobpcg (Xini, A, tol, maxiter, 1); +## [~,~,~,~,r] = lobpcg (Xini, A, [], @@(x)precfun(x), tol, maxiter, 1); +## subplot (2,2,1), semilogy (r'); hold on; +## semilogy (rnp', ':>'); +## title ('No preconditioning (top)'); axis tight; +## M(1,2) = 2; +## precfun = @@(x)M\x; % M is no longer symmetric +## [~,~,~,~,rns] = lobpcg (Xini, A, [], @@(x)precfun(x), tol, maxiter, 1); +## subplot (2,2,2), semilogy (r'); hold on; +## semilogy (rns', '--s'); +## title ('Nonsymmetric preconditioning (square)'); axis tight; +## M(1,2) = 0; +## precfun = @@(x)M\(x+10*sin(x)); % nonlinear preconditioning +## [~,~,~,~,rnl] = lobpcg (Xini, A, [], @@(x)precfun(x), tol, maxiter, 1); +## subplot (2,2,3), semilogy (r'); hold on; +## semilogy (rnl', '-.*'); +## title ('Nonlinear preconditioning (star)'); axis tight; +## M = abs (M - 3.5 * speye (n, n)); +## precfun = @@(x)M\x; +## [~,~,~,~,rs] = lobpcg (Xini, A, [], @@(x)precfun(x), tol, maxiter, 1); +## subplot (2,2,4), semilogy (r'); hold on; +## semilogy (rs', '-d'); +## title ('Selective preconditioning (diamond)'); axis tight; +## @end example +## +## @heading References +## This main function @code{lobpcg} is a version of the preconditioned conjugate +## gradient method (Algorithm 5.1) described in A. V. Knyazev, Toward the Optimal +## Preconditioned Eigensolver: +## Locally Optimal Block Preconditioned Conjugate Gradient Method, +## SIAM Journal on Scientific Computing 23 (2001), no. 2, pp. 517-541. +## @uref{http://dx.doi.org/10.1137/S1064827500366124} +## +## @heading Known bugs/features +## @itemize @bullet +## @item +## an excessively small requested tolerance may result in often restarts and +## instability. The code is not written to produce an eps-level accuracy! Use +## common sense. +## @item +## the code may be very sensitive to the number of eigenpairs computed, +## if there is a cluster of eigenvalues not completely included, cf. +## @example +## operatorA = diag ([1 1.99 2:99]); +## [blockVectorX, lambda] = lobpcg (randn (100, 1),operatorA, 1e-10, 80, 2); +## [blockVectorX, lambda] = lobpcg (randn (100, 2),operatorA, 1e-10, 80, 2); +## [blockVectorX, lambda] = lobpcg (randn (100, 3),operatorA, 1e-10, 80, 2); +## @end example +## @end itemize +## +## @heading Distribution +## The main distribution site: @uref{http://math.ucdenver.edu/~aknyazev/} +## +## A C-version of this code is a part of the @uref{http://code.google.com/p/blopex/} +## package and is directly available, e.g., in PETSc and HYPRE. +## @end deftypefn + +function [blockVectorX,lambda,varargout] = lobpcg(blockVectorX,operatorA,varargin) + %Begin + + % constants + + CONVENTIONAL_CONSTRAINTS = 1; + SYMMETRIC_CONSTRAINTS = 2; + + %Initial settings + + failureFlag = 1; + if nargin < 2 + error('BLOPEX:lobpcg:NotEnoughInputs',... + strcat('There must be at least 2 input agruments: ',... + 'blockVectorX and operatorA')); + end + if nargin > 8 + warning('BLOPEX:lobpcg:TooManyInputs',... + strcat('There must be at most 8 input agruments ',... + 'unless arguments are passed to a function')); + end + + if ~isnumeric(blockVectorX) + error('BLOPEX:lobpcg:FirstInputNotNumeric',... + 'The first input argument blockVectorX must be numeric'); + end + [n,blockSize]=size(blockVectorX); + if blockSize > n + error('BLOPEX:lobpcg:FirstInputFat',... + 'The first input argument blockVectorX must be tall, not fat'); + end + if n < 6 + error('BLOPEX:lobpcg:MatrixTooSmall',... + 'The code does not work for matrices of small sizes'); + end + + if isa(operatorA,'numeric') + nA = size(operatorA,1); + if any(size(operatorA) ~= nA) + error('BLOPEX:lobpcg:MatrixNotSquare',... + 'operatorA must be a square matrix or a string'); + end + if size(operatorA) ~= n + error('BLOPEX:lobpcg:MatrixWrongSize',... + ['The size ' int2str(size(operatorA))... + ' of operatorA is not the same as ' int2str(n)... + ' - the number of rows of blockVectorX']); + end + end + + count_string = 0; + + operatorT = []; + operatorB = []; + residualTolerance = []; + maxIterations = []; + verbosityLevel = []; + blockVectorY = []; sizeY = 0; + for j = 1:nargin-2 + if isequal(size(varargin{j}),[n,n]) + if isempty(operatorB) + operatorB = varargin{j}; + else + error('BLOPEX:lobpcg:TooManyMatrixInputs',... + strcat('Too many matrix input arguments. ',... + 'Preconditioner operatorT must be an M-function')); + end + elseif isequal(size(varargin{j},1),n) && size(varargin{j},2) < n + if isempty(blockVectorY) + blockVectorY = varargin{j}; + sizeY=size(blockVectorY,2); + else + error('BLOPEX:lobpcg:WrongConstraintsFormat',... + 'Something wrong with blockVectorY input argument'); + end + elseif ischar(varargin{j}) || isa(varargin{j},'function_handle') + if count_string == 0 + if isempty(operatorB) + operatorB = varargin{j}; + count_string = count_string + 1; + else + operatorT = varargin{j}; + end + elseif count_string == 1 + operatorT = varargin{j}; + else + warning('BLOPEX:lobpcg:TooManyStringFunctionHandleInputs',... + 'Too many string or FunctionHandle input arguments'); + end + elseif isequal(size(varargin{j}),[n,n]) + error('BLOPEX:lobpcg:WrongPreconditionerFormat',... + 'Preconditioner operatorT must be an M-function'); + elseif max(size(varargin{j})) == 1 + if isempty(residualTolerance) + residualTolerance = varargin{j}; + elseif isempty(maxIterations) + maxIterations = varargin{j}; + elseif isempty(verbosityLevel) + verbosityLevel = varargin{j}; + else + warning('BLOPEX:lobpcg:TooManyScalarInputs',... + 'Too many scalar parameters, need only three'); + end + elseif isempty(varargin{j}) + if isempty(operatorB) + count_string = count_string + 1; + elseif ~isempty(operatorT) + count_string = count_string + 1; + elseif ~isempty(blockVectorY) + error('BLOPEX:lobpcg:UnrecognizedEmptyInput',... + ['Unrecognized empty input argument number ' int2str(j+2)]); + end + else + error('BLOPEX:lobpcg:UnrecognizedInput',... + ['Input argument number ' int2str(j+2) ' not recognized.']); + end + end + + if verbosityLevel + if issparse(blockVectorX) + fprintf(['The sparse initial guess with %i colunms '... + 'and %i raws is detected \n'],n,blockSize); + else + fprintf(['The full initial guess with %i colunms '... + 'and %i raws is detected \n'],n,blockSize); + end + if ischar(operatorA) + fprintf('The main operator is detected as an M-function %s \n',... + operatorA); + elseif isa(operatorA,'function_handle') + fprintf('The main operator is detected as an M-function %s \n',... + func2str(operatorA)); + elseif issparse(operatorA) + fprintf('The main operator is detected as a sparse matrix \n'); + else + fprintf('The main operator is detected as a full matrix \n'); + end + if isempty(operatorB) + fprintf('Solving standard eigenvalue problem, not generalized \n'); + elseif ischar(operatorB) + fprintf(['The second operator of the generalized eigenproblem \n'... + 'is detected as an M-function %s \n'],operatorB); + elseif isa(operatorB,'function_handle') + fprintf(['The second operator of the generalized eigenproblem \n'... + 'is detected as an M-function %s \n'],func2str(operatorB)); + elseif issparse(operatorB) + fprintf(strcat('The second operator of the generalized',... + 'eigenproblem \n is detected as a sparse matrix \n')); + else + fprintf(strcat('The second operator of the generalized',... + 'eigenproblem \n is detected as a full matrix \n')); + end + if isempty(operatorT) + fprintf('No preconditioner is detected \n'); + elseif ischar(operatorT) + fprintf('The preconditioner is detected as an M-function %s \n',... + operatorT); + elseif isa(operatorT,'function_handle') + fprintf('The preconditioner is detected as an M-function %s \n',... + func2str(operatorT)); + end + if isempty(blockVectorY) + fprintf('No matrix of constraints is detected \n') + elseif issparse(blockVectorY) + fprintf('The sparse matrix of %i constraints is detected \n',sizeY); + else + fprintf('The full matrix of %i constraints is detected \n',sizeY); + end + if issparse(blockVectorY) ~= issparse(blockVectorX) + warning('BLOPEX:lobpcg:SparsityInconsistent',... + strcat('The sparsity formats of the initial guess and ',... + 'the constraints are inconsistent')); + end + end + + % Set defaults + + if isempty(residualTolerance) + residualTolerance = sqrt(eps)*n; + end + if isempty(maxIterations) + maxIterations = min(n,20); + end + if isempty(verbosityLevel) + verbosityLevel = 0; + end + + if verbosityLevel + fprintf('Tolerance %e and maximum number of iterations %i \n',... + residualTolerance,maxIterations) + end + + %constraints preprocessing + if isempty(blockVectorY) + constraintStyle = 0; + else + % constraintStyle = SYMMETRIC_CONSTRAINTS; % more accurate? + constraintStyle = CONVENTIONAL_CONSTRAINTS; + end + + if constraintStyle == CONVENTIONAL_CONSTRAINTS + + if isempty(operatorB) + gramY = blockVectorY'*blockVectorY; + else + if isnumeric(operatorB) + blockVectorBY = operatorB*blockVectorY; + else + blockVectorBY = feval(operatorB,blockVectorY); + end + gramY=blockVectorY'*blockVectorBY; + end + gramY=(gramY'+gramY)*0.5; + if isempty(operatorB) + blockVectorX = blockVectorX - ... + blockVectorY*(gramY\(blockVectorY'*blockVectorX)); + else + blockVectorX =blockVectorX - ... + blockVectorY*(gramY\(blockVectorBY'*blockVectorX)); + end + + elseif constraintStyle == SYMMETRIC_CONSTRAINTS + + if ~isempty(operatorB) + if isnumeric(operatorB) + blockVectorY = operatorB*blockVectorY; + else + blockVectorY = feval(operatorB,blockVectorY); + end + end + if isempty(operatorT) + gramY = blockVectorY'*blockVectorY; + else + blockVectorTY = feval(operatorT,blockVectorY); + gramY = blockVectorY'*blockVectorTY; + end + gramY=(gramY'+gramY)*0.5; + if isempty(operatorT) + blockVectorX = blockVectorX - ... + blockVectorY*(gramY\(blockVectorY'*blockVectorX)); + else + blockVectorX = blockVectorX - ... + blockVectorTY*(gramY\(blockVectorY'*blockVectorX)); + end + + end + + %Making the initial vectors (operatorB-) orthonormal + if isempty(operatorB) + %[blockVectorX,gramXBX] = qr(blockVectorX,0); + gramXBX=blockVectorX'*blockVectorX; + if ~isreal(gramXBX) + gramXBX=(gramXBX+gramXBX')*0.5; + end + [gramXBX,cholFlag]=chol(gramXBX); + if cholFlag ~= 0 + error('BLOPEX:lobpcg:ConstraintsTooTight',... + 'The initial approximation after constraints is not full rank'); + end + blockVectorX = blockVectorX/gramXBX; + else + %[blockVectorX,blockVectorBX] = orth(operatorB,blockVectorX); + if isnumeric(operatorB) + blockVectorBX = operatorB*blockVectorX; + else + blockVectorBX = feval(operatorB,blockVectorX); + end + gramXBX=blockVectorX'*blockVectorBX; + if ~isreal(gramXBX) + gramXBX=(gramXBX+gramXBX')*0.5; + end + [gramXBX,cholFlag]=chol(gramXBX); + if cholFlag ~= 0 + error('BLOPEX:lobpcg:InitialNotFullRank',... + sprintf('%s\n%s', ... + 'The initial approximation after constraints is not full rank',... + 'or/and operatorB is not positive definite')); + end + blockVectorX = blockVectorX/gramXBX; + blockVectorBX = blockVectorBX/gramXBX; + end + + % Checking if the problem is big enough for the algorithm, + % i.e. n-sizeY > 5*blockSize + % Theoretically, the algorithm should be able to run if + % n-sizeY > 3*blockSize, + % but the extreme cases might be unstable, so we use 5 instead of 3 here. + if n-sizeY < 5*blockSize + error('BLOPEX:lobpcg:MatrixTooSmall','%s\n%s', ... + 'The problem size is too small, relative to the block size.',... + 'Try using eig() or eigs() instead.'); + end + + % Preallocation + residualNormsHistory=zeros(blockSize,maxIterations); + lambdaHistory=zeros(blockSize,maxIterations+1); + condestGhistory=zeros(1,maxIterations+1); + + blockVectorBR=zeros(n,blockSize); + blockVectorAR=zeros(n,blockSize); + blockVectorP=zeros(n,blockSize); + blockVectorAP=zeros(n,blockSize); + blockVectorBP=zeros(n,blockSize); + + %Initial settings for the loop + if isnumeric(operatorA) + blockVectorAX = operatorA*blockVectorX; + else + blockVectorAX = feval(operatorA,blockVectorX); + end + + gramXAX = full(blockVectorX'*blockVectorAX); + gramXAX = (gramXAX + gramXAX')*0.5; + % eig(...,'chol') uses only the diagonal and upper triangle - + % not true in MATLAB + % Octave v3.2.3-4, eig() does not support inputting 'chol' + [coordX,gramXAX]=eig(gramXAX,eye(blockSize)); + + lambda=diag(gramXAX); %eig returns non-ordered eigenvalues on the diagonal + + if issparse(blockVectorX) + coordX=sparse(coordX); + end + + blockVectorX = blockVectorX*coordX; + blockVectorAX = blockVectorAX*coordX; + if ~isempty(operatorB) + blockVectorBX = blockVectorBX*coordX; + end + clear coordX + + condestGhistory(1)=-log10(eps)/2; %if too small cause unnecessary restarts + + lambdaHistory(1:blockSize,1)=lambda; + + activeMask = true(blockSize,1); + % currentBlockSize = blockSize; %iterate all + % + % restart=1;%steepest descent + + %The main part of the method is the loop of the CG method: begin + for iterationNumber=1:maxIterations + + % %Computing the active residuals + % if isempty(operatorB) + % if currentBlockSize > 1 + % blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ... + % blockVectorX(:,activeMask)*spdiags(lambda(activeMask),0,currentBlockSize,currentBlockSize); + % else + % blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ... + % blockVectorX(:,activeMask)*lambda(activeMask); + % %to make blockVectorR full when lambda is just a scalar + % end + % else + % if currentBlockSize > 1 + % blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ... + % blockVectorBX(:,activeMask)*spdiags(lambda(activeMask),0,currentBlockSize,currentBlockSize); + % else + % blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ... + % blockVectorBX(:,activeMask)*lambda(activeMask); + % %to make blockVectorR full when lambda is just a scalar + % end + % end + + %Computing all residuals + if isempty(operatorB) + if blockSize > 1 + blockVectorR = blockVectorAX - ... + blockVectorX*spdiags(lambda,0,blockSize,blockSize); + else + blockVectorR = blockVectorAX - blockVectorX*lambda; + %to make blockVectorR full when lambda is just a scalar + end + else + if blockSize > 1 + blockVectorR=blockVectorAX - ... + blockVectorBX*spdiags(lambda,0,blockSize,blockSize); + else + blockVectorR = blockVectorAX - blockVectorBX*lambda; + %to make blockVectorR full when lambda is just a scalar + end + end + + %Satisfying the constraints for the active residulas + if constraintStyle == SYMMETRIC_CONSTRAINTS + if isempty(operatorT) + blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ... + blockVectorY*(gramY\(blockVectorY'*... + blockVectorR(:,activeMask))); + else + blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ... + blockVectorY*(gramY\(blockVectorTY'*... + blockVectorR(:,activeMask))); + end + end + + residualNorms=full(sqrt(sum(conj(blockVectorR).*blockVectorR)')); + residualNormsHistory(1:blockSize,iterationNumber)=residualNorms; + + %index antifreeze + activeMask = full(residualNorms > residualTolerance) & activeMask; + %activeMask = full(residualNorms > residualTolerance); + %above allows vectors back into active, which causes problems with frosen Ps + %activeMask = full(residualNorms > 0); %iterate all, ignore freeze + + currentBlockSize = sum(activeMask); + if currentBlockSize == 0 + failureFlag=0; %all eigenpairs converged + break + end + + %Applying the preconditioner operatorT to the active residulas + if ~isempty(operatorT) + blockVectorR(:,activeMask) = ... + feval(operatorT,blockVectorR(:,activeMask)); + end + + if constraintStyle == CONVENTIONAL_CONSTRAINTS + if isempty(operatorB) + blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ... + blockVectorY*(gramY\(blockVectorY'*... + blockVectorR(:,activeMask))); + else + blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ... + blockVectorY*(gramY\(blockVectorBY'*... + blockVectorR(:,activeMask))); + end + end + + %Making active (preconditioned) residuals orthogonal to blockVectorX + if isempty(operatorB) + blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ... + blockVectorX*(blockVectorX'*blockVectorR(:,activeMask)); + else + blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ... + blockVectorX*(blockVectorBX'*blockVectorR(:,activeMask)); + end + + %Making active residuals orthonormal + if isempty(operatorB) + %[blockVectorR(:,activeMask),gramRBR]=... + %qr(blockVectorR(:,activeMask),0); %to increase stability + gramRBR=blockVectorR(:,activeMask)'*blockVectorR(:,activeMask); + if ~isreal(gramRBR) + gramRBR=(gramRBR+gramRBR')*0.5; + end + [gramRBR,cholFlag]=chol(gramRBR); + if cholFlag == 0 + blockVectorR(:,activeMask) = blockVectorR(:,activeMask)/gramRBR; + else + warning('BLOPEX:lobpcg:ResidualNotFullRank',... + 'The residual is not full rank.'); + break + end + else + if isnumeric(operatorB) + blockVectorBR(:,activeMask) = ... + operatorB*blockVectorR(:,activeMask); + else + blockVectorBR(:,activeMask) = ... + feval(operatorB,blockVectorR(:,activeMask)); + end + gramRBR=blockVectorR(:,activeMask)'*blockVectorBR(:,activeMask); + if ~isreal(gramRBR) + gramRBR=(gramRBR+gramRBR')*0.5; + end + [gramRBR,cholFlag]=chol(gramRBR); + if cholFlag == 0 + blockVectorR(:,activeMask) = ... + blockVectorR(:,activeMask)/gramRBR; + blockVectorBR(:,activeMask) = ... + blockVectorBR(:,activeMask)/gramRBR; + else + warning('BLOPEX:lobpcg:ResidualNotFullRankOrElse',... + strcat('The residual is not full rank or/and operatorB ',... + 'is not positive definite.')); + break + end + + end + clear gramRBR; + + if isnumeric(operatorA) + blockVectorAR(:,activeMask) = operatorA*blockVectorR(:,activeMask); + else + blockVectorAR(:,activeMask) = ... + feval(operatorA,blockVectorR(:,activeMask)); + end + + if iterationNumber > 1 + + %Making active conjugate directions orthonormal + if isempty(operatorB) + %[blockVectorP(:,activeMask),gramPBP] = qr(blockVectorP(:,activeMask),0); + gramPBP=blockVectorP(:,activeMask)'*blockVectorP(:,activeMask); + if ~isreal(gramPBP) + gramPBP=(gramPBP+gramPBP')*0.5; + end + [gramPBP,cholFlag]=chol(gramPBP); + if cholFlag == 0 + blockVectorP(:,activeMask) = ... + blockVectorP(:,activeMask)/gramPBP; + blockVectorAP(:,activeMask) = ... + blockVectorAP(:,activeMask)/gramPBP; + else + warning('BLOPEX:lobpcg:DirectionNotFullRank',... + 'The direction matrix is not full rank.'); + break + end + else + gramPBP=blockVectorP(:,activeMask)'*blockVectorBP(:,activeMask); + if ~isreal(gramPBP) + gramPBP=(gramPBP+gramPBP')*0.5; + end + [gramPBP,cholFlag]=chol(gramPBP); + if cholFlag == 0 + blockVectorP(:,activeMask) = ... + blockVectorP(:,activeMask)/gramPBP; + blockVectorAP(:,activeMask) = ... + blockVectorAP(:,activeMask)/gramPBP; + blockVectorBP(:,activeMask) = ... + blockVectorBP(:,activeMask)/gramPBP; + else + warning('BLOPEX:lobpcg:DirectionNotFullRank',... + strcat('The direction matrix is not full rank ',... + 'or/and operatorB is not positive definite.')); + break + end + end + clear gramPBP + end + + condestGmean = mean(condestGhistory(max(1,iterationNumber-10-... + round(log(currentBlockSize))):iterationNumber)); + + % restart=1; + + % The Raileight-Ritz method for [blockVectorX blockVectorR blockVectorP] + + if residualNorms > eps^0.6 + explicitGramFlag = 0; + else + explicitGramFlag = 1; %suggested by Garrett Moran, private + end + + activeRSize=size(blockVectorR(:,activeMask),2); + if iterationNumber == 1 + activePSize=0; + restart=1; + else + activePSize=size(blockVectorP(:,activeMask),2); + restart=0; + end + + gramXAR=full(blockVectorAX'*blockVectorR(:,activeMask)); + gramRAR=full(blockVectorAR(:,activeMask)'*blockVectorR(:,activeMask)); + gramRAR=(gramRAR'+gramRAR)*0.5; + + if explicitGramFlag + gramXAX=full(blockVectorAX'*blockVectorX); + gramXAX=(gramXAX'+gramXAX)*0.5; + if isempty(operatorB) + gramXBX=full(blockVectorX'*blockVectorX); + gramRBR=full(blockVectorR(:,activeMask)'*... + blockVectorR(:,activeMask)); + gramXBR=full(blockVectorX'*blockVectorR(:,activeMask)); + else + gramXBX=full(blockVectorBX'*blockVectorX); + gramRBR=full(blockVectorBR(:,activeMask)'*... + blockVectorR(:,activeMask)); + gramXBR=full(blockVectorBX'*blockVectorR(:,activeMask)); + end + gramXBX=(gramXBX'+gramXBX)*0.5; + gramRBR=(gramRBR'+gramRBR)*0.5; + + end + + for cond_try=1:2, %cond_try == 2 when restart + + if ~restart + gramXAP=full(blockVectorAX'*blockVectorP(:,activeMask)); + gramRAP=full(blockVectorAR(:,activeMask)'*... + blockVectorP(:,activeMask)); + gramPAP=full(blockVectorAP(:,activeMask)'*... + blockVectorP(:,activeMask)); + gramPAP=(gramPAP'+gramPAP)*0.5; + + if explicitGramFlag + gramA = [ gramXAX gramXAR gramXAP + gramXAR' gramRAR gramRAP + gramXAP' gramRAP' gramPAP ]; + else + gramA = [ diag(lambda) gramXAR gramXAP + gramXAR' gramRAR gramRAP + gramXAP' gramRAP' gramPAP ]; + end + + clear gramXAP gramRAP gramPAP + + if isempty(operatorB) + gramXBP=full(blockVectorX'*blockVectorP(:,activeMask)); + gramRBP=full(blockVectorR(:,activeMask)'*... + blockVectorP(:,activeMask)); + else + gramXBP=full(blockVectorBX'*blockVectorP(:,activeMask)); + gramRBP=full(blockVectorBR(:,activeMask)'*... + blockVectorP(:,activeMask)); + %or blockVectorR(:,activeMask)'*blockVectorBP(:,activeMask); + end + + if explicitGramFlag + if isempty(operatorB) + gramPBP=full(blockVectorP(:,activeMask)'*... + blockVectorP(:,activeMask)); + else + gramPBP=full(blockVectorBP(:,activeMask)'*... + blockVectorP(:,activeMask)); + end + gramPBP=(gramPBP'+gramPBP)*0.5; + gramB = [ gramXBX gramXBR gramXBP + gramXBR' gramRBR gramRBP + gramXBP' gramRBP' gramPBP ]; + clear gramPBP + else + gramB=[eye(blockSize) zeros(blockSize,activeRSize) gramXBP + zeros(blockSize,activeRSize)' eye(activeRSize) gramRBP + gramXBP' gramRBP' eye(activePSize) ]; + end + + clear gramXBP gramRBP; + + else + + if explicitGramFlag + gramA = [ gramXAX gramXAR + gramXAR' gramRAR ]; + gramB = [ gramXBX gramXBR + gramXBR' eye(activeRSize) ]; + clear gramXAX gramXBX gramXBR + else + gramA = [ diag(lambda) gramXAR + gramXAR' gramRAR ]; + gramB = eye(blockSize+activeRSize); + end + + clear gramXAR gramRAR; + + end + + condestG = log10(cond(gramB))+1; + if (condestG/condestGmean > 2 && condestG > 2 )|| condestG > 8 + %black magic - need to guess the restart + if verbosityLevel + fprintf('Restart on step %i as condestG %5.4e \n',... + iterationNumber,condestG); + end + if cond_try == 1 && ~restart + restart=1; %steepest descent restart for stability + else + warning('BLOPEX:lobpcg:IllConditioning',... + 'Gramm matrix ill-conditioned: results unpredictable'); + end + else + break + end + + end + + [gramA,gramB]=eig(gramA,gramB); + lambda=diag(gramB(1:blockSize,1:blockSize)); + coordX=gramA(:,1:blockSize); + + clear gramA gramB + + if issparse(blockVectorX) + coordX=sparse(coordX); + end + + if ~restart + blockVectorP = blockVectorR(:,activeMask)*... + coordX(blockSize+1:blockSize+activeRSize,:) + ... + blockVectorP(:,activeMask)*... + coordX(blockSize+activeRSize+1:blockSize + ... + activeRSize+activePSize,:); + blockVectorAP = blockVectorAR(:,activeMask)*... + coordX(blockSize+1:blockSize+activeRSize,:) + ... + blockVectorAP(:,activeMask)*... + coordX(blockSize+activeRSize+1:blockSize + ... + activeRSize+activePSize,:); + if ~isempty(operatorB) + blockVectorBP = blockVectorBR(:,activeMask)*... + coordX(blockSize+1:blockSize+activeRSize,:) + ... + blockVectorBP(:,activeMask)*... + coordX(blockSize+activeRSize+1:blockSize+activeRSize+activePSize,:); + end + else %use block steepest descent + blockVectorP = blockVectorR(:,activeMask)*... + coordX(blockSize+1:blockSize+activeRSize,:); + blockVectorAP = blockVectorAR(:,activeMask)*... + coordX(blockSize+1:blockSize+activeRSize,:); + if ~isempty(operatorB) + blockVectorBP = blockVectorBR(:,activeMask)*... + coordX(blockSize+1:blockSize+activeRSize,:); + end + end + + blockVectorX = blockVectorX*coordX(1:blockSize,:) + blockVectorP; + blockVectorAX=blockVectorAX*coordX(1:blockSize,:) + blockVectorAP; + if ~isempty(operatorB) + blockVectorBX=blockVectorBX*coordX(1:blockSize,:) + blockVectorBP; + end + clear coordX + %%end RR + + lambdaHistory(1:blockSize,iterationNumber+1)=lambda; + condestGhistory(iterationNumber+1)=condestG; + + if verbosityLevel + fprintf('Iteration %i current block size %i \n',... + iterationNumber,currentBlockSize); + fprintf('Eigenvalues lambda %17.16e \n',lambda); + fprintf('Residual Norms %e \n',residualNorms'); + end + end + %The main step of the method was the CG cycle: end + + %Postprocessing + + %Making sure blockVectorX's "exactly" satisfy the blockVectorY constrains?? + + %Making sure blockVectorX's are "exactly" othonormalized by final "exact" RR + if isempty(operatorB) + gramXBX=full(blockVectorX'*blockVectorX); + else + if isnumeric(operatorB) + blockVectorBX = operatorB*blockVectorX; + else + blockVectorBX = feval(operatorB,blockVectorX); + end + gramXBX=full(blockVectorX'*blockVectorBX); + end + gramXBX=(gramXBX'+gramXBX)*0.5; + + if isnumeric(operatorA) + blockVectorAX = operatorA*blockVectorX; + else + blockVectorAX = feval(operatorA,blockVectorX); + end + gramXAX = full(blockVectorX'*blockVectorAX); + gramXAX = (gramXAX + gramXAX')*0.5; + + %Raileigh-Ritz for blockVectorX, which is already operatorB-orthonormal + [coordX,gramXBX] = eig(gramXAX,gramXBX); + lambda=diag(gramXBX); + + if issparse(blockVectorX) + coordX=sparse(coordX); + end + + blockVectorX = blockVectorX*coordX; + blockVectorAX = blockVectorAX*coordX; + if ~isempty(operatorB) + blockVectorBX = blockVectorBX*coordX; + end + + %Computing all residuals + if isempty(operatorB) + if blockSize > 1 + blockVectorR = blockVectorAX - ... + blockVectorX*spdiags(lambda,0,blockSize,blockSize); + else + blockVectorR = blockVectorAX - blockVectorX*lambda; + %to make blockVectorR full when lambda is just a scalar + end + else + if blockSize > 1 + blockVectorR=blockVectorAX - ... + blockVectorBX*spdiags(lambda,0,blockSize,blockSize); + else + blockVectorR = blockVectorAX - blockVectorBX*lambda; + %to make blockVectorR full when lambda is just a scalar + end + end + + residualNorms=full(sqrt(sum(conj(blockVectorR).*blockVectorR)')); + residualNormsHistory(1:blockSize,iterationNumber)=residualNorms; + + if verbosityLevel + fprintf('Final Eigenvalues lambda %17.16e \n',lambda); + fprintf('Final Residual Norms %e \n',residualNorms'); + end + + if verbosityLevel == 2 + whos + figure(491) + semilogy((abs(residualNormsHistory(1:blockSize,1:iterationNumber-1)))'); + title('Residuals for Different Eigenpairs','fontsize',16); + ylabel('Eucledian norm of residuals','fontsize',16); + xlabel('Iteration number','fontsize',16); + %axis tight; + %axis([0 maxIterations+1 1e-15 1e3]) + set(gca,'FontSize',14); + figure(492); + semilogy(abs((lambdaHistory(1:blockSize,1:iterationNumber)-... + repmat(lambda,1,iterationNumber)))'); + title('Eigenvalue errors for Different Eigenpairs','fontsize',16); + ylabel('Estimated eigenvalue errors','fontsize',16); + xlabel('Iteration number','fontsize',16); + %axis tight; + %axis([0 maxIterations+1 1e-15 1e3]) + set(gca,'FontSize',14); + drawnow; + end + + varargout(1)={failureFlag}; + varargout(2)={lambdaHistory(1:blockSize,1:iterationNumber)}; + varargout(3)={residualNormsHistory(1:blockSize,1:iterationNumber-1)}; +end