mgsr.m

UTV工具包提供46个Matlab函数

function q = mgsr(U,kappa)
%  mgsr --> Modified Gram-Schmidt with re-orthogonalization (expansion step).
%
%  <Synopsis>
%    q = mgsr(U,kappa)
%
%  <Description>
%    Modified Gram-Schmidt with re-orthogonalization is used to expand
%    the m-by-n matrix U having orthogonal columns with a new column q,
%    which is orthogonal to the other columns of U. The parameter kappa
%    (greater than one) is used to decide if re-orthogonalization is needed;
%    kappa = 1 ensures re-orthogonalization, however, a typical value for
%    kappa is sqrt(2).
%
%  <Algorithm>
%    The algorithm relies on the fact that one re-orthogonalization
%    is always enough.
%
%  <See Also>
%    ulv_csne  --> Corrected semi-normal equations expansion (ULV).
%    urv_csne  --> Corrected semi-normal equations expansion (URV).
%    ullv_csne --> Corrected semi-normal equations expansion (ULLV).

%  <References>
%  [1] J.W. Daniel, W.B. Gragg, L. Kaufman and G.W. Stewart,
%      "Reorthogonalization and Stable Algorithms for Updating the
%      Gram-Schmidt QR Factorization", Math. Comp., 30 (1976), pp. 772--795.
%
%  <Revision>
%    Ricardo D. Fierro, California State University San Marcos
%    Per Christian Hansen, IMM, Technical University of Denmark
%    Peter S.K. Hansen, IMM, Technical University of Denmark
%
%    Last revised: June 22, 1999
%-----------------------------------------------------------------------

[m,n] = size(U);

% Vector e1.
q = eye(m,1);

% Use Modified Gram-Schmidt to orthogonalize e1 against the columns of U.
for (k=1:n)
  r_kn = U(:,k)'*q;
  q    = q - U(:,k)*r_kn;
end

% 2-norm of first solution q.
q1norm = norm(q);

% Re-orthogonalization step.
if (q1norm <= 1.0/kappa)
  for (k=1:n)
    r_kn = U(:,k)'*q;
    q    = q - U(:,k)*r_kn;
  end

  % 2-norm of second solution q.
  q2norm = norm(q);

  % If e1 lies in the range of U, then we can orthogonalize any vector to U.
  if (q2norm <= q1norm/kappa)
    q = [1:m]'; q = q/norm(q);

    for (k=1:n)
      r_kn = U(:,k)'*q;
      q    = q - U(:,k)*r_kn;
    end

    % 2-norm of third solution q.
    q3norm = norm(q);

    % Re-orthogonalization step.
    if (q3norm <= 1.0/kappa)
      for (k=1:n)
        r_kn = U(:,k)'*q;
        q    = q - U(:,k)*r_kn;
      end

      % 2-norm of 4th solution q.
      q4norm = norm(q);

      % Normalization of the expanded column of U.
      q = q/q4norm;
    else
      % Normalization of the expanded column of U.
      q = q/q3norm;
    end
  else
    % Normalization of the expanded column of U.
    q = q/q2norm;
  end
else
  % Normalization of the expanded column of U.
  q = q/q1norm;
end

%-----------------------------------------------------------------------
% End of function mgsr
%-----------------------------------------------------------------------