From a6a3e5ddc50e4dee9c51248774607288f8c8ee7b Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Sat, 1 Apr 2023 16:16:20 +0200
Subject: [PATCH] First DLRA and LAPACK test

---
 Makefile                         |  9 +++-
 local/make.inc                   |  8 ++--
 src/DLRA_mod.F90                 | 78 +++++++++++++++++++++++++++++++-
 src/stepon.F90                   |  2 +
 testcases/zpinch_example/fort.90 |  0
 5 files changed, 89 insertions(+), 8 deletions(-)
 create mode 100644 testcases/zpinch_example/fort.90

diff --git a/Makefile b/Makefile
index bd218c82..ee63f350 100644
--- a/Makefile
+++ b/Makefile
@@ -93,7 +93,7 @@ $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/parallel_mod.o
 $(OBJDIR)/ppexit.o $(OBJDIR)/prec_const_mod.o \
 $(OBJDIR)/processing_mod.o $(OBJDIR)/readinputs.o $(OBJDIR)/restarts_mod.o \
 $(OBJDIR)/solve_EM_fields.o $(OBJDIR)/species_mod.o $(OBJDIR)/stepon.o $(OBJDIR)/tesend.o \
-$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
+$(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/DLRA_mod.o
 
 # To compile the executable
  compile: $(FOBJ)
@@ -282,7 +282,8 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
    $(OBJDIR)/basic_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/grid_mod.o \
 	 $(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o \
 	 $(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/solve_EM_fields.o\
-	 $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+	 $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
+	 $(OBJDIR)/DLRA_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
 
  $(OBJDIR)/tesend.o : src/tesend.F90 \
@@ -296,3 +297,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
  $(OBJDIR)/utility_mod.o : src/utility_mod.F90  \
    $(OBJDIR)/grid_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
+
+ $(OBJDIR)/DLRA_mod.o : src/DLRA_mod.F90  \
+   $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/DLRA_mod.F90 -o $@
diff --git a/local/make.inc b/local/make.inc
index 8be74cc8..6781cd5d 100644
--- a/local/make.inc
+++ b/local/make.inc
@@ -78,10 +78,10 @@ ifdef FFTW3DIR
 endif
 
 # Add lapack local lib
-ifdef LAPACKDIR
-      LIBS += -llapack -lblas
-      LDIRS += -L$(LAPACKDIR)/lib
-endif
+#ifdef LAPACKDIR
+      LIBS  += -llapack -lblas
+#      LDIRS += -L$(LAPACKDIR)/lib
+#endif
 
 # FM library
 ifdef FMDIR
diff --git a/src/DLRA_mod.F90 b/src/DLRA_mod.F90
index be925949..34f3b779 100644
--- a/src/DLRA_mod.F90
+++ b/src/DLRA_mod.F90
@@ -1,9 +1,8 @@
 module DLRA
     USE prec_const
-    USE lapack
     implicit none
 
-    PUBLIC :: filter_singular_value_ky_pj
+    PUBLIC :: filter_singular_value_ky_pj, test_SVD
 
     CONTAINS
 
@@ -17,4 +16,79 @@ module DLRA
         ! Singular value decomposition
         ! CALL SVD(array_ky_pj,singular_values)
     END SUBROUTINE
+
+    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))
+
+        WRITE(*,*) 'Input matrix A = '
+
+        DO i = 1, m
+            WRITE(*,'(2X, 3F8.3)') (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)
+
+        ! 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
+
+        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
+      
+        DO i = 1, MIN(m,n)
+          A(i,i) = v(i)   ! Set the diagonal elements of A to the values in v
+        END DO
+      
+      END FUNCTION diagmat
+
 end module DLRA
\ No newline at end of file
diff --git a/src/stepon.F90 b/src/stepon.F90
index d9cbe1fc..89c61ab7 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -8,6 +8,7 @@ SUBROUTINE stepon
    use mpi,                   ONLY: MPI_COMM_WORLD
    USE time_integration,      ONLY: ntimelevel
    USE prec_const,            ONLY: xp
+   USE DLRA,                  ONLY: test_SVD
    IMPLICIT NONE
 
    INTEGER :: num_step, ierr
@@ -56,6 +57,7 @@ SUBROUTINE stepon
 
       !! TEST SINGULAR VALUE DECOMPOSITION
       ! CALL filter_singular_value_ky_pj(nsv,moments)
+      ! CALL test_SVD
       
       IF( nlend ) EXIT ! exit do loop
 
diff --git a/testcases/zpinch_example/fort.90 b/testcases/zpinch_example/fort.90
new file mode 100644
index 00000000..e69de29b
-- 
GitLab