program test_dgesdd_dgesvd

! gfortran -o test_dgesdd_dgesvd progdgesdd.f90 -L/usr/local/lib -llapack -lcblas -lf77blas -latlas

  implicit none
  integer, parameter:: n = 128
  integer, parameter:: lwork = 7*n**2 + 4*n
  double precision A(n,n), Acopy(n,n), U(n,n), VT(n,n), sigma(n),temp(n,n)
  double precision norm1_A, norm1_res, nothing
  integer iwork(8*n)
  double precision dwork(lwork)
  integer i, j, info, choice
  double precision dlange

  open(10, FILE = 'matgesdd.dat' , STATUS = 'old')
  do j=1,n
     do i=1,n
        read(10,*) A(i,j)
     enddo
  enddo
  close(10)
  Acopy = A

  write(*,fmt='(A,$)') ' dgesdd (choice=0) or dgesvd(choice=1), choice='
  read(*,*) choice

  if ( choice == 0 ) then
     call dgesdd('A',n,n,Acopy,n,sigma,U,n,VT,n,dwork,lwork,iwork,info)
  else
     call dgesvd('A','A',n,n,Acopy,n,sigma,U,n,VT,n,dwork,lwork,info)
  end if

  if ( info /= 0 ) then
     write(*,*) 'Problem in dgesdd, info = ', info
  else
     if ( choice == 0 ) then
        write(*,*) ' With dgesdd:'
     else
        write(*,*) ' With dgesvd:'
     endif
     write(*,*) 'Singular values:'
     do i = 1,n
        write(*,*) 'sigma(',i,') = ',sigma(i)
     enddo
     
     write(*,*) 'Test || U*diag(sigma)*VT - A ||_1 / || A ||_1 '
     !  compute || A ||_1
     norm1_A = dlange( '1', n, n, A, n, nothing)
     write(*,*) '|| A ||_1 = ', norm1_A
     
     !  compute U*diag(sigma)
     do j = 1,n
        call dscal(n, sigma(j), U(1,j), 1)
     enddo
     
     ! compute  A <- U*diag(sigma))*VT - A
     call dgemm('N','N',n,n,n,1.d0,U,n,VT,n,-1.d0,A,n)
     
     ! compute  || U*diag(sigma))*VT - A ||_1
     norm1_res = dlange( '1', n, n, A, n, nothing)
     write(*,*) '|| res ||_1 = ', norm1_res
     write(*,*)  '|| U*diag(sigma)*VT - A ||_1 / || A ||_1 = ', norm1_res/norm1_A
  endif

end program test_dgesdd_dgesvd