[n,kmp,ll,maxl,v,q,rl,snormw,prod,r0nrm]=drlcal(n,kmp,ll,maxl,v,q,rl,snormw,prod,r0nrm);

function [n,kmp,ll,maxl,v,q,rl,snormw,prod,r0nrm]=drlcal(n,kmp,ll,maxl,v,q,rl,snormw,prod,r0nrm);
%***BEGIN PROLOGUE  DRLCAL
%***SUBSIDIARY
%***PURPOSE  Internal routine for DGMRES.
%***LIBRARY   SLATEC (SLAP)
%***CATEGORY  D2A4, D2B4
%***TYPE      doubleprecision (SRLCAL-S, DRLCAL-D)
%***KEYWORDS  GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
%             NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
%***AUTHOR  Brown, Peter, (LLNL), [e-mail]
%           Hindmarsh, Alan, (LLNL), [e-mail]
%           Seager, Mark K., (LLNL), [e-mail]
%             Lawrence Livermore National Laboratory
%             PO Box 808, L-60
%             Livermore, CA 94550 [phone]
%***DESCRIPTION
%         This routine calculates the scaled residual RL from the
%         V(I)'s.
% *Usage:
%      INTEGER N, KMP, LL, MAXL
%      doubleprecision V(N,LL), Q(2*MAXL), RL(N), SNORMW, PROD, R0NORM
%
%      CALL DRLCAL(N, KMP, LL, MAXL, V, Q, RL, SNORMW, PROD, R0NRM)
%
% *Arguments:
% N      :IN       Integer
%         The order of the matrix A, and the lengths
%         of the vectors SR, SZ, R0 and Z.
% KMP    :IN       Integer
%         The number of previous V vectors the new vector VNEW
%         must be made orthogonal to. (KMP <= MAXL)
% LL     :IN       Integer
%         The current dimension of the Krylov subspace.
% MAXL   :IN       Integer
%         The maximum dimension of the Krylov subspace.
% V      :IN       doubleprecision V(N,LL)
%         The N x LL array containing the orthogonal vectors
%         V(*,1) to V(*,LL).
% Q      :IN       doubleprecision Q(2*MAXL)
%         A doubleprecision array of length 2*MAXL containing the
%         components of the Givens rotations used in the QR
%         decomposition of HES.  It is loaded in DHEQR and used in
%         DHELS.
% RL     :OUT      doubleprecision RL(N)
%         The residual vector RL.  This is either SB*(B-A*XL) if
%         not preconditioning or preconditioning on the right,
%         or SB*(M-inverse)*(B-A*XL) if preconditioning on the
%         left.
% SNORMW :IN       doubleprecision
%         Scale factor.
% PROD   :IN       doubleprecision
%         The product s1*s2*...*sl = the product of the sines of the
%         Givens rotations used in the QR factorization of
%         the Hessenberg matrix HES.
% R0NRM  :IN       doubleprecision
%         The scaled norm of initial residual R0.
%
%***SEE ALSO  DGMRES
%***ROUTINES CALLED  DCOPY, DSCAL
%***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  DRLCAL
%         The following is for optimized compilation on LLNL/LTSS Crays.
%LLL. OPTIMIZE
%     .. Scalar Arguments ..
%     .. Array Arguments ..
persistent c i i2 ip1 k llm1 llp1 s tem ; 

q_shape=size(q);q=reshape(q,1,[]);
v_shape=size(v);v=reshape([v(:).',zeros(1,ceil(numel(v)./prod([n])).*prod([n])-numel(v))],n,[]);
%     .. Local Scalars ..
if isempty(c), c=0; end;
if isempty(s), s=0; end;
if isempty(tem), tem=0; end;
if isempty(i), i=0; end;
if isempty(i2), i2=0; end;
if isempty(ip1), ip1=0; end;
if isempty(k), k=0; end;
if isempty(llm1), llm1=0; end;
if isempty(llp1), llp1=0; end;
%     .. External Subroutines ..
%***FIRST EXECUTABLE STATEMENT  DRLCAL
if( kmp==maxl )
%
%         calculate RL.  Start by copying V(*,1) into RL.
%
[n,v(sub2ind(size(v),1,1):end),dumvar3,rl]=dcopy(n,v(sub2ind(size(v),1,1):end),1,rl,1);
llm1 = fix(ll - 1);
for i = 1 : llm1;
ip1 = fix(i + 1);
i2 = fix(i.*2);
s = q(i2);
c = q(i2-1);
for k = 1 : n;
rl(k) = s.*rl(k) + c.*v(k,ip1);
end; k = fix(n+1);
end; i = fix(llm1+1);
s = q(2.*ll);
c = q(2.*ll-1)./snormw;
llp1 = fix(ll + 1);
for k = 1 : n;
rl(k) = s.*rl(k) + c.*v(k,llp1);
end; k = fix(n+1);
end;
%
%         When KMP < MAXL, RL vector already partially calculated.
%         Scale RL by R0NRM*PROD to obtain the residual RL.
%
tem = r0nrm.*prod;
[n,tem,rl]=dscal(n,tem,rl,1);
%------------- LAST LINE OF DRLCAL FOLLOWS ----------------------------
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
end
%DECK DROT