rosser_recipe

function R = rosser_recipe(classname)
%ROSSER_RECIPE constructs the classic symmetric eigenvalue test.
%   This matrix was a challenge for many matrix eigenvalue
%   algorithms. But LAPACK's DSYEV routine used in MATLAB has
%   no trouble with it. The matrix is 8-by-8 with integer
%   elements.
%
%   It has:
%       * A double eigenvalue.
%       * Three nearly equal eigenvalues.
%       * Dominant eigenvalues of opposite sign.
%       * A zero eigenvalue.
%       * A small, nonzero eigenvalue.
%
%   Modify the matrices M[1...4] to use your favorite
%   eigenvalues.
%
%   ROSSER_RECIPE(CLASSNAME) produces a matrix of class
%   CLASSNAME. CLASSNAME must be either 'single' or 'double'
%   (the default).
%
%   Created by David A. Gross. Inspired by [1].
%   Construction from [2,3].
%
%REFERENCES
%   [1] http://blogs.mathworks.com/cleve/2014/01/06/ ...
%   the-rosser-matrix/, accessed on 2014/01/07
%
%   [2] Rosser, J.B.; Lanczos, C.; Hestenes, M.R.; Karush, W.
%   Separation of close eigenvalues of a real symmetric matrix
%   (1951), J. Res. Natl. Bur. Stand., Vol. 47, No. 4, p. 291,
%   Appendix 1, https://archive.org/details/jresv47n4p291,
%   accessed on 2014/01/07
%
%   [3] T. Muir, History of Determinants III, 289 (Macmillan
%   and Co., Ltd., London, 1920), http://igm.univ-mlv.fr/ ...
%   ~al/Classiques/Muir/History_3/, accessed on 2014/01/07

%   Description Copyright 1984-2005 The MathWorks, Inc.
%   $Revision: 5.10.4.2 $  $Date: 2005/11/18 14:15:39 $

if nargin < 1, classname = 'double'; end

% make our eigenvalues in 2x2 matrices
M1 = [102  1 ;  1 -102]; % lambda = ± sqrt(102^2 + 1)
M2 = [101  1 ;  1  101]; % lambda = 101 ± 1
M3 = [  1 10 ; 10  101]; % lambda = 51 ± sqrt(51^2-1)
M4 = [ 98 14 ; 14    2]; % lambda = 100, 0

B = zeros(8);

% explode M[1...4] into an 8x8 matrix
B([1,6],[1,6]) = M1;
B([2,8],[2,8]) = M2;
B([4,5],[4,5]) = M3;
B([3,7],[3,7]) = M4;

sylvester88_A = @(a,b,c,d) [ ...
    a  b  c  d ; ...
    b -a -d  c ; ...
    c  d -a -b ; ...
    d -c  b -a ];

sylvester44 = @(a,b,c,d) [ ...
    a  b  c  d ; ...
    b -a  d -c ; ...
    c -d -a  b ; ...
    d  c -b -a ];

% make Sylvester's "penorthogonant" of determinant 10^8
P = blkdiag(sylvester88_A(2,1,1,2),sylvester44(1,-1,-2,2));

% P'*P = 10I
R = P'*B*P;

R = cast(R,classname);

end