[m,n,q,ldq,wa]=qform(m,n,q,ldq,wa);

function [m,n,q,ldq,wa]=qform(m,n,q,ldq,wa);
%***BEGIN PROLOGUE  QFORM
%***SUBSIDIARY
%***PURPOSE  Subsidiary to SNSQ and SNSQE
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (QFORM-S, DQFORM-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     This subroutine proceeds from the computed QR factorization of
%     an M by N matrix A to accumulate the M by M orthogonal matrix
%     Q from its factored form.
%
%     The subroutine statement is
%
%       subroutine QFORM(M,N,Q,LDQ,WA)
%
%     where
%
%       M is a positive integer input variable set to the number
%         of rows of A and the order of Q.
%
%       N is a positive integer input variable set to the number
%         of columns of A.
%
%       Q is an M by M array. On input the full lower trapezoid in
%         the first min(M,N) columns of Q contains the factored form.
%         On output Q has been accumulated into a square matrix.
%
%       LDQ is a positive integer input variable not less than M
%         which specifies the leading dimension of the array Q.
%
%       WA is a work array of length M.
%
%***SEE ALSO  SNSQ, SNSQE
%***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  QFORM
persistent firstCall i j jm1 k l minmn np1 one summlv temp zero ; if isempty(firstCall),firstCall=1;end; 

q_shape=size(q);q=reshape([q(:).',zeros(1,ceil(numel(q)./prod([ldq])).*prod([ldq])-numel(q))],ldq,[]);
wa_shape=size(wa);wa=reshape(wa,1,[]);
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(jm1), jm1=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(minmn), minmn=0; end;
if isempty(np1), np1=0; end;
if isempty(one), one=0; end;
if isempty(summlv), summlv=0; end;
if isempty(temp), temp=0; end;
if isempty(zero), zero=0; end;
if firstCall,   one =[1.0e0];  end;
if firstCall,  zero=[0.0e0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  QFORM
minmn = fix(min(m,n));
if( minmn>=2 )
for j = 2 : minmn;
jm1 = fix(j - 1);
for i = 1 : jm1;
q(i,j) = zero;
end; i = fix(jm1+1);
end; j = fix(minmn+1);
end;
%
%     INITIALIZE REMAINING COLUMNS TO THOSE OF THE IDENTITY MATRIX.
%
np1 = fix(n + 1);
if( m>=np1 )
for j = np1 : m;
for i = 1 : m;
q(i,j) = zero;
end; i = fix(m+1);
q(j,j) = one;
end; j = fix(m+1);
end;
%
%     ACCUMULATE Q FROM ITS FACTORED FORM.
%
for l = 1 : minmn;
k = fix(minmn - l + 1);
for i = k : m;
wa(i) = q(i,k);
q(i,k) = zero;
end; i = fix(m+1);
q(k,k) = one;
if( wa(k)~=zero )
for j = k : m;
summlv = zero;
for i = k : m;
summlv = summlv + q(i,j).*wa(i);
end; i = fix(m+1);
temp = summlv./wa(k);
for i = k : m;
q(i,j) = q(i,j) - temp.*wa(i);
end; i = fix(m+1);
end; j = fix(m+1);
end;
end; l = fix(minmn+1);
%
%     LAST CARD OF SUBROUTINE QFORM.
%
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
wa_shape=zeros(wa_shape);wa_shape(:)=wa(1:numel(wa_shape));wa=wa_shape;
end
%DECK QK15