module gramschmidt
implicit none
contains

subroutine printm(A)
real::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::v(:)
integer i
print "(10F9.4)",(v(i),i=1,size(v));
end subroutine printv

subroutine qrdec(a,r)
real::a(:,:),r(:,:)
integer::n,m,i,j
n=size(a,1)
m=size(a,2)
do i=1,m
	r(i,i)=sqrt(dot_product(a(:,i),a(:,i)))
	a(:,i)=a(:,i)/r(i,i)
	do j=i+1,m
		r(i,j)=dot_product(a(:,i),a(:,j))
		r(j,i)=0
		a(:,j)=a(:,j)-r(i,j)*a(:,i)
	end do
end do
end subroutine qrdec

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

subroutine qrsolve(A,r,b,x)
real::A(:,:),r(:,:),b(:),x(:)
real,allocatable::v(:)
integer::m,i
m=size(A,2)
v=matmul(transpose(A),b)
x(1:m)=v(1:m)
call qrbak(r,x)
end subroutine qrsolve

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

end module gramschmidt