[n,r,ldr,w,b,alpha,cosmlv,sinmlv]=dwupdt(n,r,ldr,w,b,alpha,cosmlv,sinmlv);

function [n,r,ldr,w,b,alpha,cosmlv,sinmlv]=dwupdt(n,r,ldr,w,b,alpha,cosmlv,sinmlv);
%***BEGIN PROLOGUE  DWUPDT
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DNLS1 and DNLS1E
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (RWUPDT-S, DWUPDT-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     Given an N by N upper triangular matrix R, this subroutine
%     computes the QR decomposition of the matrix formed when a row
%     is added to R. If the row is specified by the vector W, then
%     DWUPDT determines an orthogonal matrix Q such that when the
%     N+1 by N matrix composed of R augmented by W is premultiplied
%     by (Q TRANSPOSE), the resulting matrix is upper trapezoidal.
%     The orthogonal matrix Q is the product of N transformations
%
%           G(1)*G(2)* ... *G(N)
%
%     where G(I) is a Givens rotation in the (I,N+1) plane which
%     eliminates elements in the I-th plane. DWUPDT also
%     computes the product (Q TRANSPOSE)*C where C is the
%     (N+1)-vector (b,alpha). Q itself is not accumulated, rather
%     the information to recover the G rotations is supplied.
%
%     The subroutine statement is
%
%       subroutine DWUPDT(N,R,LDR,W,B,ALPHA,COS,SIN)
%
%     where
%
%       N is a positive integer input variable set to the order of R.
%
%       R is an N by N array. On input the upper triangular part of
%         R must contain the matrix to be updated. On output R
%         contains the updated triangular matrix.
%
%       LDR is a positive integer input variable not less than N
%         which specifies the leading dimension of the array R.
%
%       W is an input array of length N which must contain the row
%         vector to be added to R.
%
%       B is an array of length N. On input B must contain the
%         first N elements of the vector C. On output B contains
%         the first N elements of the vector (Q TRANSPOSE)*C.
%
%       ALPHA is a variable. On input ALPHA must contain the
%         (N+1)-st element of the vector C. On output ALPHA contains
%         the (N+1)-st element of the vector (Q TRANSPOSE)*C.
%
%       COS is an output array of length N which contains the
%         cosines of the transforming Givens rotations.
%
%       SIN is an output array of length N which contains the
%         sines of the transforming Givens rotations.
%
%     **********
%
%***SEE ALSO  DNLS1, DNLS1E
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   800301  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  DWUPDT
persistent cotanmlv firstCall i j jm1 one p25 p5 rowj tanmlv temp zero ; if isempty(firstCall),firstCall=1;end; 

r_shape=size(r);r=reshape([r(:).',zeros(1,ceil(numel(r)./prod([ldr])).*prod([ldr])-numel(r))],ldr,[]);
w_shape=size(w);w=reshape(w,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
cos_shape=size(cosmlv);cosmlv=reshape(cosmlv,1,[]);
sin_shape=size(sinmlv);sinmlv=reshape(sinmlv,1,[]);
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(jm1), jm1=0; end;
if isempty(cotanmlv), cotanmlv=0; end;
if isempty(one), one=0; end;
if isempty(p5), p5=0; end;
if isempty(p25), p25=0; end;
if isempty(rowj), rowj=0; end;
if isempty(tanmlv), tanmlv=0; end;
if isempty(temp), temp=0; end;
if isempty(zero), zero=0; end;
if firstCall,   one =[1.0d0];  end;
if firstCall,  p5 =[5.0d-1];  end;
if firstCall,  p25 =[2.5d-1];  end;
if firstCall,  zero=[0.0d0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DWUPDT
for j = 1 : n;
rowj = w(j);
jm1 = fix(j - 1);
%
%        APPLY THE PREVIOUS TRANSFORMATIONS TO
%        R(I,J), I=1,2,...,J-1, AND TO W(J).
%
if( jm1>=1 )
for i = 1 : jm1;
temp = cosmlv(i).*r(i,j) + sinmlv(i).*rowj;
rowj = -sinmlv(i).*r(i,j) + cosmlv(i).*rowj;
r(i,j) = temp;
end; i = fix(jm1+1);
end;
%
%        DETERMINE A GIVENS ROTATION WHICH ELIMINATES W(J).
%
cosmlv(j) = one;
sinmlv(j) = zero;
if( rowj~=zero )
if( abs(r(j,j))>=abs(rowj) )
tanmlv = rowj./r(j,j);
cosmlv(j) = p5./sqrt(p25+p25.*tanmlv.^2);
sinmlv(j) = cosmlv(j).*tanmlv;
else;
cotanmlv = r(j,j)./rowj;
sinmlv(j) = p5./sqrt(p25+p25.*cotanmlv.^2);
cosmlv(j) = sinmlv(j).*cotanmlv;
end;
%
%        APPLY THE CURRENT TRANSFORMATION TO R(J,J), B(J), AND ALPHA.
%
r(j,j) = cosmlv(j).*r(j,j) + sinmlv(j).*rowj;
temp = cosmlv(j).*b(j) + sinmlv(j).*alpha;
alpha = -sinmlv(j).*b(j) + cosmlv(j).*alpha;
b(j) = temp;
end;
end; j = fix(n+1);
%
%     LAST CARD OF SUBROUTINE DWUPDT.
%
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
cos_shape=zeros(cos_shape);cos_shape(:)=cosmlv(1:numel(cos_shape));cosmlv=cos_shape;
sin_shape=zeros(sin_shape);sin_shape(:)=sinmlv(1:numel(sin_shape));sinmlv=sin_shape;
end
%DECK DX4