Revision 5ff90986884ba22a58fb60d64341ab718d5004cc authored by Berend Hasselman on 05 September 2009, 00:00:00 UTC, committed by Gabor Csardi on 05 September 2009, 00:00:00 UTC
1 parent b291c25
liqrup.f
subroutine liqrup(q,ldq,n,r,ldr,u,v,wk)
integer ldq,n,ldr
double precision q(ldq,*),r(ldr,*),u(*),v(*),wk(*)
c-----------------------------------------------------------------------------
c
c Arguments
c
c Inout Q Real(ldq,n) orthogonal matrix from QR
c In ldq Integer leading dimension of Q
c In n Integer order of Q and R
c Inout R Real(ldr,n) upper triangular matrix R from QR
c In ldr Integer leading dimension of R
c In u Real(c) vector u of size n
c In v Real(c) vector v of size n
c Out wk Real(c) workspace of size n
c
c on return
c
c Q Q is the matrix with orthonormal columns in a QR
c decomposition of the matrix B = A + u*v'
c
c R R is the upper triangular matrix in a QR
c decomposition of the matrix B = A + u*v'
c
c Description
c
c The matrices Q and R are a QR decomposition of a square matrix
c A = Q*R.
c Given Q and R, qrupdt computes a QR decomposition of the rank one
c modification B = A + u*trans(v) of A. Here u and v are vectors.
c
c Source : Algorithm 686
c Fortran subroutines for updating the QR decomposition
c ACM Transactions on Mathematical software, dec. 1990
c
c-----------------------------------------------------------------------------
c Local variables and functions
integer k,i
double precision c,s,tmp
double precision ddot
c calculate wk = trans(Q)*u
do 10 i=1,n
wk(i) = ddot(n,q(1,i),1,u,1)
10 continue
c zero components wk(n),wk(n-1)...wk(2)
c and apply rotators to R and Q.
do 20 k=n-1,1,-1
call nuvgiv(wk(k),wk(k+1),c,s)
call drot(n-k+1,r(k,k),ldr,r(k+1,k),ldr,c,s)
call drot(n ,q(1,k),1 ,q(1,k+1),1 ,c,s)
20 continue
c r(1,1:n) += wk(1)*v(1:n)
call daxpy(n,wk(1),v,1,r(1,1),ldr)
c R is of upper hessenberg form. Triangularize R.
do 40 k=1,n-1
call nuvgiv(r(k,k),r(k+1,k),c,s)
call drot(n-k,r(k,k+1),ldr,r(k+1,k+1),ldr,c,s)
call drot(n ,q(1,k) ,1 ,q(1,k+1) ,1 ,c,s)
40 continue
return
end
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...