function [H, n, data] = dhdae(example)
%DHDAE: model NEPv from distance problem of dHDAE system
%
% 	[H, n, data] = dhdae(example) generates NEPv from computation of 
%	distance to singularity of dHDEAE systems. The NEPv is in the form of 
%		H(x)x = lam*x,
%	where 
%		H(x) := A1 + (x'*A2*x)*A2 + ... + (x'*Am*x)*Am 
%	is a symmetric matrix and lam is a largest eigenvalue H(x). 
% 	
%
%	INPUT: example 
%	'linear'	--- a randomly generated linear dHDAE system of size n=30;
%   	'quadratic'	---	a quadratic, lumped-parameter mass-spring damper
%   				system of size n = 30, with a randomly generated
%   				gyroscopic matrix.
%
%	OUTPUT: 
%	H	--- function handle H(x);
%	n	--- size n;
%	data--- data = {A1,A2, ..., Am}.
%
%
% -------
%	Usage:
%	ex = 'linear';
%	[H, n] = rec(ex); 
%	Hx0 = H(randn(n,1)); % size 30 symmetric matrix
%
%	ex = 'quadratic';
%	[H, n] = rec(ex); 
%	Hx0 = H(randn(n,1)); % size 30 symmetric matrix
%
%
%
% -------
% 	NEPv features: Symmetric H; Single vector X; Unitarily invariant; 
%	Largest eigenvalue; 
%
%	Ding.Lu@uky.edu. Version date Sept-20-2022.
%

rng(0); % fix random seed for reproducible testing

if strcmp(example, 'linear') 
    m = 3; 
    n = 30;
    B1 = orth(rand(n,n)); B1 = B1 *diag(rand(n,1)+1.0E-6)* B1';
    B2 = orth(rand(n,n)); B2 = B2 *diag(rand(n,1)+1.0E-6)* B2';
    J = rand(n,n); J = (J-J'); J =J/norm(J); 
    A1 = (J*J - B1*B1 - B2*B2);

    A{1} = A1; A{2} = B1; A{3} = B2;

elseif strcmp(example, 'quadratic') 
    m = 4; n = 30;

    % dynamical system
    mm = [1:n]; dd = mm; kk =mm; 
    M=zeros(n); K=zeros(n); D=zeros(n);
    for j=1:n-1
        M(j,j)=mm(j);
        K(j,j)=kk(j)+kk(j+1);
        K(j,j+1)=-kk(j+1);
        K(j+1,j)=K(j,j+1);
        D(j,j)=dd(j)+dd(j+1);
        D(j,j+1)=-dd(j+1);
        D(j+1,j)=D(j,j+1);
    end
    M(n,n)=mm(n); K(n,n)=kk(n); D(n,n)=dd(n);
    B1 = M/n; B2=D/n; B3=K/n; 
    G = randn(n); G = (G-G')/2; G = G/norm(G);

    A1 = (G*G - B1*B1 - B2*B2 -B3*B3);
    A{1} = A1; A{2} = B1; A{3} = B2; A{4} = B3;

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

% construct H(x) 
H = @(x) Hsub(x, A);
data = A;

return;

end % END OF dhdae


% *******************************************
% Local functions
% *******************************************
function H = Hsub(x, A)
H = A{1};
for i = 2:length(A)
    H = H+(x'*A{i}*x)*A{i};
end
end % END OF Hsub
