From 8d12cfd29afb400548d5bcdba8de5dc93b124e8d Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Fri, 27 Sep 2019 00:50:47 +0200 Subject: [PATCH 1/6] rewrite everything. testing needed. --- Libraries/Modules/maxent_stoch_mod.F90 | 72 ++++++++++++++++---------- 1 file changed, 44 insertions(+), 28 deletions(-) diff --git a/Libraries/Modules/maxent_stoch_mod.F90 b/Libraries/Modules/maxent_stoch_mod.F90 index 20b452a6..6a0cb491 100644 --- a/Libraries/Modules/maxent_stoch_mod.F90 +++ b/Libraries/Modules/maxent_stoch_mod.F90 @@ -23,7 +23,7 @@ Module MaxEnt_stoch_mod Implicit None Real (Kind=Kind(0.d0)), Dimension(:) :: XQMC, Xtau, Alpha_tot - Real (Kind=Kind(0.d0)), Dimension(:,:) :: COV + Real (Kind=Kind(0.d0)), Dimension(:,:), intent(in) :: COV Real (Kind=Kind(0.d0)), External :: XKER, Back_trans_Aom Real (Kind=Kind(0.d0)) :: CHISQ, OM_ST, OM_EN, Beta_1, Xmom1, Err Integer :: Nsweeps, NBins, Ngamma_1, Ndis_1, nw, nt1 @@ -34,11 +34,11 @@ Module MaxEnt_stoch_mod & io_error, io_error1, i, n, nc1 Real (Kind=Kind(0.d0)), Allocatable :: Xn_M_tot(:,:), En_M_tot(:), Xn_E_tot(:,:), En_E_tot(:), & & Xn_tot(:,:,:), En_tot(:) - Real (Kind=Kind(0.d0)), Allocatable :: G_Mean(:), Xn_m(:), Xn_e(:), Xn(:,:), Vhelp(:) + Real (Kind=Kind(0.d0)), Allocatable :: G_Mean(:), Xn_m(:), Xn_e(:), Xn(:,:) Real (Kind=Kind(0.d0)) :: En_M, X, Alpha, Acc_1, Acc_2, En, DeltaE, Ratio, D Real (Kind=Kind(0.d0)) :: Aom, om, XMAX, tau Real (Kind=Kind(0.d0)) :: CPUT, CPUTM - Integer :: ICPU_1, ICPU_2, N_P_SEC + Integer :: ICPU_1, ICPU_2, N_P_SEC, info Character (64) :: File_root, File1, File_conf, File_Aom Real (Kind=Kind(0.d0)), allocatable :: Xker_table(:,:), U(:,:), sigma(:) ! Space for moments. @@ -63,7 +63,7 @@ Module MaxEnt_stoch_mod ! Setup table for the Kernel Ndis_table = 50000 Dom = (OM_EN_1 - OM_ST_1)/dble(Ndis_table-1) - Allocate ( Xker_table(Ntau, Ndis_table) ) + Allocate ( Xker_table(Ntau, Ndis_table), xqmc1(Ntau) ) do nt = 1,Ntau do nw = 1,Ndis_table tau = xtau(nt) @@ -72,40 +72,56 @@ Module MaxEnt_stoch_mod enddo enddo ! Normalize data to have zeroth moment of unity. - xqmc = xqmc / XMOM1 + xqmc1 = xqmc / XMOM1! xqmc1 is passed via a global variable to MC() cov = cov / ((XMOM1)**2) ! Diagonalize the covariance - Allocate( U(ntau,ntau), Sigma(ntau), xqmc1(Ntau) ) + Allocate( U(ntau,ntau), Sigma(ntau) ) If ( Present(L_cov) ) then - Call Diag(cov,U,sigma) - ! Write(6,*) " Cov Used" + U = cov + Call dpotrf('U', ntau, U, ntau, info) ! cov = U^T U +! Call Diag(cov,U,sigma) + Write(6,*) " Cov Used" + alpha = 1.0 + call dtrsm('L', 'U', 'T', 'N', ntau, 1, alpha, U, ntau, xqmc1, ntau) + call dtrsm('L', 'U', 'T', 'N', ntau, Ndis_table, alpha, U, ntau, xker_table, ntau) else Write(6,*) "No Cov Used" - U = 0.d0 do nt = 1,ntau - U(nt,nt) = 1.d0 sigma(nt) = cov(nt,nt) enddo - endif - do nt = 1,ntau - sigma(nt) = sqrt(sigma(nt)) - enddo - xqmc1 = 0.d0 - do nt1 = 1,ntau + do nt = 1,ntau - xqmc1(nt1) = xqmc1(nt1) + xqmc(nt)*U(nt,nt1) + sigma(nt) = sqrt(sigma(nt)) enddo - xqmc1(nt1) = xqmc1(nt1)/sigma(nt1) - enddo - ! Transform the Kernel - allocate ( Vhelp(Ntau) ) - do nw = 1,Ndis_table - Vhelp = 0.d0 - do nt1 = 1,Ntau - Vhelp(nt1) = Vhelp(nt1) + dot_product(Xker_table(:, nw), U(:,nt1)) - enddo - Xker_table(:, nw) = Vhelp/sigma - enddo + ! xqmc1 = 1/sigma * U^T * xqmc + do nt1 = 1,ntau + xqmc1(nt1) = xqmc1(nt1)/sigma(nt1) + enddo + ! Transform the Kernel + do nw = 1,Ndis_table + Xker_table(:, nw) = Xker_table(:, nw)/sigma + enddo + endif +! do nt = 1,ntau +! sigma(nt) = sqrt(sigma(nt)) +! enddo +! ! xqmc1 = 1/sigma * U^T * xqmc +! xqmc1 = 0.d0 +! do nt1 = 1,ntau +! do nt = 1,ntau ! x_i = x_i + sum_j U^T_{i, j} x_j +! xqmc1(nt1) = xqmc1(nt1) + xqmc(nt)*U(nt,nt1) +! enddo +! xqmc1(nt1) = xqmc1(nt1)/sigma(nt1) +! enddo +! ! Transform the Kernel +! allocate ( Vhelp(Ntau) ) +! do nw = 1,Ndis_table +! Vhelp = 0.d0 +! do nt1 = 1,Ntau +! Vhelp(nt1) = Vhelp(nt1) + dot_product(Xker_table(:, nw), U(:,nt1)) +! enddo +! Xker_table(:, nw) = Vhelp/sigma +! enddo deallocate( U, Sigma ) Allocate ( G_Mean(Ntau) ) G_mean = 0.d0 -- GitLab From 4b43e73ad0530f6a47396778bee4dd832446b5d7 Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Fri, 27 Sep 2019 00:53:31 +0200 Subject: [PATCH 2/6] first fix, seems like Fakher considers cov to be overwritable. --- Libraries/Modules/maxent_stoch_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Libraries/Modules/maxent_stoch_mod.F90 b/Libraries/Modules/maxent_stoch_mod.F90 index 6a0cb491..1ae870ff 100644 --- a/Libraries/Modules/maxent_stoch_mod.F90 +++ b/Libraries/Modules/maxent_stoch_mod.F90 @@ -23,7 +23,7 @@ Module MaxEnt_stoch_mod Implicit None Real (Kind=Kind(0.d0)), Dimension(:) :: XQMC, Xtau, Alpha_tot - Real (Kind=Kind(0.d0)), Dimension(:,:), intent(in) :: COV + Real (Kind=Kind(0.d0)), Dimension(:,:), intent(inout) :: COV Real (Kind=Kind(0.d0)), External :: XKER, Back_trans_Aom Real (Kind=Kind(0.d0)) :: CHISQ, OM_ST, OM_EN, Beta_1, Xmom1, Err Integer :: Nsweeps, NBins, Ngamma_1, Ndis_1, nw, nt1 -- GitLab From f753bc7dc433036ef72173ed88c867874a52e5e5 Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Fri, 27 Sep 2019 00:55:34 +0200 Subject: [PATCH 3/6] some error handling --- Libraries/Modules/maxent_stoch_mod.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Libraries/Modules/maxent_stoch_mod.F90 b/Libraries/Modules/maxent_stoch_mod.F90 index 1ae870ff..1d0ef7d6 100644 --- a/Libraries/Modules/maxent_stoch_mod.F90 +++ b/Libraries/Modules/maxent_stoch_mod.F90 @@ -79,6 +79,12 @@ Module MaxEnt_stoch_mod If ( Present(L_cov) ) then U = cov Call dpotrf('U', ntau, U, ntau, info) ! cov = U^T U +if (info < 0) then +write (*,*) "Argument wrong: ", info +endif +if (info > 0) then +write(*,*) "negative subminor found at ", info +endif ! Call Diag(cov,U,sigma) Write(6,*) " Cov Used" alpha = 1.0 -- GitLab From f0c8d8f6de9dc934de8bfeab33f602b6247351ea Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Fri, 4 Oct 2019 23:57:40 +0200 Subject: [PATCH 4/6] finish it. tests show that although the components of the resulting vector is different, their norm is identical --- Analysis/Max_SAC.F90 | 14 ++++++++------ Libraries/Modules/maxent_stoch_mod.F90 | 25 ++++++++++++++----------- 2 files changed, 22 insertions(+), 17 deletions(-) diff --git a/Analysis/Max_SAC.F90 b/Analysis/Max_SAC.F90 index b006db5a..1df64712 100644 --- a/Analysis/Max_SAC.F90 +++ b/Analysis/Max_SAC.F90 @@ -155,7 +155,7 @@ Stop end Select Ntau_old = Ntau - Call Rescale ( XCOV, XQMC,XTAU, Ntau_st, Ntau_en, Tolerance, NTAU) + ! Call Rescale ( XCOV, XQMC,XTAU, Ntau_st, Ntau_en, Tolerance, NTAU) Write(50,"('Data has been rescaled from Ntau ', I4,' to ', I4)") NTAU_old, Ntau If ( Ntau <= 4 ) then write(6,*) 'Not enough data!' @@ -171,7 +171,9 @@ XQMC_st = XQMC XTAU_st = XTAU - + do i = 1, ntau + xcov(i,i) = xcov(i,i) + 1.0 + enddo Allocate (Alpha_tot(N_alpha) ) do nt = 1,N_alpha alpha_tot(nt) = alpha_st*(R**(nt-1)) @@ -187,7 +189,7 @@ Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_ph, Back_Trans_ph, Beta, & & Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm) endif - ! Beware: Xqmc and cov are modified in the MaxEnt_stoch call. + ! Beware: cov is modified in the MaxEnt_stoch call. Case ("PP") If (N_cov == 1 ) then Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_pp, Back_Trans_pp, Beta, & @@ -196,7 +198,7 @@ Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_pp, Back_Trans_pp, Beta, & & Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm) endif - ! Beware: Xqmc and cov are modified in the MaxEnt_stoch call. + ! Beware: cov is modified in the MaxEnt_stoch call. Case ("P") If (N_cov == 1 ) then Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_p, Back_Trans_p, Beta, & @@ -205,7 +207,7 @@ Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_p, Back_Trans_p, Beta, & & Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm) endif - ! Beware: Xqmc and cov are modified in the MaxEnt_stoch call. + ! Beware: cov is modified in the MaxEnt_stoch call. Case ("T0") If (N_cov == 1 ) then Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_T0, Back_Trans_T0, Beta, & @@ -214,7 +216,7 @@ Call MaxEnt_stoch(XQMC, Xtau, Xcov, Xmom1, XKER_T0, Back_Trans_T0, Beta, & & Alpha_tot, Ngamma, OM_ST, OM_EN, Ndis, Nsweeps, NBins, NWarm) endif - ! Beware: Xqmc and cov are modified in the MaxEnt_stoch call. + ! Beware: cov is modified in the MaxEnt_stoch call. Case default Write(6,*) "Channel not yet implemented" Stop diff --git a/Libraries/Modules/maxent_stoch_mod.F90 b/Libraries/Modules/maxent_stoch_mod.F90 index 1d0ef7d6..11c6d439 100644 --- a/Libraries/Modules/maxent_stoch_mod.F90 +++ b/Libraries/Modules/maxent_stoch_mod.F90 @@ -21,8 +21,8 @@ Module MaxEnt_stoch_mod & Ngamma_1, OM_ST, OM_EN, Ndis_1, Nsweeps, NBins, NWarm, L_cov) Implicit None - - Real (Kind=Kind(0.d0)), Dimension(:) :: XQMC, Xtau, Alpha_tot + + Real (Kind=Kind(0.d0)), Dimension(:), intent(in) :: XQMC, Xtau, Alpha_tot Real (Kind=Kind(0.d0)), Dimension(:,:), intent(inout) :: COV Real (Kind=Kind(0.d0)), External :: XKER, Back_trans_Aom Real (Kind=Kind(0.d0)) :: CHISQ, OM_ST, OM_EN, Beta_1, Xmom1, Err @@ -75,22 +75,24 @@ Module MaxEnt_stoch_mod xqmc1 = xqmc / XMOM1! xqmc1 is passed via a global variable to MC() cov = cov / ((XMOM1)**2) ! Diagonalize the covariance + Allocate( U(ntau,ntau), Sigma(ntau) ) If ( Present(L_cov) ) then U = cov Call dpotrf('U', ntau, U, ntau, info) ! cov = U^T U -if (info < 0) then -write (*,*) "Argument wrong: ", info -endif -if (info > 0) then -write(*,*) "negative subminor found at ", info -endif + if (info < 0) then + write (*,*) "Argument wrong: ", info + endif + if (info > 0) then + write(*,*) "negative subminor found at ", info + endif ! Call Diag(cov,U,sigma) - Write(6,*) " Cov Used" + Write(6,*) " Cov Used in chol - v2" alpha = 1.0 + call dtrsm('L', 'U', 'T', 'N', ntau, 1, alpha, U, ntau, xqmc1, ntau) call dtrsm('L', 'U', 'T', 'N', ntau, Ndis_table, alpha, U, ntau, xker_table, ntau) - else + else! branch tested Write(6,*) "No Cov Used" do nt = 1,ntau sigma(nt) = cov(nt,nt) @@ -666,7 +668,8 @@ endif !Implicit Integer (H-N) Implicit None Real (Kind=Kind(0.d0)), Dimension(:,:) :: Xn, Xker_table - Real (Kind=Kind(0.d0)), Dimension(:) :: Xtau, Xn_m + Real (Kind=Kind(0.d0)), Dimension(:), intent(in) :: Xtau + Real (Kind=Kind(0.d0)), Dimension(:) :: Xn_m Real (Kind=Kind(0.d0)) :: Alpha, En_m, s, ratio, A_gamma, Z_gamma, Acc_1, Acc_2 Integer :: NSweeps, nl, Lambda_max, ng1, ng2 !Local -- GitLab From cfaaf280bc00a01492c5103bc1eff3ff6959e7de Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Sat, 5 Oct 2019 14:26:17 +0200 Subject: [PATCH 5/6] remove leftover things from debugging. --- Analysis/Max_SAC.F90 | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/Analysis/Max_SAC.F90 b/Analysis/Max_SAC.F90 index 1df64712..c31c27ef 100644 --- a/Analysis/Max_SAC.F90 +++ b/Analysis/Max_SAC.F90 @@ -155,7 +155,7 @@ Stop end Select Ntau_old = Ntau - ! Call Rescale ( XCOV, XQMC,XTAU, Ntau_st, Ntau_en, Tolerance, NTAU) + Call Rescale ( XCOV, XQMC,XTAU, Ntau_st, Ntau_en, Tolerance, NTAU) Write(50,"('Data has been rescaled from Ntau ', I4,' to ', I4)") NTAU_old, Ntau If ( Ntau <= 4 ) then write(6,*) 'Not enough data!' @@ -163,17 +163,14 @@ Endif If ( nbin_qmc > 2*Ntau .and. N_cov == 0 ) Write(50,*) 'Consider using the covariance. You seem to have enough bins' If ( nbin_qmc < 2*Ntau .and. N_cov == 1 ) Write(50,*) 'You do not seem to have enough bins for a reliable estimate of the covariance ' - - + + ! Store Allocate ( XCOV_st(NTAU,NTAU), XQMC_st(NTAU),XTAU_st(NTAU) ) XCOV_st = XCOV XQMC_st = XQMC XTAU_st = XTAU - - do i = 1, ntau - xcov(i,i) = xcov(i,i) + 1.0 - enddo + Allocate (Alpha_tot(N_alpha) ) do nt = 1,N_alpha alpha_tot(nt) = alpha_st*(R**(nt-1)) @@ -412,9 +409,8 @@ Integer, allocatable :: List(:) Real (Kind=Kind(0.d0)), dimension(:,:), allocatable :: XCOV_st Real (Kind=Kind(0.d0)), dimension(:) , allocatable :: XQMC_st, XTAU_st - - - + + ! Count the number of elements ntau_new = 0 Do nt = ntau_st,ntau_en @@ -446,9 +442,6 @@ XQMC = XQMC_st XTAU = XTAU_st Deallocate (XCOV_st, XQMC_st,XTAU_st, List ) - - - + end Subroutine Rescale - -- GitLab From a2be8c6ed7b527288b29540c17d1ab5e37f9bccd Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Thu, 24 Oct 2019 20:48:53 +0200 Subject: [PATCH 6/6] fix typos --- Analysis/Max_SAC.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Analysis/Max_SAC.F90 b/Analysis/Max_SAC.F90 index c31c27ef..b5272716 100644 --- a/Analysis/Max_SAC.F90 +++ b/Analysis/Max_SAC.F90 @@ -35,7 +35,7 @@ !> ALF-project ! !> @brief -!> General wrapper for maxent. Handles particle, partile-hole, particle-particle, +!> General wrapper for maxent. Handles particle, particle-hole, particle-particle, !> channels as well as zero temperature. See documentation for details. !> ! -- GitLab