SUBROUTINE svdcmp_sp(a,w,v)
	USE nrtype; USE nrutil, ONLY : assert_eq,nrerror,outerprod
	USE nr, ONLY : pythag
	IMPLICIT NONE
	REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a
	REAL(SP), DIMENSION(:), INTENT(OUT) :: w
	REAL(SP), DIMENSION(:,:), INTENT(OUT) :: v
	INTEGER(I4B) :: i,its,j,k,l,m,n,nm
	REAL(SP) :: anorm,c,f,g,h,s,scale,x,y,z
	REAL(SP), DIMENSION(size(a,1)) :: tempm
	REAL(SP), DIMENSION(size(a,2)) :: rv1,tempn
	m=size(a,1)
	n=assert_eq(size(a,2),size(v,1),size(v,2),size(w),'svdcmp_sp')
	g=0.0
	scale=0.0
	do i=1,n
		l=i+1
		rv1(i)=scale*g
		g=0.0
		scale=0.0
		if (i <= m) then
			scale=sum(abs(a(i:m,i)))
			if (scale /= 0.0) then
				a(i:m,i)=a(i:m,i)/scale
				s=dot_product(a(i:m,i),a(i:m,i))
				f=a(i,i)
				g=-sign(sqrt(s),f)
				h=f*g-s
				a(i,i)=f-g
				tempn(l:n)=matmul(a(i:m,i),a(i:m,l:n))/h
				a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
				a(i:m,i)=scale*a(i:m,i)
			end if
		end if
		w(i)=scale*g
		g=0.0
		scale=0.0
		if ((i <= m) .and. (i /= n)) then
			scale=sum(abs(a(i,l:n)))
			if (scale /= 0.0) then
				a(i,l:n)=a(i,l:n)/scale
				s=dot_product(a(i,l:n),a(i,l:n))
				f=a(i,l)
				g=-sign(sqrt(s),f)
				h=f*g-s
				a(i,l)=f-g
				rv1(l:n)=a(i,l:n)/h
				tempm(l:m)=matmul(a(l:m,l:n),a(i,l:n))
				a(l:m,l:n)=a(l:m,l:n)+outerprod(tempm(l:m),rv1(l:n))
				a(i,l:n)=scale*a(i,l:n)
			end if
		end if
	end do
	anorm=maxval(abs(w)+abs(rv1))
	do i=n,1,-1
		if (i < n) then
			if (g /= 0.0) then
				v(l:n,i)=(a(i,l:n)/a(i,l))/g
				tempn(l:n)=matmul(a(i,l:n),v(l:n,l:n))
				v(l:n,l:n)=v(l:n,l:n)+outerprod(v(l:n,i),tempn(l:n))
			end if
			v(i,l:n)=0.0
			v(l:n,i)=0.0
		end if
		v(i,i)=1.0
		g=rv1(i)
		l=i
	end do
	do i=min(m,n),1,-1
		l=i+1
		g=w(i)
		a(i,l:n)=0.0
		if (g /= 0.0) then
			g=1.0_sp/g
			tempn(l:n)=(matmul(a(l:m,i),a(l:m,l:n))/a(i,i))*g
			a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
			a(i:m,i)=a(i:m,i)*g
		else
			a(i:m,i)=0.0
		end if
		a(i,i)=a(i,i)+1.0_sp
	end do
	do k=n,1,-1
		do its=1,30
			do l=k,1,-1
				nm=l-1
				if ((abs(rv1(l))+anorm) == anorm) exit
				if ((abs(w(nm))+anorm) == anorm) then
					c=0.0
					s=1.0
					do i=l,k
						f=s*rv1(i)
						rv1(i)=c*rv1(i)
						if ((abs(f)+anorm) == anorm) exit
						g=w(i)
						h=pythag(f,g)
						w(i)=h
						h=1.0_sp/h
						c= (g*h)
						s=-(f*h)
						tempm(1:m)=a(1:m,nm)
						a(1:m,nm)=a(1:m,nm)*c+a(1:m,i)*s
						a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
					end do
					exit
				end if
			end do
			z=w(k)
			if (l == k) then
				if (z < 0.0) then
					w(k)=-z
					v(1:n,k)=-v(1:n,k)
				end if
				exit
			end if
			if (its == 30) call nrerror('svdcmp_sp: no convergence in svdcmp')
			x=w(l)
			nm=k-1
			y=w(nm)
			g=rv1(nm)
			h=rv1(k)
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0_sp*h*y)
			g=pythag(f,1.0_sp)
			f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
			c=1.0
			s=1.0
			do j=l,nm
				i=j+1
				g=rv1(i)
				y=w(i)
				h=s*g
				g=c*g
				z=pythag(f,h)
				rv1(j)=z
				c=f/z
				s=h/z
				f= (x*c)+(g*s)
				g=-(x*s)+(g*c)
				h=y*s
				y=y*c
				tempn(1:n)=v(1:n,j)
				v(1:n,j)=v(1:n,j)*c+v(1:n,i)*s
				v(1:n,i)=-tempn(1:n)*s+v(1:n,i)*c
				z=pythag(f,h)
				w(j)=z
				if (z /= 0.0) then
					z=1.0_sp/z
					c=f*z
					s=h*z
				end if
				f= (c*g)+(s*y)
				x=-(s*g)+(c*y)
				tempm(1:m)=a(1:m,j)
				a(1:m,j)=a(1:m,j)*c+a(1:m,i)*s
				a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
			end do
			rv1(l)=0.0
			rv1(k)=f
			w(k)=x
		end do
	end do
	END SUBROUTINE svdcmp_sp

	SUBROUTINE svdcmp_dp(a,w,v)
	USE nrtype; USE nrutil, ONLY : assert_eq,nrerror,outerprod
	USE nr, ONLY : pythag
	IMPLICIT NONE
	REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: a
	REAL(DP), DIMENSION(:), INTENT(OUT) :: w
	REAL(DP), DIMENSION(:,:), INTENT(OUT) :: v
	INTEGER(I4B) :: i,its,j,k,l,m,n,nm
	REAL(DP) :: anorm,c,f,g,h,s,scale,x,y,z
	REAL(DP), DIMENSION(size(a,1)) :: tempm
	REAL(DP), DIMENSION(size(a,2)) :: rv1,tempn
	m=size(a,1)
	n=assert_eq(size(a,2),size(v,1),size(v,2),size(w),'svdcmp_dp')
	g=0.0
	scale=0.0
	do i=1,n
		l=i+1
		rv1(i)=scale*g
		g=0.0
		scale=0.0
		if (i <= m) then
			scale=sum(abs(a(i:m,i)))
			if (scale /= 0.0) then
				a(i:m,i)=a(i:m,i)/scale
				s=dot_product(a(i:m,i),a(i:m,i))
				f=a(i,i)
				g=-sign(sqrt(s),f)
				h=f*g-s
				a(i,i)=f-g
				tempn(l:n)=matmul(a(i:m,i),a(i:m,l:n))/h
				a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
				a(i:m,i)=scale*a(i:m,i)
			end if
		end if
		w(i)=scale*g
		g=0.0
		scale=0.0
		if ((i <= m) .and. (i /= n)) then
			scale=sum(abs(a(i,l:n)))
			if (scale /= 0.0) then
				a(i,l:n)=a(i,l:n)/scale
				s=dot_product(a(i,l:n),a(i,l:n))
				f=a(i,l)
				g=-sign(sqrt(s),f)
				h=f*g-s
				a(i,l)=f-g
				rv1(l:n)=a(i,l:n)/h
				tempm(l:m)=matmul(a(l:m,l:n),a(i,l:n))
				a(l:m,l:n)=a(l:m,l:n)+outerprod(tempm(l:m),rv1(l:n))
				a(i,l:n)=scale*a(i,l:n)
			end if
		end if
	end do
	anorm=maxval(abs(w)+abs(rv1))
	do i=n,1,-1
		if (i < n) then
			if (g /= 0.0) then
				v(l:n,i)=(a(i,l:n)/a(i,l))/g
				tempn(l:n)=matmul(a(i,l:n),v(l:n,l:n))
				v(l:n,l:n)=v(l:n,l:n)+outerprod(v(l:n,i),tempn(l:n))
			end if
			v(i,l:n)=0.0
			v(l:n,i)=0.0
		end if
		v(i,i)=1.0
		g=rv1(i)
		l=i
	end do
	do i=min(m,n),1,-1
		l=i+1
		g=w(i)
		a(i,l:n)=0.0
		if (g /= 0.0) then
			g=1.0_dp/g
			tempn(l:n)=(matmul(a(l:m,i),a(l:m,l:n))/a(i,i))*g
			a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
			a(i:m,i)=a(i:m,i)*g
		else
			a(i:m,i)=0.0
		end if
		a(i,i)=a(i,i)+1.0_dp
	end do
	do k=n,1,-1
		do its=1,30
			do l=k,1,-1
				nm=l-1
				if ((abs(rv1(l))+anorm) == anorm) exit
				if ((abs(w(nm))+anorm) == anorm) then
					c=0.0
					s=1.0
					do i=l,k
						f=s*rv1(i)
						rv1(i)=c*rv1(i)
						if ((abs(f)+anorm) == anorm) exit
						g=w(i)
						h=pythag(f,g)
						w(i)=h
						h=1.0_dp/h
						c= (g*h)
						s=-(f*h)
						tempm(1:m)=a(1:m,nm)
						a(1:m,nm)=a(1:m,nm)*c+a(1:m,i)*s
						a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
					end do
					exit
				end if
			end do
			z=w(k)
			if (l == k) then
				if (z < 0.0) then
					w(k)=-z
					v(1:n,k)=-v(1:n,k)
				end if
				exit
			end if
			if (its == 30) call nrerror('svdcmp_dp: no convergence in svdcmp')
			x=w(l)
			nm=k-1
			y=w(nm)
			g=rv1(nm)
			h=rv1(k)
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0_dp*h*y)
			g=pythag(f,1.0_dp)
			f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
			c=1.0
			s=1.0
			do j=l,nm
				i=j+1
				g=rv1(i)
				y=w(i)
				h=s*g
				g=c*g
				z=pythag(f,h)
				rv1(j)=z
				c=f/z
				s=h/z
				f= (x*c)+(g*s)
				g=-(x*s)+(g*c)
				h=y*s
				y=y*c
				tempn(1:n)=v(1:n,j)
				v(1:n,j)=v(1:n,j)*c+v(1:n,i)*s
				v(1:n,i)=-tempn(1:n)*s+v(1:n,i)*c
				z=pythag(f,h)
				w(j)=z
				if (z /= 0.0) then
					z=1.0_dp/z
					c=f*z
					s=h*z
				end if
				f= (c*g)+(s*y)
				x=-(s*g)+(c*y)
				tempm(1:m)=a(1:m,j)
				a(1:m,j)=a(1:m,j)*c+a(1:m,i)*s
				a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
			end do
			rv1(l)=0.0
			rv1(k)=f
			w(k)=x
		end do
	end do
	END SUBROUTINE svdcmp_dp