[r,ldr,p,k,l,z,ldz,nz,c,s,job]=schex(r,ldr,p,k,l,z,ldz,nz,c,s,job);

function [r,ldr,p,k,l,z,ldz,nz,c,s,job]=schex(r,ldr,p,k,l,z,ldz,nz,c,s,job);
%***BEGIN PROLOGUE  SCHEX
%***PURPOSE  Update the Cholesky factorization  A=TRANS(R)*R  of A
%            positive definite matrix A of order P under diagonal
%            permutations of the form TRANS(E)*A*E, where E is a
%            permutation matrix.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D7B
%***TYPE      SINGLE PRECISION (SCHEX-S, DCHEX-D, CCHEX-C)
%***KEYWORDS  CHOLESKY DECOMPOSITION, EXCHANGE, LINEAR ALGEBRA, LINPACK,
%             MATRIX, POSITIVE DEFINITE
%***AUTHOR  Stewart, G. W., (U. of Maryland)
%***DESCRIPTION
%
%     SCHEX updates the Cholesky factorization
%
%                   A = TRANS(R)*R
%
%     of a positive definite matrix A of order P under diagonal
%     permutations of the form
%
%                   TRANS(E)*A*E
%
%     where E is a permutation matrix.  Specifically, given
%     an upper triangular matrix R and a permutation matrix
%     E (which is specified by K, L, and JOB), SCHEX determines
%     an orthogonal matrix U such that
%
%                           U*R*E = RR,
%
%     where RR is upper triangular.  At the users option, the
%     transformation U will be multiplied into the array Z.
%     If A = TRANS(X)*X, so that R is the triangular part of the
%     QR factorization of X, then RR is the triangular part of the
%     QR factorization of X*E, i.e., X with its columns permuted.
%     For a less terse description of what SCHEX does and how
%     it may be applied, see the LINPACK guide.
%
%     The matrix Q is determined as the product U(L-K)*...*U(1)
%     of plane rotations of the form
%
%                           (    C(I)       S(I) )
%                           (                    ) ,
%                           (    -S(I)      C(I) )
%
%     where C(I) is real.  The rows these rotations operate on
%     are described below.
%
%     There are two types of permutations, which are determined
%     by the value of JOB.
%
%     1. Right circular shift (JOB = 1).
%
%         The columns are rearranged in the following order.
%
%                1,...,K-1,L,K,K+1,...,L-1,L+1,...,P.
%
%         U is the product of L-K rotations U(I), where U(I)
%         acts in the (L-I,L-I+1)-plane.
%
%     2. Left circular shift (JOB = 2).
%         The columns are rearranged in the following order
%
%                1,...,K-1,K+1,K+2,...,L,K,L+1,...,P.
%
%         U is the product of L-K rotations U(I), where U(I)
%         acts in the (K+I-1,K+I)-plane.
%
%     On Entry
%
%         R      REAL(LDR,P), where LDR .GE. P.
%                R contains the upper triangular factor
%                that is to be updated.  Elements of R
%                below the diagonal are not referenced.
%
%         LDR    INTEGER.
%                LDR is the leading dimension of the array R.
%
%         P      INTEGER.
%                P is the order of the matrix R.
%
%         K      INTEGER.
%                K is the first column to be permuted.
%
%         L      INTEGER.
%                L is the last column to be permuted.
%                L must be strictly greater than K.
%
%         Z      REAL(LDZ,NZ), where LDZ.GE.P.
%                Z is an array of NZ P-vectors into which the
%                transformation U is multiplied.  Z is
%                not referenced if NZ = 0.
%
%         LDZ    INTEGER.
%                LDZ is the leading dimension of the array Z.
%
%         NZ     INTEGER.
%                NZ is the number of columns of the matrix Z.
%
%         JOB    INTEGER.
%                JOB determines the type of permutation.
%                       JOB = 1  right circular shift.
%                       JOB = 2  left circular shift.
%
%     On Return
%
%         R      contains the updated factor.
%
%         Z      contains the updated matrix Z.
%
%         C      REAL(P).
%                C contains the cosines of the transforming rotations.
%
%         S      REAL(P).
%                S contains the sines of the transforming rotations.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  SROTG
%***REVISION HISTORY  (YYMMDD)
%   780814  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  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  SCHEX
persistent i ii il iu j jj km1 kp1 lm1 lmk t ; 

r_shape=size(r);r=reshape([r(:).',zeros(1,ceil(numel(r)./prod([ldr])).*prod([ldr])-numel(r))],ldr,[]);
z_shape=size(z);z=reshape([z(:).',zeros(1,ceil(numel(z)./prod([ldz])).*prod([ldz])-numel(z))],ldz,[]);
s_shape=size(s);s=reshape(s,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
%
if isempty(i), i=0; end;
if isempty(ii), ii=0; end;
if isempty(il), il=0; end;
if isempty(iu), iu=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(km1), km1=0; end;
if isempty(kp1), kp1=0; end;
if isempty(lmk), lmk=0; end;
if isempty(lm1), lm1=0; end;
if isempty(t), t=0; end;
%
%     INITIALIZE
%
%***FIRST EXECUTABLE STATEMENT  SCHEX
km1 = fix(k - 1);
kp1 = fix(k + 1);
lmk = fix(l - k);
lm1 = fix(l - 1);
%
%     PERFORM THE APPROPRIATE TASK.
%
if( job==2 )
%
%     LEFT CIRCULAR SHIFT
%
%
%        REORDER THE COLUMNS
%
for i = 1 : k;
ii = fix(lmk + i);
s(ii) = r(i,k);
end; i = fix(k+1);
for j = k : lm1;
for i = 1 : j;
r(i,j) = r(i,j+1);
end; i = fix(j+1);
jj = fix(j - km1);
s(jj) = r(j+1,j+1);
end; j = fix(lm1+1);
for i = 1 : k;
ii = fix(lmk + i);
r(i,l) = s(ii);
end; i = fix(k+1);
for i = kp1 : l;
r(i,l) = 0.0e0;
end; i = fix(l+1);
%
%        REDUCTION LOOP.
%
for j = k : p;
if( j~=k )
%
%              APPLY THE ROTATIONS.
%
iu = fix(min(j-1,l-1));
for i = k : iu;
ii = fix(i - k + 1);
t = c(ii).*r(i,j) + s(ii).*r(i+1,j);
r(i+1,j) = c(ii).*r(i+1,j) - s(ii).*r(i,j);
r(i,j) = t;
end; i = fix(iu+1);
end;
if( j<l )
jj = fix(j - k + 1);
t = s(jj);
[r(j,j),t,c(jj),s(jj)]=srotg(r(j,j),t,c(jj),s(jj));
end;
end; j = fix(p+1);
%
%        APPLY THE ROTATIONS TO Z.
%
if( nz>=1 )
for j = 1 : nz;
for i = k : lm1;
ii = fix(i - km1);
t = c(ii).*z(i,j) + s(ii).*z(i+1,j);
z(i+1,j) = c(ii).*z(i+1,j) - s(ii).*z(i,j);
z(i,j) = t;
end; i = fix(lm1+1);
end; j = fix(nz+1);
end;
else;
%
%     RIGHT CIRCULAR SHIFT.
%
%
%        REORDER THE COLUMNS.
%
for i = 1 : l;
ii = fix(l - i + 1);
s(i) = r(ii,l);
end; i = fix(l+1);
for jj = k : lm1;
j = fix(lm1 - jj + k);
for i = 1 : j;
r(i,j+1) = r(i,j);
end; i = fix(j+1);
r(j+1,j+1) = 0.0e0;
end; jj = fix(lm1+1);
if( k~=1 )
for i = 1 : km1;
ii = fix(l - i + 1);
r(i,k) = s(ii);
end; i = fix(km1+1);
end;
%
%        CALCULATE THE ROTATIONS.
%
t = s(1);
for i = 1 : lmk;
[s(i+1),t,c(i),s(i)]=srotg(s(i+1),t,c(i),s(i));
t = s(i+1);
end; i = fix(lmk+1);
r(k,k) = t;
for j = kp1 : p;
il = fix(max(1,l-j+1));
for ii = il : lmk;
i = fix(l - ii);
t = c(ii).*r(i,j) + s(ii).*r(i+1,j);
r(i+1,j) = c(ii).*r(i+1,j) - s(ii).*r(i,j);
r(i,j) = t;
end; ii = fix(lmk+1);
end; j = fix(p+1);
%
%        IF REQUIRED, APPLY THE TRANSFORMATIONS TO Z.
%
if( nz>=1 )
for j = 1 : nz;
for ii = 1 : lmk;
i = fix(l - ii);
t = c(ii).*z(i,j) + s(ii).*z(i+1,j);
z(i+1,j) = c(ii).*z(i+1,j) - s(ii).*z(i,j);
z(i,j) = t;
end; ii = fix(lmk+1);
end; j = fix(nz+1);
end;
end;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
end
%DECK SCHKW