% Solution to Problem 4.6
%
% Inputs:
%   ( d = dimension of matrix F = min(size(NumTerms)) )
%   NumTerms(d,d), where
%     NumTerms(i,j) = number of polynomials terms in F(i,j) defining Surface S
%   ( tS = total number of polynomials terms in all entries of F 
%        = sum(sum(NumTerms)) )
%   Sterms(4,tS) 
%     All terms c*x1^e1*x2^e2*x3^e3 stored sequentially, columnwise for F
%     Sterms(1,k) = exponent e1 for term k
%     Sterms(2,k) = exponent e2 for term k
%     Sterms(3,k) = exponent e3 for term k
%     Sterms(4,k) = coefficient c for term k
%   tC(3) = degree of polynomials defining components x1, x2, x3 of Curve C
%   Curve(sum(tC)+3) = coefficients of polynomials, x1 then x2 then x3,
%   from constant term up
%
d = min(size(NumTerms));
max_degree = max(tC*Sterms(1:3,:))
% Initialize array [A_0, A_1, ... , A_d] of Matrix Polynomial Coefficients
Matrix_Poly_Coeffs = zeros(d,(max_degree+1)*d);
% Initialize array of Matrix Polynomial Coefficients error bounds
dMatrix_Poly_Coeffs = zeros(d,(max_degree+1)*d);
k=1;
% For each entry of F
for j=1:d,
   for i=1:d
%     Find polynomial terms in Sterms
      terms = (k : k + NumTerms(i,j) -1);
%     Compute degree of F(i,j)
      degij = max( tC * Sterms(1:3,terms) );
%     Initialize vector of coefficients of F(i,j)
      coeffij = zeros( 1, degij + 1 );
%     Initialize vector of coefficient error bounds for F(i,j)
      dcoeffij = zeros( 1, degij + 1 );
%     For each term in polynomial p defining F(i,j)
      for r = terms,
%        Compute term and add to coeffij
%        Compute error bound on term and add to dcoeffij
         p = Sterms(4,r);
         dp = 0;
%        p = p * x1^e1
         x1 = Curve( 1 : tC(1)+1 );
         for t = 1 : Sterms(1,r),
            dp = conv(dp, abs(x1) ) + eps*conv( abs(p), abs(x1) );
             p = conv( p, x1 );
         end
%        p = p * x2^e2
         x2 = Curve( tC(1)+2 : tC(1)+tC(2)+2 );
         for t = 1 : Sterms(2,r),
            dp = conv(dp, abs(x2) ) + eps*conv( abs(p), abs(x2) );
             p = conv( p, x2 );
         end
%        p = p * x3^e3
         x3 = Curve( tC(1)+tC(2)+3 : tC(1)+tC(2)+tC(3)+3 );
         for t = 1 : Sterms(3,r),
            dp = conv(dp, abs(x3) ) + eps*conv( abs(p), abs(x3) );
             p = conv(p, x3 );
         end
%        Accumulate p into array of coefficients of F(i,j)
         coeffij = coeffij + [ p , zeros( 1, length(coeffij)-length(p) ) ];
%        Accumulate dp into array of error bounds on coefficients of F(i,j)
         dcoeffij = dcoeffij + ...
                    [ dp , zeros( 1, length(coeffij)-length(p) ) ] + ...
                    eps*abs(coeffij);
      end
%     Store F(i,j) in Matrix_Poly_Coeffs 
      loc = ( j : d : degij*d + j );
      Matrix_Poly_Coeffs( i , loc ) = coeffij;
%     Store error bound on F(i,j) in dMatrix_Poly_Coeffs 
      dMatrix_Poly_Coeffs( i , loc ) = dcoeffij;
      k = k + NumTerms(i,j);
   end
end
Matrix_Poly_Coeffs     % display matrix polynomial coefficients
disp('pause (hit enter to continue)'), pause
dMatrix_Poly_Coeffs    % display error bounds on Matrix_Poly_Coeffs
disp('pause (hit enter to continue)'), pause
dim = d*max_degree;
%
% Build matrix pencil A - lambda B whose eigenvalues we want
%
B = eye(dim);
B(1:d,1:d) = Matrix_Poly_Coeffs( : , dim+1 : dim+d );
A = zeros(dim,dim);
A(d+1:dim,1:dim-d) = eye(dim-d);
for i = 1:max_degree
  A(1:d,(i-1)*d+1:i*d) = ...
       -Matrix_Poly_Coeffs( 1:d , dim - i*d + 1 : dim+ (1-i)*d);
end
% Display pencil
A 
disp('pause (hit enter to continue)'), pause
B
disp('pause (hit enter to continue)'), pause
%
% If B(1:d,1:d) is well conditioned, convert to standard eigenproblem
% and compute error bounds
%
rcB = rcond(B(1:d,1:d));
if (rcB > 1e-8)
%    Convert to standard eigenproblem
     disp(['1/cond(A_d) = ', num2str(rcB), ' so convert to ' ...
     'standard eigenproblem for'])
     A(1:d,:) = (B(1:d,1:d))\A(1:d,:);
     A
     disp('pause (hit enter to continue)'), pause
     f1 = flops;
     D = eig(A)
     f2 = flops - f1;
     disp(['cost( D=eig(A) ) = ', num2str(f2/dim^3), ...
          ' n^3 flops, with n = ',int2str(dim)])
%    Compute error bounds for computed eigenvalues
%    Get right and left eigenvectors (Vr and Vl)
     f1 = flops;
     [Vr,Dr] = eig(A);
     [Vl,Dl] = eig(A');
     f2 = flops - f1;
     disp(['cost( left,right evectors ) = ', num2str(f2/dim^3), ...
          ' n^3 flops, with n = ',int2str(dim)])
%    Reorder eigenvalues, eigenvectors so they match
     [sDr,sIDr] = sort(diag(Dr));
     Vr = Vr(:,sIDr);
     D = sDr;
     [sDl,sIDl] = sort(diag(Dl));
     Vl = Vl(:,sIDl);
%    Compute eigenvalue condition numbers
     Econd = 1 ./ abs(diag(conj(Vl')*Vr));
%    Compute eigenvalue error bounds
     n1 = norm( dMatrix_Poly_Coeffs( : , d+1:d+dim ), 1 );
     n2 = norm( dMatrix_Poly_Coeffs( : , 1:d ), 1 );
     n3 = norm( A( 1:d, : ), 1 );
     Ebnd = Econd .*  ( n1/rcB + n2/rcB * n3 );
     disp(' Eigenvalue                  Error Bound' )
     [D, Ebnd]
     disp('pause (hit enter to continue)'), pause
%    Display intersection points
%    Get coeffs of polynomials and derivatives of polynomials
     polyx1  = fliplr(Curve( 1 : tC(1)+1 ));
     dpolyx1 = (tC(1):-1:1) .* polyx1(1:tC(1));
     polyx2 = fliplr(Curve( tC(1)+2 : tC(1)+tC(2)+2 ));
     dpolyx2 = (tC(2):-1:1) .* polyx2(1:tC(2));
     polyx3 = fliplr(Curve( tC(1)+tC(2)+3 : tC(1)+tC(2)+tC(3)+3 ));
     dpolyx3 = (tC(3):-1:1) .* polyx3(1:tC(3));
%    Evaluate polynomials and error bounds
     Cx1 = polyval( polyx1 , D );
     Cx2 = polyval( polyx2 , D );
     Cx3 = polyval( polyx3 , D );
     dCx1 = abs(polyval( dpolyx1 , D )) .* Ebnd + ...
            eps*polyval( abs(polyx1), abs(D) );
     dCx2 = abs(polyval( dpolyx2 , D )) .* Ebnd + ...
            eps*polyval( abs(polyx2), abs(D) );
     dCx3 = abs(polyval( dpolyx3 , D )) .* Ebnd + ...
            eps*polyval( abs(polyx3), abs(D) );
%    Display results for all eigenvalues, including complex and infinite
     disp('all eigenvalues, x1 coords of intersections, error bounds')
     [D,Cx1,dCx1]
     disp('all eigenvalues, x2 coords of intersections, error bounds')
     [D,Cx2,dCx2]
     disp('all eigenvalues, x3 coords of intersections, error bounds')
     [D,Cx3,dCx3]
     disp('pause (hit enter to continue)'), pause
%    Display results for real eigenvalues only
     Dreal = find(abs(imag(D)) < 1e-11*abs(real(D)));
     disp('real eigenvalues, x1 coords of intersections, error bounds')
     [D(Dreal),Cx1(Dreal),dCx1(Dreal)]
     disp('real eigenvalues, x2 coords of intersections, error bounds')
     [D(Dreal),Cx2(Dreal),dCx2(Dreal)]
     disp('real eigenvalues, x3 coords of intersections, error bounds')
     [D(Dreal),Cx3(Dreal),dCx3(Dreal)]
else
%    Solve as generalized eigenproblem (no error bound computed)
     disp(['1/cond(A_d) = ', num2str(rcB), ' so solve as ', ...
     'generalized eigenproblem'])
%    Display flop count
     f1 = flops;
     D=eig(B,A);
     f2 = flops - f1;
     D = sort(D);
     disp(['cost( D=eig(B,A) ) = ', num2str(f2/dim^3), ...
          ' n^3 flops, with n = ',int2str(dim)])
     disp('pause (hit enter to continue)'), pause
%    Display intersection points
     polyx1  = fliplr(Curve( 1 : tC(1)+1 ));
     polyx2 = fliplr(Curve( tC(1)+2 : tC(1)+tC(2)+2 ));
     polyx3 = fliplr(Curve( tC(1)+tC(2)+3 : tC(1)+tC(2)+tC(3)+3 ));
     Cx1 = polyval( polyx1 , D );
     Cx2 = polyval( polyx2 , D );
     Cx3 = polyval( polyx3 , D );
%    Display results for all eigenvalues, including complex and infinite
     disp('all eigenvalues, x1 coords of intersections')
     [D,Cx1]
     disp('all eigenvalues, x2 coords of intersections')
     [D,Cx2]
     disp('all eigenvalues, x3 coords of intersections')
     [D,Cx3]
     disp('pause (hit enter to continue)'), pause
%    Display results for real eigenvalues only
     Dreal = find(abs(imag(D)) <= 1e-11*abs(real(D)));
     disp('real eigenvalues, all coords of intersections')
     [D(Dreal),Cx1(Dreal),Cx2(Dreal),Cx3(Dreal)]
end
