| [vnew,v,hes,n,ll,ldhes,kmp,snormw]=dorth(vnew,v,hes,n,ll,ldhes,kmp,snormw); |
function [vnew,v,hes,n,ll,ldhes,kmp,snormw]=dorth(vnew,v,hes,n,ll,ldhes,kmp,snormw);
%***BEGIN PROLOGUE DORTH
%***SUBSIDIARY
%***PURPOSE Internal routine for DGMRES.
%***LIBRARY SLATEC (SLAP)
%***CATEGORY D2A4, D2B4
%***TYPE doubleprecision (SORTH-S, DORTH-D)
%***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
% NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
%***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
% Hindmarsh, Alan, (LLNL), alanh@llnl.gov
% Seager, Mark K., (LLNL), seager@llnl.gov
% Lawrence Livermore National Laboratory
% PO Box 808, L-60
% Livermore, CA 94550 (510) 423-3141
%***DESCRIPTION
% This routine orthogonalizes the vector VNEW against the
% previous KMP vectors in the V array. It uses a modified
% Gram-Schmidt orthogonalization procedure with conditional
% reorthogonalization.
%
% *Usage:
% INTEGER N, LL, LDHES, KMP
% doubleprecision VNEW(N), V(N,LL), HES(LDHES,LL), SNORMW
%
% CALL DORTH(VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
%
% *Arguments:
% VNEW :INOUT doubleprecision VNEW(N)
% On input, the vector of length N containing a scaled
% product of the Jacobian and the vector V(*,LL).
% On output, the new vector orthogonal to V(*,i0) to V(*,LL),
% where i0 = max(1, LL-KMP+1).
% V :IN doubleprecision V(N,LL)
% The N x LL array containing the previous LL
% orthogonal vectors V(*,1) to V(*,LL).
% HES :INOUT doubleprecision HES(LDHES,LL)
% On input, an LL x LL upper Hessenberg matrix containing,
% in HES(I,K), K<LL, the scaled inner products of
% A*V(*,K) and V(*,i).
% On return, column LL of HES is filled in with
% the scaled inner products of A*V(*,LL) and V(*,i).
% N :IN Integer
% The order of the matrix A, and the length of VNEW.
% LL :IN Integer
% The current order of the matrix HES.
% LDHES :IN Integer
% The leading dimension of the HES array.
% KMP :IN Integer
% The number of previous vectors the new vector VNEW
% must be made orthogonal to (KMP <= MAXL).
% SNORMW :OUT doubleprecision
% Scalar containing the l-2 norm of VNEW.
%
%***SEE ALSO DGMRES
%***ROUTINES CALLED DAXPY, DDOT, DNRM2
%***REVISION HISTORY (YYMMDD)
% 890404 DATE WRITTEN
% 890404 Previous REVISION DATE
% 890915 Made changes requested at July 1989 CML Meeting. (MKS)
% 890922 Numerous changes to prologue to make closer to SLATEC
% standard. (FNF)
% 890929 Numerous changes to reduce SP/DP differences. (FNF)
% 910411 Prologue converted to Version 4.0 format. (BAB)
% 910506 Made subsidiary to DGMRES. (FNF)
% 920511 Added complete declaration section. (WRB)
%***end PROLOGUE DORTH
% The following is for optimized compilation on LLNL/LTSS Crays.
%LLL. OPTIMIZE
% .. Scalar Arguments ..
% .. Array Arguments ..
persistent arg i i0 sumdsq tem vnrm ;
hes_shape=size(hes);hes=reshape([hes(:).',zeros(1,ceil(numel(hes)./prod([ldhes])).*prod([ldhes])-numel(hes))],ldhes,[]);
v_shape=size(v);v=reshape([v(:).',zeros(1,ceil(numel(v)./prod([n])).*prod([n])-numel(v))],n,[]);
vnew_shape=size(vnew);vnew=reshape(vnew,1,[]);
% .. Local Scalars ..
if isempty(arg), arg=0; end;
if isempty(sumdsq), sumdsq=0; end;
if isempty(tem), tem=0; end;
if isempty(vnrm), vnrm=0; end;
if isempty(i), i=0; end;
if isempty(i0), i0=0; end;
% .. External Functions ..
% .. External Subroutines ..
% .. Intrinsic Functions ..
%***FIRST EXECUTABLE STATEMENT DORTH
%
% Get norm of unaltered VNEW for later use.
%
[vnrm ,n,vnew]=dnrm2(n,vnew,1);
% -------------------------------------------------------------------
% Perform the modified Gram-Schmidt procedure on VNEW =A*V(LL).
% Scaled inner products give new column of HES.
% Projections of earlier vectors are subtracted from VNEW.
% -------------------------------------------------------------------
i0 = fix(max(1,ll-kmp+1));
for i = i0 : ll;
[hes(i,ll) ,n,v(sub2ind(size(v),1,i):end),dumvar4,vnew]=ddot(n,v(sub2ind(size(v),1,i):end),1,vnew,1);
tem = -hes(i,ll);
[n,tem,v(sub2ind(size(v),1,i):end),dumvar4,vnew]=daxpy(n,tem,v(sub2ind(size(v),1,i):end),1,vnew,1);
end; i = fix(ll+1);
% -------------------------------------------------------------------
% Compute SNORMW = norm of VNEW. If VNEW is small compared
% to its input value (in norm), then reorthogonalize VNEW to
% V(*,1) through V(*,LL). Correct if relative correction
% exceeds 1000*(unit roundoff). Finally, correct SNORMW using
% the dot products involved.
% -------------------------------------------------------------------
[snormw ,n,vnew]=dnrm2(n,vnew,1);
if( vnrm+0.001d0.*snormw~=vnrm )
hes_shape=zeros(hes_shape);hes_shape(:)=hes(1:numel(hes_shape));hes=hes_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
vnew_shape=zeros(vnew_shape);vnew_shape(:)=vnew(1:numel(vnew_shape));vnew=vnew_shape;
return;
end;
sumdsq = 0;
for i = i0 : ll;
tem = -ddot(n,v(sub2ind(size(v),1,i):end),1,vnew,1);
if( hes(i,ll)+0.001d0.*tem~=hes(i,ll) )
hes(i,ll) = hes(i,ll) - tem;
[n,tem,v(sub2ind(size(v),1,i):end),dumvar4,vnew]=daxpy(n,tem,v(sub2ind(size(v),1,i):end),1,vnew,1);
sumdsq = sumdsq + tem.^2;
end;
end; i = fix(ll+1);
if( sumdsq==0.0d0 )
hes_shape=zeros(hes_shape);hes_shape(:)=hes(1:numel(hes_shape));hes=hes_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
vnew_shape=zeros(vnew_shape);vnew_shape(:)=vnew(1:numel(vnew_shape));vnew=vnew_shape;
return;
end;
arg = max(0.0d0,snormw.^2-sumdsq);
snormw = sqrt(arg);
%
%------------- LAST LINE OF DORTH FOLLOWS ----------------------------
hes_shape=zeros(hes_shape);hes_shape(:)=hes(1:numel(hes_shape));hes=hes_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
vnew_shape=zeros(vnew_shape);vnew_shape(:)=vnew(1:numel(vnew_shape));vnew=vnew_shape;
end
%DECK DORTHR
|
|
