[m,n,a,lda,pivot,ipvt,lipvt,sigma,acnorm,wa]=qrfac(m,n,a,lda,pivot,ipvt,lipvt,sigma,acnorm,wa);

function [m,n,a,lda,pivot,ipvt,lipvt,sigma,acnorm,wa]=qrfac(m,n,a,lda,pivot,ipvt,lipvt,sigma,acnorm,wa);
%***BEGIN PROLOGUE  QRFAC
%***SUBSIDIARY
%***PURPOSE  Subsidiary to SNLS1, SNLS1E, SNSQ and SNSQE
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (QRFAC-S, DQRFAC-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     This subroutine uses Householder transformations with column
%     pivoting (optional) to compute a QR factorization of the
%     M by N matrix A. That is, QRFAC determines an orthogonal
%     matrix Q, a permutation matrix P, and an upper trapezoidal
%     matrix R with diagonal elements of nonincreasing magnitude,
%     such that A*P = Q*R. The Householder transformation for
%     column K, K = 1,2,...,MIN(M,N), is of the form
%
%                           T
%           I - (1/U(K))*U*U
%
%     where U has zeros in the first K-1 positions. The form of
%     this transformation and the method of pivoting first
%     appeared in the corresponding LINPACK subroutine.
%
%     The subroutine statement is
%
%       subroutine QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,SIGMA,ACNORM,WA)
%
%     where
%
%       M is a positive integer input variable set to the number
%         of rows of A.
%
%       N is a positive integer input variable set to the number
%         of columns of A.
%
%       A is an M by N array. On input A contains the matrix for
%         which the QR factorization is to be computed. On output
%         the strict upper trapezoidal part of A contains the strict
%         upper trapezoidal part of R, and the lower trapezoidal
%         part of A contains a factored form of Q (the non-trivial
%         elements of the U vectors described above).
%
%       LDA is a positive integer input variable not less than M
%         which specifies the leading dimension of the array A.
%
%       PIVOT is a logical input variable. If pivot is set true,
%         then column pivoting is enforced. If pivot is set false,
%         then no column pivoting is done.
%
%       IPVT is an integer output array of length LIPVT. IPVT
%         defines the permutation matrix P such that A*P = Q*R.
%         Column J of P is column IPVT(J) of the identity matrix.
%         If pivot is false, IPVT is not referenced.
%
%       LIPVT is a positive integer input variable. If PIVOT is
%             false, then LIPVT may be as small as 1. If PIVOT is
%             true, then LIPVT must be at least N.
%
%       SIGMA is an output array of length N which contains the
%         diagonal elements of R.
%
%       ACNORM is an output array of length N which contains the
%         norms of the corresponding columns of the input matrix A.
%         If this information is not needed, then ACNORM can coincide
%         with SIGMA.
%
%       WA is a work array of length N. If pivot is false, then WA
%         can coincide with SIGMA.
%
%***SEE ALSO  SNLS1, SNLS1E, SNSQ, SNSQE
%***ROUTINES CALLED  ENORM, R1MACH
%***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  QRFAC
persistent ajnorm epsmch firstCall i j jp1 k kmax minmn one p05 summlv temp zero ; if isempty(firstCall),firstCall=1;end; 

ipvt_shape=size(ipvt);ipvt=reshape(ipvt,1,[]);
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
sigma_shape=size(sigma);sigma=reshape(sigma,1,[]);
acnorm_shape=size(acnorm);acnorm=reshape(acnorm,1,[]);
wa_shape=size(wa);wa=reshape(wa,1,[]);
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(jp1), jp1=0; end;
if isempty(k), k=0; end;
if isempty(kmax), kmax=0; end;
if isempty(minmn), minmn=0; end;
if isempty(ajnorm), ajnorm=0; end;
if isempty(epsmch), epsmch=0; end;
if isempty(one), one=0; end;
if isempty(p05), p05=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,  p05 =[5.0e-2];  end;
if firstCall,  zero=[0.0e0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  QRFAC
[epsmch ]=r1mach(4);
%
%     COMPUTE THE INITIAL COLUMN NORMS AND INITIALIZE SEVERAL ARRAYS.
%
for j = 1 : n;
[acnorm(j) ,m,a(sub2ind(size(a),1,j):end)]=enorm(m,a(sub2ind(size(a),1,j):end));
sigma(j) = acnorm(j);
wa(j) = sigma(j);
if( pivot )
ipvt(j) = fix(j);
end;
end; j = fix(n+1);
%
%     REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS.
%
minmn = fix(min(m,n));
for j = 1 : minmn;
if( pivot )
%
%        BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION.
%
kmax = fix(j);
for k = j : n;
if( sigma(k)>sigma(kmax) )
kmax = fix(k);
end;
end; k = fix(n+1);
if( kmax~=j )
for i = 1 : m;
temp = a(i,j);
a(i,j) = a(i,kmax);
a(i,kmax) = temp;
end; i = fix(m+1);
sigma(kmax) = sigma(j);
wa(kmax) = wa(j);
k = fix(ipvt(j));
ipvt(j) = fix(ipvt(kmax));
ipvt(kmax) = fix(k);
end;
end;
%
%        COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE
%        J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR.
%
[ajnorm ,dumvar2,a(sub2ind(size(a),j,j):end)]=enorm(m-j+1,a(sub2ind(size(a),j,j):end));
if( ajnorm~=zero )
if( a(j,j)<zero )
ajnorm = -ajnorm;
end;
for i = j : m;
a(i,j) = a(i,j)./ajnorm;
end; i = fix(m+1);
a(j,j) = a(j,j) + one;
%
%        APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS
%        AND UPDATE THE NORMS.
%
jp1 = fix(j + 1);
if( n>=jp1 )
for k = jp1 : n;
summlv = zero;
for i = j : m;
summlv = summlv + a(i,j).*a(i,k);
end; i = fix(m+1);
temp = summlv./a(j,j);
for i = j : m;
a(i,k) = a(i,k) - temp.*a(i,j);
end; i = fix(m+1);
if( ~(~pivot || sigma(k)==zero) )
temp = a(j,k)./sigma(k);
sigma(k) = sigma(k).*sqrt(max(zero,one-temp.^2));
if( p05.*(sigma(k)./wa(k)).^2<=epsmch )
[sigma(k) ,dumvar2,a(sub2ind(size(a),jp1,k):end)]=enorm(m-j,a(sub2ind(size(a),jp1,k):end));
wa(k) = sigma(k);
end;
end;
end; k = fix(n+1);
end;
end;
sigma(j) = -ajnorm;
end; j = fix(minmn+1);
%
%     LAST CARD OF SUBROUTINE QRFAC.
%
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
sigma_shape=zeros(sigma_shape);sigma_shape(:)=sigma(1:numel(sigma_shape));sigma=sigma_shape;
acnorm_shape=zeros(acnorm_shape);acnorm_shape(:)=acnorm(1:numel(acnorm_shape));acnorm=acnorm_shape;
wa_shape=zeros(wa_shape);wa_shape(:)=wa(1:numel(wa_shape));wa=wa_shape;
end
%DECK QRSOLV