% quad_eig.m
%
% quadratic eigenproblem script
% solves the problem
%
% ( lambda^2 * M + lambda * C + K ) * x = 0
%
% where M,C,K are sparse matrices
%

% MODIFY these three lines

M = sparse( [ 2 1; 1 3] );
C = sparse( [ 4 0; 0 5] );
K = sparse( [ 6 0; 0 7] );

% MODIFY the shift value, do not use zero for singular matrices A

alpha = -1;

% MODIFY the number of eigenvalues you want to find

neig = 4;



rows = size(M,1);
cols = size(M,2);


% linearization of problem          Ax = lambda Bx
% ------------------------
% etemplates p. 283 9.4, 9.5, 9.6

A = sparse( [ zeros( rows, cols ), eye( rows, cols ); -K, -C ] );

B = sparse( [ eye( rows, cols ), zeros( rows, cols ); zeros( rows, cols ), M]);


% solve linear problem
% --------------------

% call ppp, pass in problem type 2, linear generalized Hermitian

problem = 2;
matbal = 0;

fulldual = 0;          % find eigenvalues and eigenvectors
semidual = 1;
group    = 1;
treatbd  = 1;

nb = 1;
maxit = size(A,1) * max( rows, cols );

precond = 1;    % LU factorization must be done, but L and U don't have to
L = 0;          % be provided, will be computed by ppp
U = 0;

[ neig, ritz, eigvecl, eigvecr, resl, resr, condnum, ... 
  ritzvalue, szoft, duality, exdual,omega,resl_e,resr_e,tolconv,numlsteps] ...
= ppp( problem, matbal, A, B, neig, ...
       fulldual, semidual, group, treatbd, ...
       nb, maxit, precond, alpha, L, U );


% convert eigenvectors of linear problem to eigenvectors of quadratic problem
% --------------------
%
% eigvecl, eigvecr : eigenvectors of linear problem, in columns
% evl, evr         : eigenvectors of quadratic problem, in columns

evl = eigvecl( rows + 1 : 2 * rows, : );

evr = eigvecr( 1 : rows, : );


% display results
% ---------------
% eigenvalues are in ritz array, first neig values are valid

fprintf( '\nfirst %i eigenvalues close to %f + %fi are\n', ...
         neig, real(alpha), imag(alpha) );
display( ritz( 1:neig ) );

% find quadratic problem residuals
% --------------------------------

% largest matrix norm

matrix_norm = max( [ norm(M,1) norm(C,1) norm(K,1) ] );

qresl = zeros( rows, neig );          % allocate space for residuals
qresr = zeros( rows, neig );

qresl_norms = zeros( 1, neig );       % space for residual norms
qresr_norms = zeros( 1, neig );

for index = 1:neig
  lambda = ritz( index );

  x = evr( :, index );                % right eigenvector
  qresr( :, index ) = ( lambda * lambda * M + lambda * C + K ) * x;

  qresr_norms( index ) = norm( qresr( :, index ) ) / ...
                         ( matrix_norm * norm( x ) );

  y = evl( :, index );
  qresl( :, index ) = ( y' * ( lambda * lambda * M + lambda * C + K ) )';

  qresl_norms( index ) = norm( qresl( :, index ) ) / ...
                         ( matrix_norm * norm( y ) );
end

% plot eigenvalues
% ================

subplot( 2, 1, 1 );
hold on;

title( 'eigenvalues' );

% clip axes so only region of interest appears
minx = min( [ real(alpha) real( ritz( 1:neig )' ) ] );
maxx = max( [ real(alpha) real( ritz( 1:neig )' ) ] );
miny = min( [ imag(alpha) imag( ritz( 1:neig )' ) ] );
maxy = max( [ imag(alpha) imag( ritz( 1:neig )' ) ] );

offx = ( maxx - minx ) / 8.0;
offy = ( maxy - miny ) / 8.0;

axis( [ minx-offx maxx+offx miny-offy maxy+offy ] );

% actual eigenvalues
act_eig = eig( full(A), full(B) );
plot( real( act_eig ), imag( act_eig ), 'ro' );

% approximate eigenvalues
plot( real( ritz( 1:neig ) ), imag( ritz( 1:neig ) ), 'b+' );

% plot shift value alpha
fprintf( 'alpha: %f + %fi\n', real(alpha), imag(alpha) );
plot( real( alpha ), imag( alpha ), 'g*' );

hold off;

fprintf( 'shift value:  green *\n'  );
fprintf( 'eigenvalues:  blue  +  (approximate)\n' );
fprintf( '              red   o  (actual)\n' );

% plot residual norms
% ------------------- 

subplot( 2, 1, 2 );
hold on;
title( 'residual norms' );
semilogy( qresl_norms, 'gx' );
semilogy( qresr_norms, 'bx' );
hold off;

fprintf( 'residual norms:  left  - green x\n' );
fprintf( '                 right - blue  x\n' );

