KK = 1e3; 
nnodes = 29; % number of nodes of the bridge (odd)
[ A ] = bridge_stiffness_matrix( nnodes, KK );
plot_bridge(nnodes);
spy(A);

spy(A'-A);

K_A = condest( A )
%   K_A =
%      Inf

n = size( A, 2 );
nr = n - 3;
ir = [ 3 : n - 1 ];
Ar = A( ir, ir );
K_Ar = condest( Ar )
%   K_Ar =
%      1.8838e+04

node_force = 1 : 2 : nnodes;
b = zeros( n, 1 );
b( 2 * node_force, 1 ) = -1;
br = b( ir, 1 );
dr = Ar \ br;
d = zeros( n, 1 );
d( ir, 1 ) = dr;
[ fig_h ] = plot_bridge( nnodes, d );

% To run this part, you need to have the functions forward_substitutions
% and backward_substitutions in the Matlab path
[ Lr, Ur, Pr ] = lu( Ar );
for l = 1 : 10    
    b = 2 * ( rand( n, 1 ) - 0.5 );  % random force in [-1,1]
    br = b( ir, 1 );
    [ yr ] = forward_substitutions( Lr, Pr * br );
    [ dr ] = backward_substitutions( Ur, yr );
    d = zeros( n, 1 );
    d( ir, 1 ) = dr;
    [ fig_h ] = plot_bridge( nnodes, d );
    pause( 1 );
end

node_force = 1 : 2 : nnodes;
b = zeros( n, 1 );
b( 2 * node_force, 1 ) = -1;
br = b( ir, 1 );
dr_gmres = gmres( Ar, br, [], 5e-5, nr );
%   gmres converged at iteration 54 to a solution with relative residual 1.8e-05.
dr_cg = pcg( Ar, br, 5e-5, 1000, [] );
%   pcg converged at iteration 56 to a solution with relative residual 9.8e-06.