module qrgivens
implicit none
contains

subroutine printm(A)
real*4::A(:,:)
integer i,j
do i=1,size(A,1); print "(10F9.4)",(A(i,j),j=1,size(A,2)); end do
end subroutine printm

subroutine printv(v)
real*4::v(:)
integer i
print "(10F9.4)",(v(i),i=1,size(v));
end subroutine printv

subroutine qrdec(A)
real*4::A(:,:),theta,s,c,xq,xp
integer::n,m,p,q,k
n=size(A,1)
m=size(A,2)
do q=1,m
do p=q+1,n
	theta=atan2(A(p,q),A(q,q))
	do k=q,m
		xq=A(q,k)
		xp=A(p,k)
		A(p,k)= cos(theta)*xp-sin(theta)*xq
		A(q,k)= sin(theta)*xp+cos(theta)*xq
	end do
	A(p,q)=theta
end do
end do
end subroutine qrdec

subroutine qrgetr(A,R)
real*4::A(:,:),R(:,:)
integer i,j
R(:,:)=0
do i=1,size(A,2)
do j=i,size(A,2)
	R(i,j)=A(i,j)
end do
end do
end subroutine qrgetr

subroutine qrgetq(A,Q)
real*4::A(:,:),Q(:,:)
real*4,allocatable::e(:)
integer i,k
allocate(e(size(A,1)))
do i=1,size(A,1)
	e(:)=0
	e(i)=1
	call qrqt(A,e)
	Q(i,:)=e(:)
end do
end subroutine qrgetq

subroutine qrbak(A,x)
real*4::A(:,:),x(:),s
integer::m,k,i
m=size(A,2)
do i=m,1,-1
	s=0
	do k=i+1,m
		s=s+A(i,k)*x(k)
	end do
	x(i)=(x(i)-s)/A(i,i)
end do
end subroutine qrbak

subroutine qrqt(A,v)
real*4::A(:,:),v(:),vp,vq,theta
integer::n,m,p,q
n=size(A,1)
m=size(A,2)
do q=1,m
do p=q+1,n
	theta=A(p,q)
	vp=v(p)
	vq=v(q)
	v(p)=+cos(theta)*vp-sin(theta)*vq
	v(q)=+sin(theta)*vp+cos(theta)*vq
end do
end do
end subroutine qrqt

function qrsolve(A,b) result (x)
real*4::A(:,:),b(:),x(size(A,2))
real*4,allocatable::v(:)
integer::n,m,i
n=size(A,1)
m=size(A,2)
allocate(v(n))
v(:)=b(:)
call qrqt(A,v)
x(1:m)=v(1:m)
call qrbak(A,x)
end function qrsolve

function qrinv(a) result (b)
real*4 a(:,:),b(size(a,2),size(a,2))
real*4 e(size(a,1))
integer i
e(:)=0
do i=1,size(a,2)
	e(i)=1
	b(:,i)=qrsolve(A,e)
	e(i)=0
end do
end function qrinv

end module qrgivens