Skip to content
Snippets Groups Projects
Commit 79f69ff8 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

fixed the test

parent e76a04a1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 187 additions and 42 deletions
src/DLRA_mod.F90
+ 187
− 42
View file @ 79f69ff8
......@@ -17,42 +17,56 @@ module DLRA
! CALL SVD(array_ky_pj,singular_values)
END SUBROUTINE
SUBROUTINE test_SVD
SUBROUTINE test_svd
#ifdef TEST_SVD
! Program to perform Singular Value Decomposition (SVD)
! using LAPACK library
! Specify the dimensions of the input matrix A
INTEGER, PARAMETER :: m = 3, n = 2
INTEGER, PARAMETER :: lda = m
! Declare the input matrix A
REAL, DIMENSION(lda,n) :: A
! Specify the dimensions of the output matrices
INTEGER, PARAMETER :: ldu = m, ldvt = n
INTEGER, PARAMETER :: lwork = 5*n
REAL, DIMENSION(ldu,m) :: U
REAL, DIMENSION(n) :: S
REAL, DIMENSION(ldvt,n) :: VT
REAL, DIMENSION(lwork) :: work
INTEGER :: info, i,j
COMPLEX, DIMENSION(m,n) :: A
! OUTPUT
COMPLEX, DIMENSION(m,m) :: U
REAL, DIMENSION(MIN(m,n)) :: S
COMPLEX, DIMENSION(n,n) :: VT
! local variables
INTEGER :: lda, ldu, ldvt, info, lwork, i, j
COMPLEX, DIMENSION(:), ALLOCATABLE :: work
REAL, DIMENSION(:), ALLOCATABLE :: rwork
! Set the leading dimensions for the input and output arrays
lda = MAX(1, m)
ldu = MAX(1, m)
ldvt = MAX(1, n)
! Define the input matrix A
A = RESHAPE((/ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 /), SHAPE(A))
A = RESHAPE((/ (1.0,0.1), (2.0,0.2), (3.0,0.3), (4.0,0.4), (5.0,0.5), (6.0,0.6) /), SHAPE(A))
! Print the input matrix A
WRITE(*,*) 'Input matrix A = '
DO i = 1, m
WRITE(*,'(2X, 3F8.3)') (A(i,j), j=1,n)
WRITE(*,*) ('(',REAL(A(i,j)), AIMAG(A(i,j)),')', j=1,n)
END DO
! Compute the SVD of A using the LAPACK subroutine SGESVD
CALL SGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, info)
ALLOCATE(work(5*n), rwork(5*n))
! Compute the optimal workspace size
lwork = -1
CALL CGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
! Allocate memory for the workspace arrays
lwork = CEILING(REAL(work(1)), KIND=SELECTED_REAL_KIND(1, 6))
! Compute the SVD of A using the LAPACK subroutine CGESVD
CALL CGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
! Print the results
WRITE(*,*) 'U = '
DO i = 1, m
WRITE(*,'(6F8.3)') (U(i,j), j=1,m)
WRITE(*,*) ('(',REAL(U(i,j)), AIMAG(U(i,j)),')', j=1,m)
END DO
WRITE(*,*)
WRITE(*,*) 'S = '
......@@ -60,35 +74,166 @@ module DLRA
WRITE(*,*)
WRITE(*,*) 'VT = '
DO i = 1, n
WRITE(*,'(6F8.3)') (VT(i,j), j=1,n)
WRITE(*,*) ('(',REAL(VT(i,j)), AIMAG(VT(i,j)),')', j=1,n)
END DO
! Reconstruct A from its SVD
A = MATMUL(U, MATMUL(diagmat(S,m,n), TRANSPOSE(VT)))
A = MATMUL(U, MATMUL(diagmat(S,m,n), VT))
! Print the reconstructed matrix A
WRITE(*,*) 'Reconstructed matrix A = '
DO i = 1, m
WRITE(*,'(2X, 3F8.3)') (A(i,j), j=1,n)
WRITE(*,*) ('(',REAL(A(i,j)), AIMAG(A(i,j)),')', j=1,n)
END DO
print*, "this was a test of the SVD using LAPACK. End run."
stop
END SUBROUTINE test_SVD
FUNCTION diagmat(v, m, n) RESULT(A)
REAL, DIMENSION(:), INTENT(IN) :: v
INTEGER, INTENT(IN) :: m, n
REAL, DIMENSION(m,n) :: A
INTEGER :: i, j
A = 0.0 ! Initialize A to a zero matrix
#endif
END SUBROUTINE test_svd
DO i = 1, MIN(m,n)
A(i,i) = v(i) ! Set the diagonal elements of A to the values in v
END DO
! SUBROUTINE test_svd
! ! Program to perform Singular Value Decomposition (SVD)
! ! using LAPACK library
! ! Specify the dimensions of the input matrix A
! INTEGER, PARAMETER :: m = 3, n = 2
! INTEGER, PARAMETER :: lda = m
! ! Declare the input matrix A
! REAL, DIMENSION(lda,n) :: A
! ! Specify the dimensions of the output matrices
! INTEGER, PARAMETER :: ldu = m, ldvt = n
! INTEGER, PARAMETER :: lwork = 5*n
! REAL, DIMENSION(ldu,m) :: U
! REAL, DIMENSION(n) :: S
! REAL, DIMENSION(ldvt,n) :: VT
! REAL, DIMENSION(lwork) :: work
! INTEGER :: info,i,j
! ! Define the input matrix A
! A = RESHAPE((/ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 /), SHAPE(A))
! ! Compute the SVD of A using the LAPACK subroutine SGESVD
! CALL SGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, info)
! ! Print the results
! WRITE(*,*) 'U = '
! DO i = 1, m
! WRITE(*,'(6F8.3)') (U(i,j), j=1,m)
! END DO
! WRITE(*,*)
! WRITE(*,*) 'S = '
! WRITE(*,'(2F8.3)') (S(i), i=1,n)
! WRITE(*,*)
! WRITE(*,*) 'VT = '
! DO i = 1, n
! WRITE(*,'(6F8.3)') (VT(i,j), j=1,n)
! END DO
! ! Reconstruct A from its SVD
! A = MATMUL(U, MATMUL(diagmat(S,m,n), TRANSPOSE(VT)))
! ! Print the reconstructed matrix A
! WRITE(*,*) 'Reconstructed matrix A = '
! DO i = 1, m
! WRITE(*,'(2X, 3F8.3)') (A(i,j), j=1,n)
! END DO
! stop
! END SUBROUTINE test_svd
! SUBROUTINE test_svd
! IMPLICIT NONE
! INTEGER, PARAMETER :: m = 3, n = 2
! COMPLEX(xp), DIMENSION(m, n) :: A
! COMPLEX(xp), DIMENSION(m, m) :: U
! COMPLEX(xp), DIMENSION(n, n) :: VT
! REAL(xp), DIMENSION(MIN(m,n)) :: S
! INTEGER :: i, j
! ! Initialize A
! A = RESHAPE((/ (1.0_xp, 2.0_xp), (3.0_xp, 4.0_xp), (5.0_xp, 6.0_xp),&
! (1.5_xp, 2.5_xp), (3.5_xp, 4.5_xp), (5.5_xp, 6.5_xp) /), [m, n])
! ! Print input
! WRITE(*,*) "A = "
! DO i = 1, m
! WRITE(*,"(3F8.3)") (REAL(A(i,j)), AIMAG(A(i,j)), j=1,n)
! END DO
! WRITE(*,*)
! ! Call the SVD subroutine
! CALL svd(A, U, S, VT, m, n)
! ! Print the resultas
! WRITE(*,*) "U = "
! DO i = 1, m
! WRITE(*,"(3F8.3)") (REAL(U(i,j)), AIMAG(U(i,j)), j=1,m)
! END DO
! WRITE(*,*)
! WRITE(*,*) "S = ", S
! WRITE(*,*)
! WRITE(*,*) "VT = "
! DO i = 1, n
! WRITE(*,"(3F8.3)") (REAL(VT(i,j)), AIMAG(VT(i,j)), j=1,n)
! END DO
! END SUBROUTINE test_svd
SUBROUTINE svd(A, U, S, VT, m, n)
IMPLICIT NONE
! INPUT
COMPLEX(xp), DIMENSION(m,n), INTENT(IN) :: A
INTEGER, INTENT(IN) :: m, n
! OUTPUT
COMPLEX(xp), DIMENSION(m,m), INTENT(OUT) :: U
REAL(xp), DIMENSION(MIN(m,n)), INTENT(OUT) :: S
COMPLEX(xp), DIMENSION(n,n), INTENT(OUT) :: VT
! local variables
INTEGER :: lda, ldu, ldvt, info, lwork
COMPLEX(xp), DIMENSION(:), ALLOCATABLE :: work
REAL(xp), DIMENSION(:), ALLOCATABLE :: rwork
#ifdef LAPACKDIR
! Set the leading dimensions for the input and output arrays
lda = MAX(1, m)
ldu = MAX(1, m)
ldvt = MAX(1, n)
! Compute the optimal workspace size
lwork = -1
CALL ZGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
! Allocate memory for the workspace arrays
lwork = CEILING(REAL(work(1)), KIND=SELECTED_REAL_KIND(1, 6))
ALLOCATE(work(lwork), rwork(5*MIN(m,n)))
! Compute the SVD
CALL ZGESVD('A', 'A', m, n, A, lda, S, U, ldu, VT, ldvt, work, lwork, rwork, info)
! Free the workspace arrays
DEALLOCATE(work, rwork)
#endif
END SUBROUTINE svd
END FUNCTION diagmat
FUNCTION diagmat(S, m, n) RESULT(D)
IMPLICIT NONE
! INPUT
REAL, DIMENSION(:), INTENT(IN) :: S
INTEGER, INTENT(IN) :: m, n
! OUTPUT
COMPLEX, DIMENSION(m,n) :: D
! Local variables
INTEGER :: i, j
! Initialize the output array to zero
D = 0.0
! Fill the diagonal elements of the output array with the input vector S
DO i = 1, MIN(m,n)
D(i,i) = S(i)
END DO
END FUNCTION diagmat
end module DLRA
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment