diff --git a/Analysis/Max_SAC.F90 b/Analysis/Max_SAC.F90 index b8ccc2cad06ce4176a06bb9a840954f5359c8693..4f7b89fdd2345557e93eb34924d793e5fdecca5c 100644 --- a/Analysis/Max_SAC.F90 +++ b/Analysis/Max_SAC.F90 @@ -35,8 +35,8 @@ !> ALF-project ! !> @brief -!> General wrapper for maxent. Handles particle, partile-hole, particle-particle, -!> channels as well as zero temperature. See documentation for details. +!> General wrapper for maxent. Handles particle, partile-hole, particle-particle, +!> channels as well as zero temperature. See documentation for details. !> ! !-------------------------------------------------------------------- @@ -170,7 +170,6 @@ XQMC_st = XQMC XTAU_st = XTAU - Allocate (Alpha_tot(N_alpha) ) do nt = 1,N_alpha alpha_tot(nt) = alpha_st*(R**(nt-1)) @@ -186,7 +185,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, & @@ -195,7 +194,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, & @@ -204,7 +203,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, & @@ -410,8 +409,6 @@ 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 @@ -443,8 +440,4 @@ XQMC = XQMC_st XTAU = XTAU_st Deallocate (XCOV_st, XQMC_st,XTAU_st, List ) - - - - end Subroutine Rescale diff --git a/Libraries/Modules/maxent_stoch_mod.F90 b/Libraries/Modules/maxent_stoch_mod.F90 index 20b452a658c9a5846d6c33969a010f6e58cb5f29..11c6d4398fa610e5760ae1ff752e1530bae00a6b 100644 --- a/Libraries/Modules/maxent_stoch_mod.F90 +++ b/Libraries/Modules/maxent_stoch_mod.F90 @@ -21,9 +21,9 @@ 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(:,:) :: COV + + 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 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,64 @@ 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" - else + 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 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! branch tested 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 @@ -644,7 +668,8 @@ Module MaxEnt_stoch_mod !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