[x,ldx,n,p,qraux,jpvt,work,job]=dqrdc(x,ldx,n,p,qraux,jpvt,work,job);

function [x,ldx,n,p,qraux,jpvt,work,job]=dqrdc(x,ldx,n,p,qraux,jpvt,work,job);
persistent j jj jp l lp1 lup maxj maxnrm negj nrmxl pl pu swapj t tt ; 

if isempty(jj), jj=0; end;
%***BEGIN PROLOGUE  DQRDC
%***PURPOSE  Use Householder transformations to compute the QR
%            factorization of an N by P matrix.  Column pivoting is a
%            users option.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D5
%***TYPE      doubleprecision (SQRDC-S, DQRDC-D, CQRDC-C)
%***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX, ORTHOGONAL TRIANGULAR,
%             QR DECOMPOSITION
%***AUTHOR  Stewart, G. W., (U. of Maryland)
%***DESCRIPTION
%
%     DQRDC uses Householder transformations to compute the QR
%     factorization of an N by P matrix X.  Column pivoting
%     based on the 2-norms of the reduced columns may be
%     performed at the user's option.
%
%     On Entry
%
%        X       doubleprecision(LDX,P), where LDX .GE. N.
%                X contains the matrix whose decomposition is to be
%                computed.
%
%        LDX     INTEGER.
%                LDX is the leading dimension of the array X.
%
%        N       INTEGER.
%                N is the number of rows of the matrix X.
%
%        P       INTEGER.
%                P is the number of columns of the matrix X.
%
%        JPVT    INTEGER(P).
%                JPVT contains integers that control the selection
%                of the pivot columns.  The K-th column X(K) of X
%                is placed in one of three classes according to the
%                value of JPVT(K).
%
%                   If JPVT(K) .GT. 0, then X(K) is an initial
%                                      column.
%
%                   If JPVT(K) .EQ. 0, then X(K) is a free column.
%
%                   If JPVT(K) .LT. 0, then X(K) is a final column.
%
%                Before the decomposition is computed, initial columns
%                are moved to the beginning of the array X and final
%                columns to the end.  Both initial and final columns
%                are frozen in place during the computation and only
%                free columns are moved.  At the K-th stage of the
%                reduction, if X(K) is occupied by a free column
%                it is interchanged with the free column of largest
%                reduced norm.  JPVT is not referenced if
%                JOB .EQ. 0.
%
%        WORK    doubleprecision(P).
%                WORK is a work array.  WORK is not referenced if
%                JOB .EQ. 0.
%
%        JOB     INTEGER.
%                JOB is an integer that initiates column pivoting.
%                If JOB .EQ. 0, no pivoting is done.
%                If JOB .NE. 0, pivoting is done.
%
%     On Return
%
%        X       X contains in its upper triangle the upper
%                triangular matrix R of the QR factorization.
%                Below its diagonal X contains information from
%                which the orthogonal part of the decomposition
%                can be recovered.  Note that if pivoting has
%                been requested, the decomposition is not that
%                of the original matrix X but that of X
%                with its columns permuted as described by JPVT.
%
%        QRAUX   doubleprecision(P).
%                QRAUX contains further information required to recover
%                the orthogonal part of the decomposition.
%
%        JPVT    JPVT(K) contains the index of the column of the
%                original matrix that has been interchanged into
%                the K-th column, if pivoting was requested.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  DAXPY, DDOT, DNRM2, DSCAL, DSWAP
%***REVISION HISTORY  (YYMMDD)
%   780814  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DQRDC
jpvt_shape=size(jpvt);jpvt=reshape(jpvt,1,[]);
x_shape=size(x);x=reshape([x(:).',zeros(1,ceil(numel(x)./prod([ldx])).*prod([ldx])-numel(x))],ldx,[]);
qraux_shape=size(qraux);qraux=reshape(qraux,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
%
if isempty(j), j=0; end;
if isempty(jp), jp=0; end;
if isempty(l), l=0; end;
if isempty(lp1), lp1=0; end;
if isempty(lup), lup=0; end;
if isempty(maxj), maxj=0; end;
if isempty(pl), pl=0; end;
if isempty(pu), pu=0; end;
if isempty(maxnrm), maxnrm=0; end;
if isempty(tt), tt=0; end;
if isempty(nrmxl), nrmxl=0; end;
if isempty(t), t=0; end;
if isempty(negj), negj=false; end;
if isempty(swapj), swapj=false; end;
%
%***FIRST EXECUTABLE STATEMENT  DQRDC
pl = 1;
pu = 0;
if( job~=0 )
%
%        PIVOTING HAS BEEN REQUESTED.  REARRANGE THE COLUMNS
%        ACCORDING TO JPVT.
%
for j = 1 : p;
swapj = jpvt(j)>0;
negj = jpvt(j)<0;
jpvt(j) = fix(j);
if( negj )
jpvt(j) = fix(-j);
end;
if( swapj )
if( j~=pl )
[n,dumvar2,dumvar3,dumvar4]=dswap(n,x(sub2ind(size(x),1,pl):end),1,x(sub2ind(size(x),1,j):end),1);      x(sub2ind(size(x),1,pl):end)=dumvar2; x(sub2ind(size(x),1,j):end)=dumvar4; 
end;
jpvt(j) = fix(jpvt(pl));
jpvt(pl) = fix(j);
pl = fix(pl + 1);
end;
end; j = fix(p+1);
pu = fix(p);
for jj = 1 : p;
j = fix(p - jj + 1);
if( jpvt(j)<0 )
jpvt(j) = fix(-jpvt(j));
if( j~=pu )
[n,dumvar2,dumvar3,dumvar4]=dswap(n,x(sub2ind(size(x),1,pu):end),1,x(sub2ind(size(x),1,j):end),1);      x(sub2ind(size(x),1,pu):end)=dumvar2; x(sub2ind(size(x),1,j):end)=dumvar4; 
jp = fix(jpvt(pu));
jpvt(pu) = fix(jpvt(j));
jpvt(j) = fix(jp);
end;
pu = fix(pu - 1);
end;
end; jj = fix(p+1);
end;
%
%     COMPUTE THE NORMS OF THE FREE COLUMNS.
%
if( pu>=pl )
for j = pl : pu;
[qraux(j) ,n,x(sub2ind(size(x),1,j):end)]=dnrm2(n,x(sub2ind(size(x),1,j):end),1);
work(j) = qraux(j);
end; j = fix(pu+1);
end;
%
%     PERFORM THE HOUSEHOLDER REDUCTION OF X.
%
lup = fix(min(n,p));
for l = 1 : lup;
if( l>=pl && l<pu )
%
%           LOCATE THE COLUMN OF LARGEST NORM AND BRING IT
%           INTO THE PIVOT POSITION.
%
maxnrm = 0.0d0;
maxj = fix(l);
for j = l : pu;
if( qraux(j)>maxnrm )
maxnrm = qraux(j);
maxj = fix(j);
end;
end; j = fix(pu+1);
if( maxj~=l )
[n,dumvar2,dumvar3,dumvar4]=dswap(n,x(sub2ind(size(x),1,l):end),1,x(sub2ind(size(x),1,maxj):end),1);      x(sub2ind(size(x),1,l):end)=dumvar2; x(sub2ind(size(x),1,maxj):end)=dumvar4; 
qraux(maxj) = qraux(l);
work(maxj) = work(l);
jp = fix(jpvt(maxj));
jpvt(maxj) = fix(jpvt(l));
jpvt(l) = fix(jp);
end;
end;
qraux(l) = 0.0d0;
if( l~=n )
%
%           COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L.
%
[nrmxl ,dumvar2,x(sub2ind(size(x),l,l):end)]=dnrm2(n-l+1,x(sub2ind(size(x),l,l):end),1);
if( nrmxl~=0.0d0 )
if( x(l,l)~=0.0d0 )
nrmxl = (abs(nrmxl).*sign(x(l,l)));
end;
[dumvar1,dumvar2,x(sub2ind(size(x),l,l):end)]=dscal(n-l+1,1.0d0./nrmxl,x(sub2ind(size(x),l,l):end),1);
x(l,l) = 1.0d0 + x(l,l);
%
%              APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS,
%              UPDATING THE NORMS.
%
lp1 = fix(l + 1);
if( p>=lp1 )
for j = lp1 : p;
t = -ddot(n-l+1,x(sub2ind(size(x),l,l):end),1,x(sub2ind(size(x),l,j):end),1)./x(l,l);
[dumvar1,t,dumvar3,dumvar4,dumvar5]=daxpy(n-l+1,t,x(sub2ind(size(x),l,l):end),1,x(sub2ind(size(x),l,j):end),1);      x(sub2ind(size(x),l,l):end)=dumvar3; x(sub2ind(size(x),l,j):end)=dumvar5; 
if( j>=pl && j<=pu )
if( qraux(j)~=0.0d0 )
tt = 1.0d0 -(abs(x(l,j))./qraux(j)).^2;
tt = max(tt,0.0d0);
t = tt;
tt = 1.0d0 + 0.05d0.*tt.*(qraux(j)./work(j)).^2;
if( tt==1.0d0 )
[qraux(j) ,dumvar2,x(sub2ind(size(x),l+1,j):end)]=dnrm2(n-l,x(sub2ind(size(x),l+1,j):end),1);
work(j) = qraux(j);
else;
qraux(j) = qraux(j).*sqrt(t);
end;
end;
end;
end; j = fix(p+1);
end;
%
%              SAVE THE TRANSFORMATION.
%
qraux(l) = x(l,l);
x(l,l) = -nrmxl;
end;
end;
end; l = fix(lup+1);
jpvt_shape=zeros(jpvt_shape);jpvt_shape(:)=jpvt(1:numel(jpvt_shape));jpvt=jpvt_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
qraux_shape=zeros(qraux_shape);qraux_shape(:)=qraux(1:numel(qraux_shape));qraux=qraux_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK DQRFAC