! Performs Gram-Schmidt process
! Replaces X(n, m) by an orthonormal matrix U(n, m)
! "m" vectors with "n" coordinates
subroutine Gram_Schmidt(n, m, X, tol)

implicit real(8) (a-h, o-z)

complex(8) X(n,m)
complex(8) p(n)	! Projection vector =SUM[ <x(k+1),u(i)>*u(i) ]
complex(8) z(n)	! Vector normal to u(1),...,u(k): z=x(k+1) -p
complex(8) u(n)	! working vector = current "u"

! Normalize vector #1:
z =X(:,1)
z =z/vNorm2(n,z)
X(:,1) = z

Jvector_Do: DO j=2, m
! Calculate "j"th vector = ortonormal to 1, 2, (j-1).
k_Do: do k=1, j-1
p =(0.,0.)
	i_Do: do i=1, k			! Calculate p =SUM:
	u =X(:, i)
	p =p +dot_product(u, X(:,j))*u
	! End "p"
	enddo i_Do
enddo k_Do

z =X(:, j) -p

if(vNorm2(n,z) >=TOL) then
X(:, j) =z/vNorm2(n,z)
else
X(:, j) =(0.,0.)		!z
endif
ENDDO Jvector_Do

!write(1,*) 'G-S: X-orthonormal:'
!do i =1,n
!write(1,*) X(i,1:m)
!enddo
END subroutine Gram_Schmidt
