function B = jacobimul(A,i,j,c,s)
%jacobimul Computes the Jacobi rotation for an iteration
%of the Jacobi algorithm that approximates the eigenvalues
%of a real symmetric matrix.
%
%   B = jacobimul(A,i,j,c,s) computes the product J(i,j,c,s)'*A*J(i,j,c,s),
%   where A is an n x n real symmetric matrix and c, s, i, j are
%   Givens matrix parameters. Upon return, B = J'*A*J.
   
%use givensmul to compute J(i,j,c,s)'*A
B = givensmul(A,i,j,c,-s);
%now form A*J(i,j,c,s)
% a = B(:,i);
% b = B(:,j);
% B(:,i) = c*a - s*b;
% B(:,j) = s*a + c*b;
B = givensmulp(B,i,j,c,s);
