caxcpy(n,a,x,incx,y,incy)
	 purpose:      constant times a conjugated vector plus a vector.
	 arguments:
	 n (in)        vector length
	 a (in)        complex factor
	 x (in)        added vector
	 incx (in)     x increments
	 y (io)        accumulator vector
	 incy (in)     y increments
	
cch1dn(n,R,ldr,u,rw,info)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine downdates R -> R1 so that
	               R1'*R1 = A - u*u'
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 u (io)        the vector determining the rank-1 update
	               on exit, u contains the reflector sines
	               used to transform R to R1.
	 rw (out)      cosine parts of reflectors.
	
	 info (out)    on exit, error code:
	                info = 0: success.
	                info = 1: update violates positive-definiteness.
	                info = 2: R is singular.
	
cch1up(n,R,ldr,u,w)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A + u*u'
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 u (io)        the vector determining the rank-1 update
	               on exit, u contains the rotation sines
	               used to transform R to R1.
	 w (out)       cosine parts of rotations.
	
cchdex(n,R,ldr,j,rw)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1].
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R.
	 R (io)        on entry, the original upper trapezoidal matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= n.
	 j (in)        the position of the deleted row/column.
	 rw (out)      a real workspace vector of size n.
	
cchinx(n,R,ldr,j,u,rw,info)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u,
	               jj = [1:j-1,j+1:n+1].
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R.
	 R (io)        on entry, the original upper trapezoidal matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= n+1.
	 j (in)        the position of the inserted row/column
	 u (io)        on entry, the inserted row/column.
	               on exit, u is destroyed.
	 rw (out)      real workspace vector of size n+1.
	 info (out)    on exit, error code:
	                info = 0: success.
	                info = 1: update violates positive-definiteness.
	                info = 2: R is singular.
	                info = 3: diagonal element of u is not real.
	
cchshx(n,R,ldr,i,j,w,rw)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A(p,p), where p is the permutation
	               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
	               [1:j-1,shift(j:i,+1),i+1:n] if j < i.
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 i (in)        the first index determining the range (see above).
	 j (in)        the second index determining the range (see above).
	 w (o)         a workspace vector of size n.
	 rw (o)        a real workspace vector of size n.
	
cgqvec(m,n,Q,ldq,u)
	 purpose:      given an unitary m-by-n matrix Q, n < m, generates
	               a vector u such that Q'*u = 0 and norm(u) = 1.
	 arguments:
	 m (in)        number of rows of matrix Q.
	 n (in)        number of columns of matrix Q.
	 Q (in)        the unitary matrix Q.
	 ldq (in)      leading dimension of Q.
	 u (out)       the generated vector.
	
clu1up(m,n,L,ldl,R,ldr,u,v)
	 purpose:      updates an LU factorization after rank-1 modification
	               i.e., given an m-by-k lower-triangular matrix L with unit
	               diagonal and a k-by-n upper-trapezoidal matrix R,
	               where k = min(m,n),
	               this subroutine updates L -> L1 and R -> R1 so that
	               L is again lower unit triangular, R upper trapezoidal,
	               and L1*R1 = L*R + u*v.'.
	               (complex version)
	 arguments:
	 m (in)        order of the matrix L.
	 n (in)        number of columns of the matrix U.
	 L (io)        on entry, the unit lower triangular matrix L.
	               on exit, the updated matrix L1.
	 ldl (in)      the leading dimension of L. ldl >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= min(m,n).
	 u (io)        the left m-vector. On exit, if k < m, u is destroyed.
	 v (io)        the right n-vector. On exit, v is destroyed.
	
	 REMARK:       Algorithm is due to
	               J. Bennett: Triangular factors of modified matrices,
	                           Numerische Mathematik, 7 (1965)
	
clup1up(m,n,L,ldl,R,ldr,p,u,v,w)
	 purpose:      updates a row-pivoted LU factorization after rank-1 modification
	               i.e., given an m-by-k lower-triangular matrix L with unit
	               diagonal, a k-by-n upper-trapezoidal matrix R, and a
	               permutation matrix P, where k = min(m,n),
	               this subroutine updates L -> L1, R -> R1 and P -> P1 so that
	               L is again lower unit triangular, R upper trapezoidal,
	               P permutation and P1'*L1*R1 = P'*L*R + u*v.'.
	               (complex version)
	 arguments:
	 m (in)        order of the matrix L.
	 n (in)        number of columns of the matrix U.
	 L (io)        on entry, the unit lower triangular matrix L.
	               on exit, the updated matrix L1.
	 ldl (in)      the leading dimension of L. ldl >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= min(m,n).
	 p (in)        the permutation vector representing P
	 u (in)        the left m-vector.
	 v (in)        the right n-vector.
	 w (work)      a workspace vector of size m.
	
	 REMARK:       Algorithm is due to
	               A. Kielbasinski, H. Schwetlick, Numerische Lineare
	               Algebra, Verlag Harri Deutsch, 1988
	
cqhqr(m,n,R,ldr,c,s)
	 purpose:      given an m-by-n upper Hessenberg matrix R, this
	               subroutine updates R to upper trapezoidal form
	               using min(m-1,n) Givens rotations.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix R
	 n (in)        number of columns of the matrix R
	 R (io)        on entry, the upper Hessenberg matrix R
	               on exit, the updated upper trapezoidal matrix
	 ldr (in)      leading dimension of R, >= m
	 c(out)        rotation cosines, size at least min(m-1,n)
	 s(out)        rotation sines, size at least min(m-1,n)
	
cqr1up(m,n,k,Q,ldq,R,ldr,u,v,w,rw)
	 purpose:      updates a QR factorization after rank-1 modification
	               i.e., given a m-by-k unitary Q and m-by-n upper
	               trapezoidal R, an m-vector u and n-vector v,
	               this subroutine updates Q -> Q1 and R -> R1 so that
	               Q1*R1 = Q*R + u*v', and Q1 is again unitary
	               and R1 upper trapezoidal.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n < m (economical form).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      the leading dimension of Q. ldq >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R..
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= k.
	 u (io)        the left m-vector. On exit, if k < m, u is destroyed.
	 v (io)        the right n-vector. On exit, v is destroyed.
	 w (out)       a workspace vector of size k.
	 rw (out)      a real workspace vector of size k.
	
cqrdec(m,n,k,Q,ldq,R,ldr,j,rw)
	 purpose:      updates a QR factorization after deleting
	               a column.
	               i.e., given an m-by-k unitary matrix Q, an k-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 remains unitary, R1 is upper
	               trapezoidal, and Q1*R1 = [A(:,1:j-1) A(:,j+1:n)],
	               where A = Q*R.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n < m (economical form,
	               basis dimension will decrease).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= k.
	 j (in)        the position of the deleted column in R.
	               1 <= j <= n.
	 rw (o)        a real workspace vector of size k-j.
	
cqrder(m,n,Q,ldq,R,ldr,j,w,rw)
	 purpose:      updates a QR factorization after deleting a row.
	               i.e., given an m-by-m unitary matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:m, this subroutine updates Q ->Q1 and an R -> R1
	               so that Q1 is again unitary, R1 upper trapezoidal,
	               and Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R.
	               (complex version)
	
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 Q (io)        on entry, the unitary matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= m.
	 j (in)        the position of the deleted row.
	 w (out)       a workspace vector of size m.
	 rw (out)      a real workspace vector of size m.
	
cqrinc(m,n,k,Q,ldq,R,ldr,j,x,rw)
	 purpose:      updates a QR factorization after inserting a new
	               column.
	               i.e., given an m-by-k unitary matrix Q, an m-by-n upper
	               trapezoidal matrix R and index j in the range 1:n+1,
	               this subroutine updates the matrix Q -> Q1 and R -> R1
	               so that Q1 is again unitary, R1 upper trapezoidal, and
	               Q1*R1 = [A(:,1:j-1); x; A(:,j:n)], where A = Q*R.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n <= m (economical form,
	               basis dimension will increase).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= min(m,n+1).
	 j (in)        the position of the new column in R1
	 x (in)        the column being inserted
	 rw (out)      a real workspace vector of size k.
	
cqrinr(m,n,Q,ldq,R,ldr,j,x,rw)
	 purpose:      updates a QR factorization after inserting a new
	               row.
	               i.e., given an m-by-m unitary matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:m+1, this subroutine updates Q -> Q1  and R -> R1
	               so that Q1 is again unitary, R1 upper trapezoidal,
	               and Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 Q (io)        on entry, the unitary matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m+1.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= m+1.
	 j (in)        the position of the new row in R1
	 x (io)        on entry, the row being added
	               on exit, x is destroyed.
	 rw (out)      a real workspace vector of size min(m,n).
	
cqrot(dir,m,n,Q,ldq,c,s)
	 purpose:      Apply a sequence of inv. rotations from right
	
	 arguments:
	 dir (in)      if 'B' or 'b', rotations are applied from backwards
	               if 'F' or 'f', from forwards.
	 m (in)        number of rows of matrix Q
	 n (in)        number of columns of the matrix Q
	 Q (io)        on entry, the matrix Q
	               on exit, the updated matrix Q1
	 ldq (in)      the leading dimension of Q
	 c (in)        n-1 rotation cosines
	 s (in)        n-1 rotation sines
	
cqrqh(m,n,R,ldr,c,s)
	 purpose:      brings an upper trapezoidal matrix R into upper
	               Hessenberg form using min(m-1,n) Givens rotations.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix R
	 n (in)        number of columns of the matrix R
	 R (io)        on entry, the upper Hessenberg matrix R
	               on exit, the updated upper trapezoidal matrix
	 ldr (in)      leading dimension of R, >= m
	 c(in)         rotation cosines, size at least min(m-1,n)
	 s(in)         rotation sines, size at least min(m-1,n)
	
cqrshc(m,n,k,Q,ldq,R,ldr,i,j,w,rw)
	 purpose:      updates a QR factorization after circular shift of
	               columns.
	               i.e., given an m-by-k unitary matrix Q, an k-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 is again unitary, R1 upper
	               trapezoidal, and
	               Q1*R1 = A(:,p), where A = Q*R and p is the permutation
	               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
	               [1:j-1,shift(j:i,+1),i+1:n] if j < i.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q1, and rows of R1. Must be
	               either k = m (full Q) or k = n <= m (economical form).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= k.
	 i (in)        the first index determining the range (see above)
	 j (in)        the second index determining the range (see above)
	 w (o)         a workspace vector of size k.
	 rw (o)        a real workspace vector of size k.
	
cqrtv1(n,u,w)
	 purpose:      generates a sequence of n-1 Givens rotations that
	               eliminate all but the first element of a vector u.
	 arguments:
	 n (in)        the length of the vector u
	 u (io)        on entry, the vector u.
	               on exit, u(2:n) contains the rotation sines, u(1)
	               contains the remaining element.
	 w (o)         on exit, w contains the rotation cosines.
	
dch1dn(n,R,ldr,u,w,info)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine downdates R -> R1 so that
	               R1'*R1 = A - u*u'
	               (real version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 u (io)        the vector determining the rank-1 update
	               on exit, u contains the reflector sines
	               used to transform R to R1.
	 w (out)       cosine parts of reflectors.
	
	 info (out)    on exit, error code:
	                info = 0: success.
	                info = 1: update violates positive-definiteness.
	                info = 2: R is singular.
	
dch1up(n,R,ldr,u,w)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A + u*u'
	               (real version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 u (io)        the vector determining the rank-1 update
	               on exit, u contains the rotation sines
	               used to transform R to R1.
	 w (out)       cosine parts of rotations.
	
dchdex(n,R,ldr,j,w)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1].
	               (real version)
	 arguments:
	 n (in)        the order of matrix R.
	 R (io)        on entry, the original upper trapezoidal matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= n.
	 j (in)        the position of the deleted row/column.
	 w (out)       a workspace vector of size n.
	
dchinx(n,R,ldr,j,u,w,info)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u,
	               jj = [1:j-1,j+1:n+1].
	               (real version)
	 arguments:
	 n (in)        the order of matrix R.
	 R (io)        on entry, the original upper trapezoidal matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= n+1.
	 j (in)        the position of the inserted row/column
	 u (io)        on entry, the inserted row/column.
	               on exit, u is destroyed.
	 w (out)       workspace vector of size n+1.
	 info (out)    on exit, error code:
	                info = 1: update violates positive-definiteness.
	                info = 2: R is singular.
	
dchshx(n,R,ldr,i,j,w)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A(p,p), where p is the permutation
	               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
	               [1:j-1,shift(j:i,+1),i+1:n] if j < i.
	               (real version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 i (in)        the first index determining the range (see above).
	 j (in)        the second index determining the range (see above).
	 w (o)         a workspace vector of size 2*n.
	
dgqvec(m,n,Q,ldq,u)
	 purpose:      given an orthogonal m-by-n matrix Q, n < m, generates
	               a vector u such that Q'*u = 0 and norm(u) = 1.
	 arguments:
	 m (in)        number of rows of matrix Q.
	 n (in)        number of columns of matrix Q.
	 Q (in)        the orthogonal matrix Q.
	 ldq (in)      leading dimension of Q.
	 u (out)       the generated vector.
	
dlu1up(m,n,L,ldl,R,ldr,u,v)
	 purpose:      updates an LU factorization after rank-1 modification
	               i.e., given an m-by-k lower-triangular matrix L with unit
	               diagonal and a k-by-n upper-trapezoidal matrix R,
	               where k = min(m,n),
	               this subroutine updates L -> L1 and R -> R1 so that
	               L is again lower unit triangular, R upper trapezoidal,
	               and L1*R1 = L*R + u*v.'.
	               (real version)
	 arguments:
	 m (in)        order of the matrix L.
	 n (in)        number of columns of the matrix U.
	 L (io)        on entry, the unit lower triangular matrix L.
	               on exit, the updated matrix L1.
	 ldl (in)      the leading dimension of L. ldl >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= min(m,n).
	 u (io)        the left m-vector. On exit, if k < m, u is destroyed.
	 v (io)        the right n-vector. On exit, v is destroyed.
	
	 REMARK:       Algorithm is due to
	               J. Bennett: Triangular factors of modified matrices,
	                           Numerische Mathematik, 7 (1965)
	
dlup1up(m,n,L,ldl,R,ldr,p,u,v,w)
	 purpose:      updates a row-pivoted LU factorization after rank-1 modification
	               i.e., given an m-by-k lower-triangular matrix L with unit
	               diagonal, a k-by-n upper-trapezoidal matrix R, and a
	               permutation matrix P, where k = min(m,n),
	               this subroutine updates L -> L1, R -> R1 and P -> P1 so that
	               L is again lower unit triangular, R upper trapezoidal,
	               P permutation and P1'*L1*R1 = P'*L*R + u*v.'.
	               (real version)
	 arguments:
	 m (in)        order of the matrix L.
	 n (in)        number of columns of the matrix U.
	 L (io)        on entry, the unit lower triangular matrix L.
	               on exit, the updated matrix L1.
	 ldl (in)      the leading dimension of L. ldl >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= min(m,n).
	 p (in)        the permutation vector representing P
	 u (in)        the left m-vector.
	 v (in)        the right n-vector.
	 w (work)      a workspace vector of size m.
	
	 REMARK:       Algorithm is due to
	               A. Kielbasinski, H. Schwetlick, Numerische Lineare
	               Algebra, Verlag Harri Deutsch, 1988
	
dqhqr(m,n,R,ldr,c,s)
	 purpose:      given an m-by-n upper Hessenberg matrix R, this
	               subroutine updates R to upper trapezoidal form
	               using min(m-1,n) Givens rotations.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix R
	 n (in)        number of columns of the matrix R
	 R (io)        on entry, the upper Hessenberg matrix R
	               on exit, the updated upper trapezoidal matrix
	 ldr (in)      leading dimension of R, >= m
	 c(out)        rotation cosines, size at least min(m-1,n)
	 s(out)        rotation sines, size at least min(m-1,n)
	
dqr1up(m,n,k,Q,ldq,R,ldr,u,v,w)
	 purpose:      updates a QR factorization after rank-1 modification
	               i.e., given a m-by-k orthogonal Q and m-by-n upper
	               trapezoidal R, an m-vector u and n-vector v,
	               this subroutine updates Q -> Q1 and R -> R1 so that
	               Q1*R1 = Q*R + u*v', and Q1 is again orthonormal
	               and R1 upper trapezoidal.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n < m (economical form).
	 Q (io)        on entry, the orthogonal m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      the leading dimension of Q. ldq >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R..
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= k.
	 u (io)        the left m-vector. On exit, if k < m, u is destroyed.
	 v (io)        the right n-vector. On exit, v is destroyed.
	 w (out)       a workspace vector of size 2*k
	
dqrdec(m,n,k,Q,ldq,R,ldr,j,w)
	 purpose:      updates a QR factorization after deleting
	               a column.
	               i.e., given an m-by-k orthogonal matrix Q, an k-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 remains orthogonal, R1 is upper
	               trapezoidal, and Q1*R1 = [A(:,1:j-1) A(:,j+1:n)],
	               where A = Q*R.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n < m (economical form,
	               basis dimension will decrease).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= k.
	 j (in)        the position of the deleted column in R.
	               1 <= j <= n.
	 w (o)         a workspace vector of size k-j.
	
dqrder(m,n,Q,ldq,R,ldr,j,w)
	 purpose:      updates a QR factorization after deleting a row.
	               i.e., given an m-by-m orthogonal matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:m, this subroutine updates Q ->Q1 and an R -> R1
	               so that Q1 is again orthogonal, R1 upper trapezoidal,
	               and Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R.
	               (real version)
	
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 Q (io)        on entry, the orthogonal matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= m.
	 j (in)        the position of the deleted row.
	 w (out)       a workspace vector of size 2*m.
	
dqrinc(m,n,k,Q,ldq,R,ldr,j,x,w)
	 purpose:      updates a QR factorization after inserting a new
	               column.
	               i.e., given an m-by-k orthogonal matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 is again orthogonal, R1 upper
	               trapezoidal, and Q1*R1 = [A(:,1:j-1); x; A(:,j:n)],
	               where A = Q*R.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n <= m (economical form,
	               basis dimension will increase).
	 Q (io)        on entry, the orthogonal m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= min(m,n+1).
	 j (in)        the position of the new column in R1
	 x (in)        the column being inserted
	 w (out)       a workspace vector of size k.
	
dqrinr(m,n,Q,ldq,R,ldr,j,x,w)
	 purpose:      updates a QR factorization after inserting a new
	               row.
	               i.e., given an m-by-m unitary matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:m+1, this subroutine updates Q -> Q1  and R -> R1
	               so that Q1 is again unitary, R1 upper trapezoidal,
	               and Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 Q (io)        on entry, the unitary matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m+1.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= m+1.
	 j (in)        the position of the new row in R1
	 x (io)        on entry, the row being added
	               on exit, x is destroyed.
	 w (out)       a workspace vector of size min(m,n).
	
dqrot(dir,m,n,Q,ldq,c,s)
	 purpose:      Apply a sequence of inv. rotations from right
	
	 arguments:
	 dir (in)      if 'B' or 'b', rotations are applied from backwards
	               if 'F' or 'f', from forwards.
	 m (in)        number of rows of matrix Q
	 n (in)        number of columns of the matrix Q
	 Q (io)        on entry, the matrix Q
	               on exit, the updated matrix Q1
	 ldq (in)      the leading dimension of Q
	 c (in)        n-1 rotation cosines
	 s (in)        n-1 rotation sines
	
dqrqh(m,n,R,ldr,c,s)
	 purpose:      brings an upper trapezoidal matrix R into upper
	               Hessenberg form using min(m-1,n) Givens rotations.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix R
	 n (in)        number of columns of the matrix R
	 R (io)        on entry, the upper Hessenberg matrix R
	               on exit, the updated upper trapezoidal matrix
	 ldr (in)      leading dimension of R, >= m
	 c(in)         rotation cosines, size at least min(m-1,n)
	 s(in)         rotation sines, size at least min(m-1,n)
	
dqrshc(m,n,k,Q,ldq,R,ldr,i,j,w)
	 purpose:      updates a QR factorization after circular shift of
	               columns.
	               i.e., given an m-by-k orthogonal matrix Q, an k-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 is again orthogonal, R1 upper
	               trapezoidal, and
	               Q1*R1 = A(:,p), where A = Q*R and p is the permutation
	               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
	               [1:j-1,shift(j:i,+1),i+1:n] if j < i.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q1, and rows of R1. Must be
	               either k = m (full Q) or k = n <= m (economical form).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= k.
	 i (in)        the first index determining the range (see above)
	 j (in)        the second index determining the range (see above)
	 w (o)         a workspace vector of size 2*k.
	
dqrtv1(n,u,w)
	 purpose:      generates a sequence of n-1 Givens rotations that
	               eliminate all but the first element of a vector u.
	 arguments:
	 n (in)        the length of the vector u
	 u (io)        on entry, the vector u.
	               on exit, u(2:n) contains the rotation sines, u(1)
	               contains the remaining element.
	 w (o)         on exit, w contains the rotation cosines.
	
sch1dn(n,R,ldr,u,w,info)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine downdates R -> R1 so that
	               R1'*R1 = A - u*u'
	               (real version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 u (io)        the vector determining the rank-1 update
	               on exit, u contains the reflector sines
	               used to transform R to R1.
	 w (out)       cosine parts of reflectors.
	
	 info (out)    on exit, error code:
	                info = 0: success.
	                info = 1: update violates positive-definiteness.
	                info = 2: R is singular.
	
sch1up(n,R,ldr,u,w)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A + u*u'
	               (real version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 u (io)        the vector determining the rank-1 update
	               on exit, u contains the rotation sines
	               used to transform R to R1.
	 w (out)       cosine parts of rotations.
	
schdex(n,R,ldr,j,w)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1].
	               (real version)
	 arguments:
	 n (in)        the order of matrix R.
	 R (io)        on entry, the original upper trapezoidal matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= n.
	 j (in)        the position of the deleted row/column.
	 w (out)       a workspace vector of size n.
	
schinx(n,R,ldr,j,u,w,info)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u,
	               jj = [1:j-1,j+1:n+1].
	               (real version)
	 arguments:
	 n (in)        the order of matrix R.
	 R (io)        on entry, the original upper trapezoidal matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= n+1.
	 j (in)        the position of the inserted row/column
	 u (io)        on entry, the inserted row/column.
	               on exit, u is destroyed.
	 w (out)       workspace vector of size n+1.
	 info (out)    on exit, error code:
	                info = 1: update violates positive-definiteness.
	                info = 2: R is singular.
	
schshx(n,R,ldr,i,j,w)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A(p,p), where p is the permutation
	               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
	               [1:j-1,shift(j:i,+1),i+1:n] if j < i.
	               (real version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 i (in)        the first index determining the range (see above).
	 j (in)        the second index determining the range (see above).
	 w (o)         a workspace vector of size 2*n.
	
sgqvec(m,n,Q,ldq,u)
	 purpose:      given an orthogonal m-by-n matrix Q, n < m, generates
	               a vector u such that Q'*u = 0 and norm(u) = 1.
	 arguments:
	 m (in)        number of rows of matrix Q.
	 n (in)        number of columns of matrix Q.
	 Q (in)        the orthogonal matrix Q.
	 ldq (in)      leading dimension of Q.
	 u (out)       the generated vector.
	
slu1up(m,n,L,ldl,R,ldr,u,v)
	 purpose:      updates an LU factorization after rank-1 modification
	               i.e., given an m-by-k lower-triangular matrix L with unit
	               diagonal and a k-by-n upper-trapezoidal matrix R,
	               where k = min(m,n),
	               this subroutine updates L -> L1 and R -> R1 so that
	               L is again lower unit triangular, R upper trapezoidal,
	               and L1*R1 = L*R + u*v.'.
	               (real version)
	 arguments:
	 m (in)        order of the matrix L.
	 n (in)        number of columns of the matrix U.
	 L (io)        on entry, the unit lower triangular matrix L.
	               on exit, the updated matrix L1.
	 ldl (in)      the leading dimension of L. ldl >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= min(m,n).
	 u (io)        the left m-vector. On exit, if k < m, u is destroyed.
	 v (io)        the right n-vector. On exit, v is destroyed.
	
	 REMARK:       Algorithm is due to
	               J. Bennett: Triangular factors of modified matrices,
	                           Numerische Mathematik, 7 (1965)
	
slup1up(m,n,L,ldl,R,ldr,p,u,v,w)
	 purpose:      updates a row-pivoted LU factorization after rank-1 modification
	               i.e., given an m-by-k lower-triangular matrix L with unit
	               diagonal, a k-by-n upper-trapezoidal matrix R, and a
	               permutation matrix P, where k = min(m,n),
	               this subroutine updates L -> L1, R -> R1 and P -> P1 so that
	               L is again lower unit triangular, R upper trapezoidal,
	               P permutation and P1'*L1*R1 = P'*L*R + u*v.'.
	               (real version)
	 arguments:
	 m (in)        order of the matrix L.
	 n (in)        number of columns of the matrix U.
	 L (io)        on entry, the unit lower triangular matrix L.
	               on exit, the updated matrix L1.
	 ldl (in)      the leading dimension of L. ldl >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= min(m,n).
	 p (in)        the permutation vector representing P
	 u (in)        the left m-vector.
	 v (in)        the right n-vector.
	 w (work)      a workspace vector of size m.
	
	 REMARK:       Algorithm is due to
	               A. Kielbasinski, H. Schwetlick, Numerische Lineare
	               Algebra, Verlag Harri Deutsch, 1988
	
sqhqr(m,n,R,ldr,c,s)
	 purpose:      given an m-by-n upper Hessenberg matrix R, this
	               subroutine updates R to upper trapezoidal form
	               using min(m-1,n) Givens rotations.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix R
	 n (in)        number of columns of the matrix R
	 R (io)        on entry, the upper Hessenberg matrix R
	               on exit, the updated upper trapezoidal matrix
	 ldr (in)      leading dimension of R, >= m
	 c(out)        rotation cosines, size at least min(m-1,n)
	 s(out)        rotation sines, size at least min(m-1,n)
	
sqr1up(m,n,k,Q,ldq,R,ldr,u,v,w)
	 purpose:      updates a QR factorization after rank-1 modification
	               i.e., given a m-by-k orthogonal Q and m-by-n upper
	               trapezoidal R, an m-vector u and n-vector v,
	               this subroutine updates Q -> Q1 and R -> R1 so that
	               Q1*R1 = Q*R + u*v', and Q1 is again orthonormal
	               and R1 upper trapezoidal.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n < m (economical form).
	 Q (io)        on entry, the orthogonal m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      the leading dimension of Q. ldq >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R..
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= k.
	 u (io)        the left m-vector. On exit, if k < m, u is destroyed.
	 v (io)        the right n-vector. On exit, v is destroyed.
	 w (out)       a workspace vector of size 2*k
	
sqrdec(m,n,k,Q,ldq,R,ldr,j,w)
	 purpose:      updates a QR factorization after deleting
	               a column.
	               i.e., given an m-by-k orthogonal matrix Q, an k-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 remains orthogonal, R1 is upper
	               trapezoidal, and Q1*R1 = [A(:,1:j-1) A(:,j+1:n)],
	               where A = Q*R.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n < m (economical form,
	               basis dimension will decrease).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= k.
	 j (in)        the position of the deleted column in R.
	               1 <= j <= n.
	 w (o)         a workspace vector of size k-j.
	
sqrder(m,n,Q,ldq,R,ldr,j,w)
	 purpose:      updates a QR factorization after deleting a row.
	               i.e., given an m-by-m orthogonal matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:m, this subroutine updates Q ->Q1 and an R -> R1
	               so that Q1 is again orthogonal, R1 upper trapezoidal,
	               and Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R.
	               (real version)
	
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 Q (io)        on entry, the orthogonal matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= m.
	 j (in)        the position of the deleted row.
	 w (out)       a workspace vector of size 2*m.
	
sqrinc(m,n,k,Q,ldq,R,ldr,j,x,w)
	 purpose:      updates a QR factorization after inserting a new
	               column.
	               i.e., given an m-by-k orthogonal matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 is again orthogonal, R1 upper
	               trapezoidal, and Q1*R1 = [A(:,1:j-1); x; A(:,j:n)],
	               where A = Q*R.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n <= m (economical form,
	               basis dimension will increase).
	 Q (io)        on entry, the orthogonal m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= min(m,n+1).
	 j (in)        the position of the new column in R1
	 x (in)        the column being inserted
	 w (out)       a workspace vector of size k.
	
sqrinr(m,n,Q,ldq,R,ldr,j,x,w)
	 purpose:      updates a QR factorization after inserting a new
	               row.
	               i.e., given an m-by-m unitary matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:m+1, this subroutine updates Q -> Q1  and R -> R1
	               so that Q1 is again unitary, R1 upper trapezoidal,
	               and Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 Q (io)        on entry, the unitary matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m+1.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= m+1.
	 j (in)        the position of the new row in R1
	 x (io)        on entry, the row being added
	               on exit, x is destroyed.
	 w (out)       a workspace vector of size min(m,n).
	
sqrot(dir,m,n,Q,ldq,c,s)
	 purpose:      Apply a sequence of inv. rotations from right
	
	 arguments:
	 dir (in)      if 'B' or 'b', rotations are applied from backwards
	               if 'F' or 'f', from forwards.
	 m (in)        number of rows of matrix Q
	 n (in)        number of columns of the matrix Q
	 Q (io)        on entry, the matrix Q
	               on exit, the updated matrix Q1
	 ldq (in)      the leading dimension of Q
	 c (in)        n-1 rotation cosines
	 s (in)        n-1 rotation sines
	
sqrqh(m,n,R,ldr,c,s)
	 purpose:      brings an upper trapezoidal matrix R into upper
	               Hessenberg form using min(m-1,n) Givens rotations.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix R
	 n (in)        number of columns of the matrix R
	 R (io)        on entry, the upper Hessenberg matrix R
	               on exit, the updated upper trapezoidal matrix
	 ldr (in)      leading dimension of R, >= m
	 c(in)         rotation cosines, size at least min(m-1,n)
	 s(in)         rotation sines, size at least min(m-1,n)
	
sqrshc(m,n,k,Q,ldq,R,ldr,i,j,w)
	 purpose:      updates a QR factorization after circular shift of
	               columns.
	               i.e., given an m-by-k orthogonal matrix Q, an k-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 is again orthogonal, R1 upper
	               trapezoidal, and
	               Q1*R1 = A(:,p), where A = Q*R and p is the permutation
	               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
	               [1:j-1,shift(j:i,+1),i+1:n] if j < i.
	               (real version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q1, and rows of R1. Must be
	               either k = m (full Q) or k = n <= m (economical form).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= k.
	 i (in)        the first index determining the range (see above)
	 j (in)        the second index determining the range (see above)
	 w (o)         a workspace vector of size 2*k.
	
sqrtv1(n,u,w)
	 purpose:      generates a sequence of n-1 Givens rotations that
	               eliminate all but the first element of a vector u.
	 arguments:
	 n (in)        the length of the vector u
	 u (io)        on entry, the vector u.
	               on exit, u(2:n) contains the rotation sines, u(1)
	               contains the remaining element.
	 w (o)         on exit, w contains the rotation cosines.
	
zaxcpy(n,a,x,incx,y,incy)
	 purpose:      constant times a conjugated vector plus a vector.
	 arguments:
	 n (in)        vector length
	 a (in)        complex factor
	 x (in)        added vector
	 incx (in)     x increments
	 y (io)        accumulator vector
	 incy (in)     y increments
	
zch1dn(n,R,ldr,u,rw,info)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine downdates R -> R1 so that
	               R1'*R1 = A - u*u'
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 u (io)        the vector determining the rank-1 update
	               on exit, u contains the reflector sines
	               used to transform R to R1.
	 rw (out)      cosine parts of reflectors.
	
	 info (out)    on exit, error code:
	                info = 0: success.
	                info = 1: update violates positive-definiteness.
	                info = 2: R is singular.
	
zch1up(n,R,ldr,u,w)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A + u*u'
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 u (io)        the vector determining the rank-1 update
	               on exit, u contains the rotation sines
	               used to transform R to R1.
	 w (out)       cosine parts of rotations.
	
zchdex(n,R,ldr,j,rw)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A(jj,jj), where jj = [1:j-1,j+1:n+1].
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R.
	 R (io)        on entry, the original upper trapezoidal matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= n.
	 j (in)        the position of the deleted row/column.
	 rw (out)      a real workspace vector of size n.
	
zchinx(n,R,ldr,j,u,rw,info)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a hermitian positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A1, A1(jj,jj) = A, A(j,:) = u', A(:,j) = u,
	               jj = [1:j-1,j+1:n+1].
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R.
	 R (io)        on entry, the original upper trapezoidal matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= n+1.
	 j (in)        the position of the inserted row/column
	 u (io)        on entry, the inserted row/column.
	               on exit, u is destroyed.
	 rw (out)      real workspace vector of size n+1.
	 info (out)    on exit, error code:
	                info = 1: update violates positive-definiteness.
	                info = 2: R is singular.
	                info = 3: diagonal element of u is not real.
	
zchshx(n,R,ldr,i,j,w,rw)
	 purpose:      given an upper triangular matrix R that is a Cholesky
	               factor of a symmetric positive definite matrix A, i.e.
	               A = R'*R, this subroutine updates R -> R1 so that
	               R1'*R1 = A(p,p), where p is the permutation
	               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
	               [1:j-1,shift(j:i,+1),i+1:n] if j < i.
	               (complex version)
	 arguments:
	 n (in)        the order of matrix R
	 R (io)        on entry, the upper triangular matrix R
	               on exit, the updated matrix R1
	 ldr (in)      leading dimension of R. ldr >= n.
	 i (in)        the first index determining the range (see above).
	 j (in)        the second index determining the range (see above).
	 w (o)         a workspace vector of size n.
	 rw (o)        a real workspace vector of size n.
	
zgqvec(m,n,Q,ldq,u)
	 purpose:      given an unitary m-by-n matrix Q, n < m, generates
	               a vector u such that Q'*u = 0 and norm(u) = 1.
	 arguments:
	 m (in)        number of rows of matrix Q.
	 n (in)        number of columns of matrix Q.
	 Q (in)        the unitary matrix Q.
	 ldq (in)      leading dimension of Q.
	 u (out)       the generated vector.
	
zlu1up(m,n,L,ldl,R,ldr,u,v)
	 purpose:      updates an LU factorization after rank-1 modification
	               i.e., given an m-by-k lower-triangular matrix L with unit
	               diagonal and a k-by-n upper-trapezoidal matrix R,
	               where k = min(m,n),
	               this subroutine updates L -> L1 and R -> R1 so that
	               L is again lower unit triangular, R upper trapezoidal,
	               and L1*R1 = L*R + u*v.'.
	               (complex version)
	 arguments:
	 m (in)        order of the matrix L.
	 n (in)        number of columns of the matrix U.
	 L (io)        on entry, the unit lower triangular matrix L.
	               on exit, the updated matrix L1.
	 ldl (in)      the leading dimension of L. ldl >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= min(m,n).
	 u (io)        the left m-vector. On exit, if k < m, u is destroyed.
	 v (io)        the right n-vector. On exit, v is destroyed.
	
	 REMARK:       Algorithm is due to
	               J. Bennett: Triangular factors of modified matrices,
	                           Numerische Mathematik, 7 (1965)
	
zlup1up(m,n,L,ldl,R,ldr,p,u,v,w)
	 purpose:      updates a row-pivoted LU factorization after rank-1 modification
	               i.e., given an m-by-k lower-triangular matrix L with unit
	               diagonal, a k-by-n upper-trapezoidal matrix R, and a
	               permutation matrix P, where k = min(m,n),
	               this subroutine updates L -> L1, R -> R1 and P -> P1 so that
	               L is again lower unit triangular, R upper trapezoidal,
	               P permutation and P1'*L1*R1 = P'*L*R + u*v.'.
	               (real version)
	 arguments:
	 m (in)        order of the matrix L.
	 n (in)        number of columns of the matrix U.
	 L (io)        on entry, the unit lower triangular matrix L.
	               on exit, the updated matrix L1.
	 ldl (in)      the leading dimension of L. ldl >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= min(m,n).
	 p (in)        the permutation vector representing P
	 u (in)        the left m-vector.
	 v (in)        the right n-vector.
	 w (work)      a workspace vector of size m.
	
	 REMARK:       Algorithm is due to
	               A. Kielbasinski, H. Schwetlick, Numerische Lineare
	               Algebra, Verlag Harri Deutsch, 1988
	
zqhqr(m,n,R,ldr,c,s)
	 purpose:      given an m-by-n upper Hessenberg matrix R, this
	               subroutine updates R to upper trapezoidal form
	               using min(m-1,n) Givens rotations.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix R
	 n (in)        number of columns of the matrix R
	 R (io)        on entry, the upper Hessenberg matrix R
	               on exit, the updated upper trapezoidal matrix
	 ldr (in)      leading dimension of R, >= m
	 c(out)        rotation cosines, size at least min(m-1,n)
	 s(out)        rotation sines, size at least min(m-1,n)
	
zqr1up(m,n,k,Q,ldq,R,ldr,u,v,w,rw)
	 purpose:      updates a QR factorization after rank-1 modification
	               i.e., given a m-by-k unitary Q and m-by-n upper
	               trapezoidal R, an m-vector u and n-vector v,
	               this subroutine updates Q -> Q1 and R -> R1 so that
	               Q1*R1 = Q*R + u*v', and Q1 is again unitary
	               and R1 upper trapezoidal.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n < m (economical form).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      the leading dimension of Q. ldq >= m.
	 R (io)        on entry, the upper trapezoidal m-by-n matrix R..
	               on exit, the updated matrix R1.
	 ldr (in)      the leading dimension of R. ldr >= k.
	 u (io)        the left m-vector. On exit, if k < m, u is destroyed.
	 v (io)        the right n-vector. On exit, v is destroyed.
	 w (out)       a workspace vector of size k.
	 rw (out)      a real workspace vector of size k.
	
zqrdec(m,n,k,Q,ldq,R,ldr,j,rw)
	 purpose:      updates a QR factorization after deleting
	               a column.
	               i.e., given an m-by-k unitary matrix Q, an k-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 remains unitary, R1 is upper
	               trapezoidal, and Q1*R1 = [A(:,1:j-1) A(:,j+1:n)],
	               where A = Q*R.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n < m (economical form,
	               basis dimension will decrease).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= k.
	 j (in)        the position of the deleted column in R.
	               1 <= j <= n.
	 rw (o)        a real workspace vector of size k-j.
	
zqrder(m,n,Q,ldq,R,ldr,j,w,rw)
	 purpose:      updates a QR factorization after deleting a row.
	               i.e., given an m-by-m unitary matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:m, this subroutine updates Q ->Q1 and an R -> R1
	               so that Q1 is again unitary, R1 upper trapezoidal,
	               and Q1*R1 = [A(1:j-1,:); A(j+1:m,:)], where A = Q*R.
	               (complex version)
	
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 Q (io)        on entry, the unitary matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= m.
	 j (in)        the position of the deleted row.
	 w (out)       a workspace vector of size m.
	 rw (out)      a real workspace vector of size m.
	
zqrinc(m,n,k,Q,ldq,R,ldr,j,x,rw)
	 purpose:      updates a QR factorization after inserting a new
	               column.
	               i.e., given an m-by-k unitary matrix Q, an m-by-n upper
	               trapezoidal matrix R and index j in the range 1:n+1,
	               this subroutine updates the matrix Q -> Q1 and R -> R1
	               so that Q1 is again unitary, R1 upper trapezoidal, and
	               Q1*R1 = [A(:,1:j-1); x; A(:,j:n)], where A = Q*R.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q, and rows of R. Must be
	               either k = m (full Q) or k = n <= m (economical form,
	               basis dimension will increase).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= min(m,n+1).
	 j (in)        the position of the new column in R1
	 x (in)        the column being inserted
	 rw (out)      a real workspace vector of size k.
	
zqrinr(m,n,Q,ldq,R,ldr,j,x,rw)
	 purpose:      updates a QR factorization after inserting a new
	               row.
	               i.e., given an m-by-m unitary matrix Q, an m-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:m+1, this subroutine updates Q -> Q1  and R -> R1
	               so that Q1 is again unitary, R1 upper trapezoidal,
	               and Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 Q (io)        on entry, the unitary matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m+1.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= m+1.
	 j (in)        the position of the new row in R1
	 x (io)        on entry, the row being added
	               on exit, x is destroyed.
	 rw (out)      a real workspace vector of size min(m,n).
	
zqrot(dir,m,n,Q,ldq,c,s)
	 purpose:      Apply a sequence of inv. rotations from right
	
	 arguments:
	 dir (in)      if 'B' or 'b', rotations are applied from backwards
	               if 'F' or 'f', from forwards.
	 m (in)        number of rows of matrix Q
	 n (in)        number of columns of the matrix Q
	 Q (io)        on entry, the matrix Q
	               on exit, the updated matrix Q1
	 ldq (in)      the leading dimension of Q
	 c (in)        n-1 rotation cosines
	 s (in)        n-1 rotation sines
	
zqrqh(m,n,R,ldr,c,s)
	 purpose:      brings an upper trapezoidal matrix R into upper
	               Hessenberg form using min(m-1,n) Givens rotations.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix R
	 n (in)        number of columns of the matrix R
	 R (io)        on entry, the upper Hessenberg matrix R
	               on exit, the updated upper trapezoidal matrix
	 ldr (in)      leading dimension of R, >= m
	 c(in)         rotation cosines, size at least min(m-1,n)
	 s(in)         rotation sines, size at least min(m-1,n)
	
zqrshc(m,n,k,Q,ldq,R,ldr,i,j,w,rw)
	 purpose:      updates a QR factorization after circular shift of
	               columns.
	               i.e., given an m-by-k unitary matrix Q, an k-by-n
	               upper trapezoidal matrix R and index j in the range
	               1:n+1, this subroutine updates the matrix Q -> Q1 and
	               R -> R1 so that Q1 is again unitary, R1 upper
	               trapezoidal, and
	               Q1*R1 = A(:,p), where A = Q*R and p is the permutation
	               [1:i-1,shift(i:j,-1),j+1:n] if i < j  or
	               [1:j-1,shift(j:i,+1),i+1:n] if j < i.
	               (complex version)
	 arguments:
	 m (in)        number of rows of the matrix Q.
	 n (in)        number of columns of the matrix R.
	 k (in)        number of columns of Q1, and rows of R1. Must be
	               either k = m (full Q) or k = n <= m (economical form).
	 Q (io)        on entry, the unitary m-by-k matrix Q.
	               on exit, the updated matrix Q1.
	 ldq (in)      leading dimension of Q. ldq >= m.
	 R (io)        on entry, the original matrix R.
	               on exit, the updated matrix R1.
	 ldr (in)      leading dimension of R. ldr >= k.
	 i (in)        the first index determining the range (see above)
	 j (in)        the second index determining the range (see above)
	 w (o)         a workspace vector of size k.
	 rw (o)        a real workspace vector of size k.
	
zqrtv1(n,u,w)
	 purpose:      generates a sequence of n-1 Givens rotations that
	               eliminate all but the first element of a vector u.
	 arguments:
	 n (in)        the length of the vector u
	 u (io)        on entry, the vector u.
	               on exit, u(2:n) contains the rotation sines, u(1)
	               contains the remaining element.
	 w (o)         on exit, w contains the rotation cosines.
	
