function [H, n, data] = numrd(example)
%numrd: model NEPv from numerical radius computation.
%
% 	[H, n, data] = numrd(example) generates NEPv from computation of 
%	the numerical radius of a square matrix B = A1+1i*A2 with Hermitian A1
% 	and A2 of size n.  The NEPv is in the form of 
%		H(x)x = lam*x,
%	where 
%		H(x) := (x'*A1*x)*A1 + (x'*A2*x)*A2
%	is a Hermitian matrix and lam is a largest eigenvalue H(x). 
% 	
%
%	INPUT: example 
%	'rnd' 	---	an NEPv of size n = 4, with 3 stable eigenvectors;
%   'grcar' ---	an NEPv of size n = 30, with 2 stable eigenvectors 
%				and 2 non-stable eigenvectors. 
%
%	OUTPUT: 
%	H	--- function handle H(x);
%	n	--- size n;
%	data--- data = {A1,A2}.
%
%
% -------
%	Usage:
%	ex = 'rnd';
%	[H, n] = rec(ex); 
%	Hx0 = H(randn(n,1)); % size 4*4
%
%	ex = 'grcar';
%	[H, n] = rec(ex); 
%	Hx0 = H(randn(n,1)); % size 30*20
%
%
% -------
% 	NEPv features: Hermitian H; Single vector X; Unitarily invariant; 
%	Largest eigenvalue; Multiple solutions; Non-stable eigenvector;
%
%	Ding.Lu@uky.edu. Version date July-20-2022.
%

if strcmp(example, 'rnd')
    n = 4;
	B = [ 
		 0.6 + 0.6i  -0.2 + 2.5i  -1.9 - 0.2i  -0.3 + 2.5i; 
		-0.1 + 2.3i  -0.3 - 2.6i  -1.3 + 0.4i  -1.2 + 1.3i;
	   	-2.0 + 0.0i  -1.6 + 0.6i  -2.1 - 0.4i   1.3 + 1.2i;
	   	-0.1 + 2.0i  -1.6 + 1.4i   1.5 + 1.0i  -0.1 - 2.3i];

elseif strcmp(example, 'grcar')
    n = 30;
    B = gallery('grcar', 30)* exp(1i*pi/5); % rotated grcar matrix

else
	disp("Error: Input example = 'sonar' or 'iono'.");
	return
end

A1 = (B + B')/2;
A2 = (B - B')/(2i);

% construct H(x) 
H = @(x) (x'*A1*x)*A1 + (x'*A2*x)*A2;
data = {A1, A2};

return;

end % END OF numrd


