From c862e8a88d29fa199e1376727a65653ba9b19acc Mon Sep 17 00:00:00 2001 From: ares201005 Date: Fri, 22 Sep 2023 16:03:34 -0400 Subject: [PATCH] C interface for genz module C interface for implicit fermi module C interface for the reset of sp2 module Corresponding tests Bug fixes --- scripts/transform.py | 13 +- src/prg_c_interface.F90 | 525 ++++++++++++++++++++++++++++++--- src/prg_progress_mod.h | 376 ++++++++++++++++++++--- src/prg_pulaycomponent_mod.F90 | 15 +- tests/CMakeLists.txt | 16 + tests/src/main.F90 | 6 +- tests/src/main.c | 494 ++++++++++++++++++++++++++++++- 7 files changed, 1339 insertions(+), 106 deletions(-) diff --git a/scripts/transform.py b/scripts/transform.py index 581fa244..6bee0172 100644 --- a/scripts/transform.py +++ b/scripts/transform.py @@ -37,6 +37,9 @@ def get_public(fortran_code): def transform_fortran_file(filename): + + subroutine_re = re.compile(r'subroutine\s+([a-zA-Z0-9_]+)\s*\(([\s\S]*)\)', re.DOTALL) + # use to find the argument names of each subroutine argument_re = re.compile(r'(real\(8\)|real\(dp\)|real\(PREC\)|integer|logical|type\(bml_matrix_t\)|character\(\d+\)|character\(len=\*\))(,.+?)?\s+::\s+([a-zA-Z0-9_,\s\(\):]+)') @@ -65,6 +68,7 @@ def transform_fortran_file(filename): "type(bml_matrix_t)": "bml_matrix_t*", "integer": "int", "character(len=*)":"char*", + "character(100)":"char*", "character(10)":"char*", "character(20)":"char*", "character(50)":"char*", @@ -92,6 +96,11 @@ def transform_fortran_file(filename): line = line.replace("\n","").strip() line = line + content[k+1].strip() line = line.replace("&", "") + if "&" in content[k+1]: + line = line.replace("\n","").strip() + line = line + content[k+2].strip() + line = line.replace("&", "") + match = subroutine_re.match(line.strip()) if match: subroutine_name = match.group(1) @@ -103,8 +112,8 @@ def transform_fortran_file(filename): ispublic = subroutine_name in public_subroutines - print(f"\ndebug-zy: a subroutine {subroutine_name} is found!", ispublic) - print(f"debug-zy: its arguments are: {argument_names}") + print(f"\nA subroutine {subroutine_name} is found!", ispublic) + print(f"Its arguments are: {argument_names}") # skip the subroutine if it's not public if not ispublic: continue is_inside_subroutine = True diff --git a/src/prg_c_interface.F90 b/src/prg_c_interface.F90 index cea5f5ee..e559ca27 100644 --- a/src/prg_c_interface.F90 +++ b/src/prg_c_interface.F90 @@ -8,15 +8,18 @@ module prg_c_interface use prg_ewald_mod use prg_dos_mod use prg_genz_mod + use prg_graphsolver_mod use prg_normalize_mod use prg_system_mod use prg_openfiles_mod + use prg_parallel_mod use prg_progress_mod use prg_pulaycomponent_mod use prg_sp2_fermi_mod use prg_sp2_mod use prg_timer_mod use bml_types_m + use prg_implicit_fermi_mod implicit none @@ -34,6 +37,16 @@ module prg_c_interface contains + function string_c2f(c_str) result(f_str) + character(c_char), dimension(:), intent(in) :: c_str + character(len=size(c_str)) :: f_str + integer :: i + do i = 1, size(c_str) + if (c_str(i) == c_null_char) exit + f_str(i:i) = c_str(i) + end do + end function string_c2f + ! C wrapper subroutine subroutine prg_version_c() bind(C, name="prg_version") call prg_version() @@ -282,7 +295,7 @@ end subroutine prg_get_eigenvalues_c subroutine prg_check_idempotency_c(mat_bml_c, threshold, idempotency)& bind(C, name="prg_check_idempotency") real(c_double), value :: threshold - real(c_double), value :: idempotency + real(c_double), intent(inout) :: idempotency type(c_ptr), value :: mat_bml_c type(bml_matrix_t) :: mat_bml mat_bml%ptr = mat_bml_c @@ -333,12 +346,10 @@ subroutine Canon_DM_PRT_c(P1, H1, Nocc, T, Q, e, mu0, m, HDIM) bind(C, name="Can call Canon_DM_PRT(P1, H1, Nocc, T, Q, e, mu0, m, HDIM) end subroutine Canon_DM_PRT_c - !------------------------------------------------ ! end of prg_densitymatrix_mod !------------------------------------------------ - !------------------------------------------------ ! Beginning of prg_charges_mod !------------------------------------------------ @@ -422,6 +433,125 @@ end subroutine prg_get_hscf_v2_c ! End of prg_charges_mod !------------------------------------------------ + + !------------------------------------------------ + ! prg_implicit_fermi_mod + !------------------------------------------------ + subroutine prg_implicit_fermi_save_inverse_c(Inv_bml_c, h_bml_c, p_bml_c, nsteps, nocc, mu,& + beta, occErrLimit, threshold, tol, SCF_IT, occiter, totns)& + bind(C, name="prg_implicit_fermi_save_inverse") + type(c_ptr), value :: h_bml_c + type(bml_matrix_t) :: h_bml + type(c_ptr), value :: p_bml_c + type(bml_matrix_t) :: p_bml + type(c_ptr), target :: Inv_bml_c(nsteps) + type(bml_matrix_t) :: Inv_bml(nsteps) + integer(c_int), value :: nsteps + integer(c_int), value :: SCF_IT + real(c_double), value :: nocc + real(c_double), value :: threshold + real(c_double), value :: tol + real(c_double), value :: occErrLimit + real(c_double), value :: beta + real(c_double), intent(inout) :: mu + integer(c_int), value :: occiter + integer(c_int), value :: totns + integer :: i + h_bml%ptr = h_bml_c + p_bml%ptr = p_bml_c + + ! Use the transfer intrinsic to move data from C pointer array to Fortran derived type array. + Inv_bml = transfer(Inv_bml_c, Inv_bml) + call prg_implicit_fermi_save_inverse(Inv_bml, h_bml, p_bml, nsteps, nocc, mu, beta,& + occErrLimit, threshold, tol, SCF_IT, occiter, totns) + + end subroutine prg_implicit_fermi_save_inverse_c + + subroutine prg_implicit_fermi_c(h_bml_c, p_bml_c, nsteps, k, nocc, mu, beta, method, osteps,& + occErrLimit, threshold, tol) bind(C, name="prg_implicit_fermi") + type(c_ptr), value :: h_bml_c + type(bml_matrix_t) :: h_bml + type(c_ptr), value :: p_bml_c + type(bml_matrix_t) :: p_bml + integer(c_int), value :: osteps + integer(c_int), value :: nsteps + integer(c_int), value :: method + integer(c_int), value :: k + real(c_double), value :: nocc + real(c_double), value :: threshold + real(c_double), value :: tol + real(c_double), value :: occErrLimit + real(c_double), value :: beta + real(c_double), value :: mu + h_bml%ptr = h_bml_c + p_bml%ptr = p_bml_c + call prg_implicit_fermi(h_bml, p_bml, nsteps, k, nocc, mu, beta, method, osteps,& + occErrLimit, threshold, tol) + end subroutine prg_implicit_fermi_c + + subroutine prg_implicit_fermi_zero_c(h_bml_c, p_bml_c, nsteps, mu, method, threshold, tol)& + bind(C, name="prg_implicit_fermi_zero") + type(c_ptr), value :: h_bml_c + type(bml_matrix_t) :: h_bml + type(c_ptr), value :: p_bml_c + type(bml_matrix_t) :: p_bml + integer(c_int), value :: nsteps + integer(c_int), value :: method + real(c_double), value :: mu + real(c_double), value :: threshold + real(c_double), value :: tol + h_bml%ptr = h_bml_c + p_bml%ptr = p_bml_c + call prg_implicit_fermi_zero(h_bml, p_bml, nsteps, mu, method, threshold, tol) + end subroutine prg_implicit_fermi_zero_c + + subroutine prg_implicit_fermi_first_order_response_c(H0_bml_c, H1_bml_c, P0_bml_c, P1_bml_c,& + Inv_bml_c, nsteps, mu0, beta, nocc, threshold)& + bind(C, name="prg_implicit_fermi_first_order_response") + type(c_ptr), value :: H0_bml_c + type(bml_matrix_t) :: H0_bml + type(c_ptr), value :: H1_bml_c + type(bml_matrix_t) :: H1_bml + type(c_ptr), target :: Inv_bml_c(nsteps) + type(bml_matrix_t) :: Inv_bml(nsteps) + type(c_ptr), value :: P0_bml_c + type(bml_matrix_t) :: P0_bml + type(c_ptr), value :: P1_bml_c + type(bml_matrix_t) :: P1_bml + real(c_double), value :: mu0 + real(c_double), value :: threshold + real(c_double), value :: beta + real(c_double), value :: nocc + integer(c_int), value :: nsteps + H0_bml%ptr = H0_bml_c + H1_bml%ptr = H1_bml_c + P0_bml%ptr = P0_bml_c + P1_bml%ptr = P1_bml_c + Inv_bml = transfer(Inv_bml_c, Inv_bml) + call prg_implicit_fermi_first_order_response(H0_bml, H1_bml, P0_bml, P1_bml, Inv_bml, nsteps,& + mu0, beta, nocc, threshold) + end subroutine prg_implicit_fermi_first_order_response_c + + subroutine prg_test_density_matrix_c(ham_bml_c, p_bml_c, beta, mu, nocc, osteps, occErrLimit, threshold) bind(C, name="prg_test_density_matrix") + type(c_ptr), value :: ham_bml_c + type(bml_matrix_t) :: ham_bml + type(c_ptr), value :: p_bml_c + type(bml_matrix_t) :: p_bml + real(c_double), value :: beta + real(c_double), value :: nocc + real(c_double), value :: occErrLimit + real(c_double), value :: threshold + real(c_double), value :: mu + integer(c_int), value :: osteps + ham_bml%ptr = ham_bml_c + p_bml%ptr = p_bml_c + call prg_test_density_matrix(ham_bml, p_bml, beta, mu, nocc, osteps, occErrLimit, threshold) + end subroutine prg_test_density_matrix_c + + !------------------------------------------------ + ! End of prg_implicit_fermi_mod + !------------------------------------------------ + !------------------------------------------------ ! prg_chebyshev_mod !------------------------------------------------ @@ -635,6 +765,83 @@ subroutine Ewald_Real_Space_c(COULOMBV, FCOUL, I, RX, RY, RZ, LBox, DELTAQ, U, E COULACC, TIMERATIO, nnRx, nnRy, nnRz, nrnnlist, nnType, HDIM, Max_Nr_Neigh) end subroutine Ewald_Real_Space_c + !------------------------------------------------ + ! prg_genz_mod + !------------------------------------------------ + + subroutine prg_init_ZSPmat_c(igenz, zk1_bml_c, zk2_bml_c, zk3_bml_c, zk4_bml_c,& + zk5_bml_c, zk6_bml_c, norb, bml_type, bml_element_type) & + bind(C, name="prg_init_ZSPmat") + integer(c_int), value :: norb + integer(c_int), value :: igenz + character(c_char), value :: bml_type + character(c_char), value :: bml_element_type + type(c_ptr), value :: zk1_bml_c + type(bml_matrix_t) :: zk1_bml + type(c_ptr), value :: zk2_bml_c + type(bml_matrix_t) :: zk2_bml + type(c_ptr), value :: zk3_bml_c + type(bml_matrix_t) :: zk3_bml + type(c_ptr), value :: zk4_bml_c + type(bml_matrix_t) :: zk4_bml + type(c_ptr), value :: zk5_bml_c + type(bml_matrix_t) :: zk5_bml + type(c_ptr), value :: zk6_bml_c + type(bml_matrix_t) :: zk6_bml + zk1_bml%ptr = zk1_bml_c + zk2_bml%ptr = zk2_bml_c + zk3_bml%ptr = zk3_bml_c + zk4_bml%ptr = zk4_bml_c + zk5_bml%ptr = zk5_bml_c + zk6_bml%ptr = zk6_bml_c + call prg_init_ZSPmat(igenz, zk1_bml, zk2_bml, zk3_bml, zk4_bml, zk5_bml, zk6_bml, norb, bml_type, bml_element_type) + end subroutine prg_init_ZSPmat_c + + subroutine prg_genz_sp_initialz0_c(smat_bml_c, zmat_bml_c, norb, mdim, bml_type_f, threshold)& + bind(C, name="prg_genz_sp_initialz0") + character(c_char), value :: bml_type_f + integer(c_int), value :: mdim + integer(c_int), value :: norb + real(c_double), value :: threshold + type(c_ptr), value :: zmat_bml_c + type(bml_matrix_t) :: zmat_bml + type(c_ptr), value :: smat_bml_c + type(bml_matrix_t) :: smat_bml + zmat_bml%ptr = zmat_bml_c + smat_bml%ptr = smat_bml_c + call prg_genz_sp_initialz0(smat_bml, zmat_bml, norb, mdim, bml_type_f, threshold) + end subroutine prg_genz_sp_initialz0_c + + subroutine prg_genz_sp_initial_zmat_c(smat_bml_c, zmat_bml_c, norb, mdim, bml_type_f, threshold)& + bind(C, name="prg_genz_sp_initial_zmat") + character(c_char), value :: bml_type_f + integer(c_int), value :: mdim + integer(c_int), value :: norb + real(c_double), value :: threshold + type(c_ptr), value :: zmat_bml_c + type(bml_matrix_t) :: zmat_bml + type(c_ptr), value :: smat_bml_c + type(bml_matrix_t) :: smat_bml + zmat_bml%ptr = zmat_bml_c + smat_bml%ptr = smat_bml_c + call prg_genz_sp_initial_zmat(smat_bml, zmat_bml, norb, mdim, bml_type_f, threshold) + end subroutine prg_genz_sp_initial_zmat_c + + subroutine prg_genz_sp_ref_c(smat_bml_c, zmat_bml_c, nref, norb, bml_type, threshold)& + bind(C, name="prg_genz_sp_ref") + integer(c_int), value :: norb + integer(c_int), value :: nref + character(c_char), value :: bml_type + real(c_double), value :: threshold + type(c_ptr), value :: smat_bml_c + type(c_ptr), value :: zmat_bml_c + type(bml_matrix_t) :: smat_bml + type(bml_matrix_t) :: zmat_bml + smat_bml%ptr = smat_bml_c + zmat_bml%ptr = zmat_bml_c + call prg_genz_sp_ref(smat_bml, zmat_bml, nref, norb, bml_type, threshold) + end subroutine prg_genz_sp_ref_c + !------------------------------------------------ ! prg_normalize_mod !------------------------------------------------ @@ -679,27 +886,40 @@ subroutine prg_normalize_cheb_c(h_bml_c, mu, emin, emax, alpha, scaledmu)& call prg_normalize_cheb(h_bml, mu, emin, emax, alpha, scaledmu) end subroutine prg_normalize_cheb_c + !------------------------------------------------ + ! prg_sys_mod + !------------------------------------------------ + + subroutine prg_get_nameandext_c(fullfilename, filename, ext) bind(C, name="prg_get_nameandext") + character(c_char), value :: fullfilename + character(c_char), value :: filename + character(c_char), value :: ext + call prg_get_nameandext(fullfilename, filename, ext) + end subroutine prg_get_nameandext_c + !------------------------------------------------ ! prg_sp2_mod !------------------------------------------------ - subroutine prg_sp2_fermi_init_c(h_bml_c, nsteps, nocc, tscale, threshold, occErrLimit, traceLimit, & - x_bml_c, mu, beta, h1, hN, sgnlist) bind(C, name="prg_sp2_fermi_init") + subroutine prg_sp2_fermi_init_c(h_bml_c, nsteps, nocc, tscale, threshold, occErrLimit,& + traceLimit, x_bml_c, mu, beta, h1, hN, sgnlist, verbose) & + bind(C, name="prg_sp2_fermi_init") type(c_ptr), value :: h_bml_c - type(bml_matrix_t) :: h_bml type(c_ptr), value :: x_bml_c + type(bml_matrix_t) :: h_bml type(bml_matrix_t) :: x_bml integer(c_int), value :: nsteps - integer(c_int), target :: sgnlist(nsteps) + integer(c_int), intent(inout) :: sgnlist(nsteps) + integer(c_int), value :: verbose real(c_double), value :: nocc real(c_double), value :: tscale real(c_double), value :: threshold real(c_double), value :: occErrLimit real(c_double), value :: traceLimit - real(c_double), value :: mu - real(c_double), value :: beta - real(c_double), value :: h1 - real(c_double), value :: hN + real(c_double), intent(inout) :: mu + real(c_double), intent(inout) :: beta + real(c_double), intent(inout) :: h1 + real(c_double), intent(inout) :: hN h_bml%ptr = h_bml_c x_bml%ptr = x_bml_c call prg_sp2_fermi_init(h_bml, nsteps, nocc, tscale, threshold, occErrLimit, traceLimit,& @@ -713,18 +933,18 @@ subroutine prg_sp2_fermi_init_norecs_c(h_bml_c, nsteps, nocc, tscale, threshold, type(bml_matrix_t) :: h_bml type(c_ptr), value :: x_bml_c type(bml_matrix_t) :: x_bml - integer(c_int), value :: nsteps - integer(c_int), target :: sgnlist(nsteps) + integer(c_int), intent(inout) :: nsteps + integer(c_int), intent(inout) :: sgnlist(nsteps) real(c_double), value :: nocc real(c_double), value :: tscale real(c_double), value :: threshold real(c_double), value :: occErrLimit real(c_double), value :: traceLimit - real(c_double), value :: mu - real(c_double), value :: beta - real(c_double), value :: h1 - real(c_double), value :: hN - integer(c_int), optional :: verbose + real(c_double), intent(inout) :: mu + real(c_double), intent(inout) :: beta + real(c_double), intent(inout) :: h1 + real(c_double), intent(inout) :: hN + integer(c_int), value :: verbose h_bml%ptr = h_bml_c x_bml%ptr = x_bml_c call prg_sp2_fermi_init_norecs(h_bml, nsteps, nocc, tscale, threshold, occErrLimit,& @@ -744,10 +964,10 @@ subroutine prg_sp2_fermi_c(h_bml_c, osteps, nsteps, nocc, mu, beta, h1, hN, sgnl real(c_double), value :: threshold real(c_double), value :: eps real(c_double), value :: traceLimit - real(c_double), value :: beta - real(c_double), value :: h1 - real(c_double), value :: hN - real(c_double), value :: mu + real(c_double), intent(inout) :: beta + real(c_double), intent(inout) :: h1 + real(c_double), intent(inout) :: hN + real(c_double), intent(inout) :: mu h_bml%ptr = h_bml_c x_bml%ptr = x_bml_c call prg_sp2_fermi(h_bml, osteps, nsteps, nocc, mu, beta, h1, hN, sgnlist, threshold,& @@ -770,6 +990,82 @@ subroutine prg_sp2_entropy_function_c(mu, h1, hN, nsteps, sgnlist, GG, ee)& deallocate(GG_tmp, ee_tmp) end subroutine prg_sp2_entropy_function_c + !------------------------------------------------ + ! prg_pulaycomponent_mod + !------------------------------------------------ + subroutine prg_PulayComponent0_c(rho_bml_c, ham_bml_c, pcm_bml_c, threshold, M, & + verbose) bind(C, name="prg_PulayComponent0") + type(c_ptr), value :: rho_bml_c + type(bml_matrix_t) :: rho_bml + type(c_ptr), value :: ham_bml_c + type(bml_matrix_t) :: ham_bml + type(c_ptr), value :: pcm_bml_c + type(bml_matrix_t) :: pcm_bml + integer(c_int), value :: M + integer(c_int), value :: verbose + real(c_double), value :: threshold + rho_bml%ptr = rho_bml_c + ham_bml%ptr = ham_bml_c + pcm_bml%ptr = pcm_bml_c + call prg_PulayComponent0(rho_bml, ham_bml, pcm_bml, threshold, M, verbose) + end subroutine prg_PulayComponent0_c + + subroutine prg_PulayComponentT_c(rho_bml_c, ham_bml_c, zmat_bml_c, pcm_bml_c, threshold, M,& + bml_type, verbose) bind(C, name="prg_PulayComponentT") + type(c_ptr), value :: rho_bml_c + type(bml_matrix_t) :: rho_bml + type(c_ptr), value :: ham_bml_c + type(bml_matrix_t) :: ham_bml + type(c_ptr), value :: zmat_bml_c + type(bml_matrix_t) :: zmat_bml + type(c_ptr), value :: pcm_bml_c + type(bml_matrix_t) :: pcm_bml + integer(c_int), value :: M + integer(c_int), value :: verbose + real(c_double), value :: threshold + character(c_char), value :: bml_type + rho_bml%ptr = rho_bml_c + ham_bml%ptr = ham_bml_c + zmat_bml%ptr = zmat_bml_c + pcm_bml%ptr = pcm_bml_c + call prg_PulayComponentT(rho_bml, ham_bml, zmat_bml, pcm_bml, threshold, M, bml_type, verbose) + end subroutine prg_PulayComponentT_c + + subroutine prg_get_pulayforce_c(nats, zmat_bml_c, ham_bml_c, rho_bml_c, dSx_bml_c, dSy_bml_c,& + dSz_bml_c, hindex_out, fpul_out, threshold) bind(C, name="prg_get_pulayforce") + integer(c_int), value :: nats + type(c_ptr), value :: dSx_bml_c + type(bml_matrix_t) :: dSx_bml + type(c_ptr), value :: dSy_bml_c + type(bml_matrix_t) :: dSy_bml + type(c_ptr), value :: dSz_bml_c + type(bml_matrix_t) :: dSz_bml + type(c_ptr), value :: rho_bml_c + type(bml_matrix_t) :: rho_bml + type(c_ptr), value :: ham_bml_c + type(bml_matrix_t) :: ham_bml + type(c_ptr), value :: zmat_bml_c + type(bml_matrix_t) :: zmat_bml + integer(c_int) :: hindex_out(2,nats) + real(c_double) :: fpul_out(3,nats) + real(c_double), allocatable :: fpul(:,:) + integer(c_int), allocatable :: hindex(:,:) + real(c_double), value :: threshold + dSx_bml%ptr = dSx_bml_c + dSy_bml%ptr = dSy_bml_c + dSz_bml%ptr = dSz_bml_c + rho_bml%ptr = rho_bml_c + ham_bml%ptr = ham_bml_c + zmat_bml%ptr = zmat_bml_c + + allocate(FPUL(3,nats), hindex(2,nats)) + call prg_get_pulayforce(nats, zmat_bml, ham_bml, rho_bml, dSx_bml, dSy_bml, dSz_bml, hindex, FPUL, threshold) + FPUL_OUT = FPUL + hindex_out = hindex + deallocate(fpul, hindex) + + end subroutine prg_get_pulayforce_c + !------------------------------------------------ ! prg_sp2_mod !------------------------------------------------ @@ -836,18 +1132,18 @@ subroutine prg_sp2_alg2_genseq_c(h_bml_c, rho_bml_c, threshold, bndfil, minsp2it idemtol, pp, icount, vv, verbose) bind(C, name="prg_sp2_alg2_genseq") integer(c_int), value :: minsp2iter integer(c_int), value :: maxsp2iter - integer(c_int), value :: icount - integer(c_int), target :: pp(:) + integer(c_int), intent(inout) :: icount real(c_double), value :: threshold real(c_double), value :: bndfil real(c_double), value :: idemtol - real(c_double), target :: vv(:) + integer(c_int), target :: pp(maxsp2iter) + real(c_double), target :: vv(maxsp2iter) character(c_char), value :: sp2conv + integer(c_int), value :: verbose type(c_ptr),value :: h_bml_c type(bml_matrix_t) :: h_bml type(c_ptr),value :: rho_bml_c type(bml_matrix_t) :: rho_bml - integer(c_int) :: verbose h_bml%ptr = h_bml_c rho_bml%ptr = rho_bml_c call prg_sp2_alg2_genseq(h_bml, rho_bml, threshold, bndfil, minsp2iter, maxsp2iter, sp2conv, idemtol,& @@ -856,7 +1152,7 @@ end subroutine prg_sp2_alg2_genseq_c subroutine prg_sp2_alg2_seq_c(h_bml_c, rho_bml_c, threshold, pp, icount, vv, verbose)& bind(C, name="prg_sp2_alg2_seq") - integer(c_int), value :: icount + integer(c_int), intent(inout) :: icount integer(c_int), target :: pp(:) real(c_double), value :: threshold real(c_double), target :: vv(:) @@ -864,25 +1160,27 @@ subroutine prg_sp2_alg2_seq_c(h_bml_c, rho_bml_c, threshold, pp, icount, vv, ver type(bml_matrix_t) :: h_bml type(c_ptr),value :: rho_bml_c type(bml_matrix_t) :: rho_bml - integer(c_int) :: verbose + integer(c_int), value :: verbose h_bml%ptr = h_bml_c rho_bml%ptr = rho_bml_c call prg_sp2_alg2_seq(h_bml, rho_bml, threshold, pp, icount, vv, verbose) end subroutine prg_sp2_alg2_seq_c - subroutine prg_prg_sp2_alg2_seq_inplace_c(rho_bml_c, threshold, pp, icount, vv, mineval, maxeval,& - verbose) bind(C, name="prg_prg_sp2_alg2_seq_inplace") - integer(c_int), value :: icount - integer(c_int), target :: pp(:) + subroutine prg_prg_sp2_alg2_seq_inplace_c(rho_bml_c, threshold, ppsize, pp, icount, vv, mineval, maxeval) & + !verbose) + bind(C, name="prg_prg_sp2_alg2_seq_inplace") + integer(c_int), intent(inout) :: icount + integer(c_int), value :: ppsize + integer(c_int), target :: pp(ppsize) real(c_double), value :: threshold - real(c_double), target :: vv(:) - real(c_double) :: mineval - real(c_double) :: maxeval + real(c_double), target :: vv(ppsize) + real(c_double), value :: mineval + real(c_double), value :: maxeval type(c_ptr),value :: rho_bml_c type(bml_matrix_t) :: rho_bml - integer(c_int) :: verbose + !integer(c_int), value :: verbose rho_bml%ptr = rho_bml_c - call prg_prg_sp2_alg2_seq_inplace(rho_bml, threshold, pp, icount, vv, mineval, maxeval,verbose) + call prg_prg_sp2_alg2_seq_inplace(rho_bml, threshold, pp, icount, vv, mineval, maxeval) !verbose) end subroutine prg_prg_sp2_alg2_seq_inplace_c subroutine prg_sp2_alg1_c(h_bml_c, rho_bml_c, threshold, bndfil, minsp2iter, maxsp2iter,& @@ -905,9 +1203,29 @@ subroutine prg_sp2_alg1_c(h_bml_c, rho_bml_c, threshold, bndfil, minsp2iter, max idemtol, verbose) end subroutine prg_sp2_alg1_c + subroutine prg_sp2_alg1_genseq_c(h_bml_c, rho_bml_c, threshold, bndfil, minsp2iter, & + maxsp2iter, sp2conv, idemtol, pp, icount, vv) bind(C, name="prg_sp2_alg1_genseq") + integer(c_int), value :: minsp2iter + integer(c_int), value :: maxsp2iter + integer(c_int), intent(inout) :: icount + real(c_double), value :: threshold + real(c_double), value :: bndfil + real(c_double), value :: idemtol + real(c_double), target :: vv(maxsp2iter) + integer(c_int), target :: pp(maxsp2iter) + character(c_char), value :: sp2conv + type(c_ptr),value :: h_bml_c + type(bml_matrix_t) :: h_bml + type(c_ptr),value :: rho_bml_c + type(bml_matrix_t) :: rho_bml + h_bml%ptr = h_bml_c + rho_bml%ptr = rho_bml_c + call prg_sp2_alg1_genseq(h_bml, rho_bml, threshold, bndfil, minsp2iter, maxsp2iter, sp2conv, idemtol, pp, icount, vv) + end subroutine prg_sp2_alg1_genseq_c + subroutine prg_sp2_alg1_seq_c(h_bml_c, rho_bml_c, threshold, pp, icount, vv)& bind(C, name="prg_sp2_alg1_seq") - integer(c_int), value :: icount + integer(c_int), intent(inout) :: icount integer(c_int), target :: pp(:) real(c_double), value :: threshold real(c_double), target :: vv(:) @@ -920,23 +1238,25 @@ subroutine prg_sp2_alg1_seq_c(h_bml_c, rho_bml_c, threshold, pp, icount, vv)& call prg_sp2_alg1_seq(h_bml, rho_bml, threshold, pp, icount, vv) end subroutine prg_sp2_alg1_seq_c - subroutine prg_prg_sp2_alg1_seq_inplace_c(rho_bml_c, threshold, pp, icount, vv, mineval, maxeval)& + subroutine prg_prg_sp2_alg1_seq_inplace_c(rho_bml_c, threshold, ppsize, pp, icount, vv, mineval, maxeval)& bind(C, name="prg_prg_sp2_alg1_seq_inplace") - integer(c_int), value :: icount - integer(c_int), target :: pp(:) + integer(c_int), value :: ppsize + integer(c_int), intent(inout) :: icount + integer(c_int), target :: pp(ppsize) real(c_double), value :: threshold - real(c_double), target :: vv(:) + real(c_double), target :: vv(ppsize) real(c_double), value :: mineval real(c_double), value :: maxeval type(c_ptr),value :: rho_bml_c type(bml_matrix_t) :: rho_bml rho_bml%ptr = rho_bml_c call prg_prg_sp2_alg1_seq_inplace(rho_bml, threshold, pp, icount, vv, mineval, maxeval) + end subroutine prg_prg_sp2_alg1_seq_inplace_c subroutine prg_sp2_submatrix_c(ham_bml_c, rho_bml_c, threshold, pp, icount, vv, mineval, maxeval,& core_size) bind(C, name="prg_sp2_submatrix") - integer(c_int), value :: icount + integer(c_int), intent(inout) :: icount integer(c_int), target :: pp(:) integer(c_int), value :: core_size real(c_double), value :: threshold @@ -954,7 +1274,7 @@ end subroutine prg_sp2_submatrix_c subroutine prg_sp2_submatrix_inplace_c(rho_bml_c, threshold, pp, icount, vv, mineval, maxeval,& core_size) bind(C, name="prg_sp2_submatrix_inplace") - integer(c_int), value :: icount + integer(c_int), intent(inout) :: icount integer(c_int), target :: pp(:) integer(c_int), value :: core_size real(c_double), value :: threshold @@ -1021,5 +1341,122 @@ subroutine prg_open_file_to_read_c(io, name) bind(C, name="prg_open_file_to_read call prg_open_file_to_read(io, name) end subroutine prg_open_file_to_read_c + !------------------------------------------------ + ! prg_parallel_mod + !------------------------------------------------ + subroutine prg_initParallel_c() bind(C, name="prg_initParallel") + call prg_initParallel() + end subroutine prg_initParallel_c + + subroutine prg_shutdownParallel_c() bind(C, name="prg_shutdownParallel") + call prg_shutdownParallel() + end subroutine prg_shutdownParallel_c + + subroutine prg_barrierParallel_c() bind(C, name="prg_barrierParallel") + call prg_barrierParallel() + end subroutine prg_barrierParallel_c + + subroutine sendReceiveParallel_c(sendBuf, sendLen, dest, recvBuf, recvLen, source, nreceived) bind(C, name="sendReceiveParallel") + integer(c_int), value :: sendLen + integer(c_int), value :: recvLen + real(c_double), target :: sendBuf(sendLen) + real(c_double), intent(out) :: recvBuf(recvLen) + integer(c_int), value :: dest + integer(c_int), value :: source + integer(c_int), value :: nreceived + call sendReceiveParallel(sendBuf, sendLen, dest, recvBuf, recvLen, source, nreceived) + end subroutine sendReceiveParallel_c + + subroutine isendParallel_c(sendBuf, sendLen, dest) bind(C, name="isendParallel") + real(c_double),target :: sendBuf(sendLen) + integer(c_int),value :: sendLen + integer(c_int),value :: dest + call isendParallel(sendBuf, sendLen, dest) + end subroutine isendParallel_c + + subroutine sendParallel_c(sendBuf, sendLen, dest) bind(C, name="sendParallel") + real(c_double),target :: sendBuf(sendLen) + integer(c_int),value :: sendLen + integer(c_int),value :: dest + call sendParallel(sendBuf, sendLen, dest) + end subroutine sendParallel_c + + subroutine prg_iprg_recvParallel_c(recvBuf, recvLen, rind) bind(C, name="prg_iprg_recvParallel") + real(c_double), target :: recvBuf(recvLen) + integer(c_int),value :: recvLen + integer(c_int), value :: rind + call prg_iprg_recvParallel(recvBuf, recvLen, rind) + end subroutine prg_iprg_recvParallel_c + + subroutine prg_recvParallel_c(recvBuf, recvLen) bind(C, name="prg_recvParallel") + real(c_double), target :: recvBuf(recvLen) + integer(c_int),value :: recvLen + call prg_recvParallel(recvBuf, recvLen) + end subroutine prg_recvParallel_c + + subroutine sumIntParallel_c(sendBuf, recvBuf, icount) bind(C, name="sumIntParallel") + integer(c_int),target :: sendBuf(icount) + integer(c_int), target :: recvBuf(icount) + integer(c_int),value :: icount + call sumIntParallel(sendBuf, recvBuf, icount) + end subroutine sumIntParallel_c + + subroutine sumRealParallel_c(sendBuf, recvBuf, icount) bind(C, name="sumRealParallel") + real(c_double),target :: sendBuf(icount) + real(c_double),intent(out) :: recvBuf(icount) + integer(c_int),value :: icount + call sumRealParallel(sendBuf, recvBuf, icount) + end subroutine sumRealParallel_c + + subroutine maxIntParallel_c(sendBuf, recvBuf, icount) bind(C, name="maxIntParallel") + integer(c_int),target :: sendBuf(icount) + integer(c_int),intent(out) :: recvBuf(icount) + integer(c_int),value :: icount + call maxIntParallel(sendBuf, recvBuf, icount) + end subroutine maxIntParallel_c + + subroutine maxRealParallel_c(sendBuf, recvBuf, icount) bind(C, name="maxRealParallel") + real(c_double),target :: sendBuf(icount) + real(c_double),intent(out) :: recvBuf(icount) + integer(c_int),value :: icount + call maxRealParallel(sendBuf, recvBuf, icount) + end subroutine maxRealParallel_c + + subroutine minIntParallel_c(sendBuf, recvBuf, icount) bind(C, name="minIntParallel") + integer(c_int),target :: sendBuf(icount) + integer(c_int),intent(out) :: recvBuf(icount) + integer(c_int),value :: icount + call minIntParallel(sendBuf, recvBuf, icount) + end subroutine minIntParallel_c + + subroutine minRealParallel_c(sendBuf, recvBuf, icount) bind(C, name="minRealParallel") + real(c_double),target :: sendBuf(icount) + real(c_double),intent(out) :: recvBuf(icount) + integer(c_int),value :: icount + call minRealParallel(sendBuf, recvBuf, icount) + end subroutine minRealParallel_c + + subroutine prg_minRealReduce_c(rvalue) bind(C, name="prg_minRealReduce") + real(c_double),value :: rvalue + call prg_minRealReduce(rvalue) + end subroutine prg_minRealReduce_c + + subroutine prg_maxRealReduce_c(rvalue) bind(C, name="prg_maxRealReduce") + real(c_double),value :: rvalue + call prg_maxRealReduce(rvalue) + end subroutine prg_maxRealReduce_c + + subroutine prg_maxIntReduce2_c(value1, value2) bind(C, name="prg_maxIntReduce2") + integer(c_int),value :: value1 + integer(c_int),value :: value2 + call prg_maxIntReduce2(value1, value2) + end subroutine prg_maxIntReduce2_c + + subroutine prg_sumIntReduce2_c(value1, value2) bind(C, name="prg_sumIntReduce2") + integer(c_int), value :: value1 + integer(c_int), value :: value2 + call prg_sumIntReduce2(value1, value2) + end subroutine prg_sumIntReduce2_c + end module prg_c_interface diff --git a/src/prg_progress_mod.h b/src/prg_progress_mod.h index a8c12a29..bb3e7563 100644 --- a/src/prg_progress_mod.h +++ b/src/prg_progress_mod.h @@ -122,7 +122,7 @@ extern "C" void prg_check_idempotency( bml_matrix_t * mat_bml, double threshold, - double idempotency); + double *idempotency); void prg_toEigenspace( bml_matrix_t * mat_bml, @@ -185,8 +185,99 @@ extern "C" int mdimin, double threshold); +//------prg_implicit_fermi_mod headers---------------------------------- + + void prg_implicit_fermi_save_inverse( + bml_matrix_t ** Inv_bml, + bml_matrix_t * h_bml, + bml_matrix_t * p_bml, + int nsteps, + double nocc, + double *mu, + double beta, + double occErrLimit, + double threshold, + double tol, + int SCF_IT, + int occiter, + int totns); + + void prg_implicit_fermi( + bml_matrix_t * h_bml, + bml_matrix_t * p_bml, + int nsteps, + int k, + double nocc, + double mu, + double beta, + int method, + int osteps, + double occErrLimit, + double threshold, + double tol); + + void prg_implicit_fermi_zero( + bml_matrix_t * h_bml, + bml_matrix_t * p_bml, + int nsteps, + double mu, + int method, + double threshold, + double tol); + + void prg_implicit_fermi_first_order_response( + bml_matrix_t * H0_bml, + bml_matrix_t * H1_bml, + bml_matrix_t * P0_bml, + bml_matrix_t * P1_bml, + bml_matrix_t ** Inv_bml, + int nsteps, + double mu0, + double beta, + double nocc, + double threshold); + + void prg_implicit_fermi_response( + bml_matrix_t * H0_bml, + bml_matrix_t * H1_bml, + bml_matrix_t * H2_bml, + bml_matrix_t * H3_bml, + bml_matrix_t * P0_bml, + bml_matrix_t * P1_bml, + bml_matrix_t * P2_bml, + bml_matrix_t * P3_bml, + int nsteps, + double mu0, + double mu, + double beta, + double nocc, + double occ_tol, + double lin_tol, + int order, + double threshold); + + void prg_finite_diff( + bml_matrix_t * H0_bml, + bml_matrix_t * H_list, + double mu0, + double mu_list, + double beta, + int order, + double lambda, + double h, + double threshold); + + void prg_test_density_matrix( + bml_matrix_t * ham_bml, + bml_matrix_t * p_bml, + double beta, + double mu, + double nocc, + int osteps, + double occErrLimit, + double threshold); + //------prg_chebyshev_mod headers -------------------------------------- -//void prg_parse_cheb(char* filename); void prg_build_density_cheb( bml_matrix_t * ham_bml, @@ -340,7 +431,8 @@ extern "C" int j2); // prg_sp2_fermi headers - void prg_sp2_fermi_init_norecs( + + void prg_sp2_fermi_init( bml_matrix_t * h_bml, int nsteps, double nocc, @@ -349,11 +441,26 @@ extern "C" double occErrLimit, double traceLimit, bml_matrix_t * x_bml, - double mu, - double beta, - double h1, - double hN, - int sgnlist, + double *mu, + double *beta, + double *h1, + double *hN, + int *sgnlist); + + void prg_sp2_fermi_init_norecs( + bml_matrix_t * h_bml, + int *nsteps, + double nocc, + double tscale, + double threshold, + double occErrLimit, + double traceLimit, + bml_matrix_t * x_bml, + double *mu, + double *beta, + double *h1, + double *hN, + int *sgnlist, int verbose); void prg_sp2_fermi( @@ -361,11 +468,11 @@ extern "C" int osteps, int nsteps, double nocc, - double mu, - double beta, - double h1, - double hN, - int sgnlist, + double *mu, + double *beta, + double *h1, + double *hN, + int *sgnlist, double threshold, double eps, double traceLimit, @@ -425,32 +532,32 @@ extern "C" int maxsp2iter, char *sp2conv, double idemtol, - int pp, - int icount, - double vv, + int *pp, + int *icount, + double *vv, int verbose); void prg_sp2_alg2_seq( bml_matrix_t * h_bml, bml_matrix_t * rho_bml, double threshold, - int pp, - int icount, - double vv, + int *pp, + int *icount, + double *vv, int verbose); - void prg_prg_sp2_alg2_seq_inplace( + void prg_sp2_alg1( bml_matrix_t * h_bml, bml_matrix_t * rho_bml, double threshold, - int pp, - int icount, - double vv, - double mineval, - double maxeval, + double bndfil, + int minsp2iter, + int maxsp2iter, + char *sp2conv, + double idemtol, int verbose); - void prg_sp2_alg1( + void prg_sp2_alg1_genseq( bml_matrix_t * h_bml, bml_matrix_t * rho_bml, double threshold, @@ -459,38 +566,53 @@ extern "C" int maxsp2iter, char *sp2conv, double idemtol, - int verbose); + int *pp, + int *icount, + double *vv); void prg_sp2_alg1_seq( bml_matrix_t * h_bml, bml_matrix_t * rho_bml, double threshold, - int pp, - int icount, - double vv); + int *pp, + int *icount, + double *vv); + + void prg_prg_sp2_alg2_seq_inplace( + bml_matrix_t * rho_bml, + double threshold, + int ppsize, + int *pp, + int *icount, + double *vv, + double mineval, + double maxeval); void prg_prg_sp2_alg1_seq_inplace( + bml_matrix_t * rho_bml, double threshold, - int pp, - int icount, - double vv, + int ppsize, + int *pp, + int *icount, + double *vv, double mineval, double maxeval); void prg_sp2_submatrix( double threshold, - int pp, - int icount, - double vv, + int *pp, + int *icount, + double *vv, double mineval, double maxeval, int core_size); void prg_sp2_submatrix_inplace( + bml_matrix_t * rho_bml, double threshold, - int pp, - int icount, - double vv, + int *pp, + int *icount, + double *vv, double mineval, double maxeval, int core_size); @@ -521,6 +643,182 @@ extern "C" void prg_print_date_and_time( char *tag); + void prg_get_pulayforce( + int nats, + bml_matrix_t * zmat_bml, + bml_matrix_t * ham_bml, + bml_matrix_t * rho_bml, + bml_matrix_t * dSx_bml, + bml_matrix_t * dSy_bml, + bml_matrix_t * dSz_bml, + int **hindex, + double **FPUL, + double threshold); + + +//-----prg_genz_mod headers ------------------------------------------- + //void prg_parse_ZSP(char* filename); + + void prg_buildzdiag( + bml_matrix_t * smat_bml, + bml_matrix_t * zmat_bml, + double threshold, + int mdimin, + char *bml_type, + int verbose); + + void prg_init_ZSPmat( + int igenz, + bml_matrix_t * zk1_bml, + bml_matrix_t * zk2_bml, + bml_matrix_t * zk3_bml, + bml_matrix_t * zk4_bml, + bml_matrix_t * zk5_bml, + bml_matrix_t * zk6_bml, + int norb, + char *bml_type, + char *bml_element_type); + + void prg_genz_sp_initialz0( + bml_matrix_t * smat_bml, + bml_matrix_t * zmat_bml, + int norb, + int mdim, + char *bml_type_f, + double threshold); + + void prg_genz_sp_initial_zmat( + bml_matrix_t * smat_bml, + bml_matrix_t * zmat_bml, + int norb, + int mdim, + char *bml_type_f, + double threshold); + + void prg_genz_sp_ref( + bml_matrix_t * smat_bml, + bml_matrix_t * zmat_bml, + int nref, + int norb, + char *bml_type, + double threshold); + +// prg_parallel_mod + + void prg_initParallel( + ); + + void prg_shutdownParallel( + ); + + void prg_barrierParallel( + ); + + void sendReceiveParallel( + double sendBuf, + int sendLen, + int dest, + double recvBuf, + int recvLen, + int source); + + void isendParallel( + double *sendBuf, + int sendLen, + int dest); + + void sendParallel( + double *sendBuf, + int sendLen, + int dest); + + void prg_iprg_recvParallel( + double *recvBuf, + int recvLen, + int rind); + + void prg_recvParallel( + double *recvBuf, + int recvLen); + + void sumIntParallel( + int *sendBuf, + int *recvBuf, + int icount); + + void sumRealParallel( + double *sendBuf, + double *recvBuf, + int icount); + + void maxIntParallel( + int *sendBuf, + int *recvBuf, + int icount); + + void maxRealParallel( + double *sendBuf, + double *recvBuf, + int icount); + + void minIntParallel( + double *sendBuf, + double *recvBuf, + int icount); + + void minRealParallel( + double *sendBuf, + double *recvBuf, + int icount); + + void prg_minRealReduce( + double rvalue); + + void prg_maxRealReduce( + double rvalue); + + void prg_maxIntReduce2( + int value1, + int value2); + + void prg_sumIntReduce2( + int value1, + int value2); + + void prg_wait( + ); + +// --- prg_pauly + + void prg_PulayComponent0( + bml_matrix_t * rho_bml, + bml_matrix_t * ham_bml, + bml_matrix_t * pcm_bml, + double threshold, + int M, + int verbose); + +//-----prg_graph_mod headers ------------------------------------------- + +//-----prg_graphsolver_mod headers ------------------------------------- + void prg_build_densityGP_T0( + bml_matrix_t * ham_bml, + bml_matrix_t * g_bml, + bml_matrix_t * rho_bml, + double threshold, + double bndfil, + double Ef, + int nparts, + int verbose); + + void prg_build_zmatGP( + bml_matrix_t * over_bml, + bml_matrix_t * g_bml, + bml_matrix_t * zmat_bml, + double threshold, + int nparts, + int verbose); + #ifdef __cplusplus } diff --git a/src/prg_pulaycomponent_mod.F90 b/src/prg_pulaycomponent_mod.F90 index 6a152a23..335bbb15 100644 --- a/src/prg_pulaycomponent_mod.F90 +++ b/src/prg_pulaycomponent_mod.F90 @@ -28,7 +28,8 @@ module prg_PulayComponent_mod !! \todo M and bml_type will have to be removed from the input parameter. !! subroutine prg_PulayComponent0(rho_bml,ham_bml,pcm_bml,threshold,M,& - &bml_type,verbose) + &verbose) + !&bml_type,verbose) implicit none @@ -39,23 +40,25 @@ subroutine prg_PulayComponent0(rho_bml,ham_bml,pcm_bml,threshold,M,& integer, intent(in) :: M integer :: nOrb, verbose real(dp), intent(in) :: threshold - character(20), intent(in) :: bml_type + !character(20), intent(in) :: bml_type + character(20) :: bml_type if(verbose.eq.2) write(*,*)"In prg_PulayComponent0 ..." nOrb = bml_get_N(rho_bml) + bml_type = bml_get_deep_type(rho_bml) call bml_zero_matrix(bml_type,bml_element_real,dp,nOrb ,nOrb,aux_bml) if(bml_get_N(pcm_bml).le.0)then !If pcm is not allocated call bml_zero_matrix(bml_type,bml_element_real,dp,nOrb,nOrb,pcm_bml) else - call bml_deallocate(pcm_bml) !If pcm is allocated then we set it to 0 - call bml_zero_matrix(bml_type,bml_element_real,dp,nOrb,nOrb,pcm_bml) + ! re-allocate will break the c pointer to pcm (for the c wrapper) + !call bml_deallocate(pcm_bml) !If pcm is allocated then we set it to 0 + !call bml_zero_matrix(bml_type,bml_element_real,dp,nOrb,nOrb,pcm_bml) + call bml_scale(0.0_dp, pcm_bml) endif - call bml_zero_matrix(bml_type,bml_element_real,dp,nOrb ,nOrb,pcm_bml) - call bml_multiply(rho_bml, ham_bml, aux_bml, 1.0d0, 1.0d0,threshold) !D*H call bml_multiply(aux_bml,rho_bml,pcm_bml, 1.0d0, 1.0d0,threshold) !(D*H)*D diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 230118ef..07599e2e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -108,11 +108,27 @@ progress_test(prg_density_c) progress_test(prg_density_T_c) progress_test(prg_density_T_fermi_c) progress_test(prg_density_T_fermi_grad_c) +# implicit fermi +progress_test(prg_implicit_fermi_c) +progress_test(prg_implicit_fermi_save_inverse_c) +# sp2 progress_test(prg_sp2_basic_c) progress_test(prg_sp2_alg1_dense_c) progress_test(prg_sp2_alg2_dense_c) progress_test(prg_sp2_alg1_ellpack_c) progress_test(prg_sp2_alg2_ellpack_c) +progress_test(prg_sp2_alg2_ellpack_poly_c) +progress_test(prg_sp2_alg1_seq_dense_c) +progress_test(prg_sp2_alg2_seq_dense_c) +progress_test(prg_sp2_alg1_seq_ellpack_c) +progress_test(prg_sp2_alg2_seq_ellpack_c) +progress_test(prg_sp2_alg1_seq_inplace_dense_c) +progress_test(prg_sp2_alg2_seq_inplace_dense_c) +progress_test(prg_sp2_alg1_seq_inplace_ellpack_c) +progress_test(prg_sp2_alg2_seq_inplace_ellpack_c) +progress_test(prg_sp2_fermi_dense_c) +progress_test(prg_pulaycomponent0_c) +# graph # LATTE modules progress_test(load_tbparms_latte) diff --git a/tests/src/main.F90 b/tests/src/main.F90 index 1eb58e16..d26c7365 100644 --- a/tests/src/main.F90 +++ b/tests/src/main.F90 @@ -156,7 +156,6 @@ program main call bml_scale(0.5_dp, rho_bml) call prg_check_idempotency(rho_bml,threshold,idempotency) write(*,*)"Idempotency for prg_build_density_T0",idempotency - write(*,*)"Fermi level:",mu if(idempotency.gt.1.0D-5)then write(*,*) "Idempotency is too high", idempotency error stop @@ -585,7 +584,7 @@ program main icount = 0 call prg_timer_start(sp2_timer) - call prg_sp2_alg2_genseq(ham_bml, rho_bml, threshold, bndfil, & + call prg_sp2_alg1_genseq(ham_bml, rho_bml, threshold, bndfil, & minsp2iter, maxsp2iter, sp2conv, sp2tol, & pp, icount, vv) call prg_timer_stop(sp2_timer) @@ -602,6 +601,7 @@ program main call bml_scale(0.5_dp, rho_bml) call prg_check_idempotency(rho_bml,threshold,idempotency) + write(6,*) "Idempotency for prg_sp2_alg1_seq_inplace_dense: ", idempotency if(idempotency.gt.idempotency_tol)then write(*,*) "Idempotency is too high", idempotency error stop @@ -981,7 +981,7 @@ program main call prg_timer_start(zdiag_timer) call prg_PulayComponent0(rho_bml,ham_bml,pcm_bml,threshold,mdim,& - &bml_type,verbose) + &verbose) call prg_timer_stop(zdiag_timer) write(*,*)"Pulay Component error ", error_calc diff --git a/tests/src/main.c b/tests/src/main.c index 4f60d50b..fe09a26d 100644 --- a/tests/src/main.c +++ b/tests/src/main.c @@ -1,6 +1,8 @@ +#include #include "string.h" #include "bml.h" +#include "bml_types.h" #include "prg_progress_mod.h" #include #include @@ -24,8 +26,9 @@ main( bml_matrix_t *zmat_bml, *zk1_bml, *zk2_bml; bml_matrix_t *zk3_bml, *zk4_bml, *zk5_bml, *zk6_bml; bml_matrix_t *nonortho_ham_bml; - bml_matrix_t *over_bml, *orthox_bml, *g_bml; - bml_matrix_t *aux_bml, *pcm_bml, *pcm_ref_bml, *inv_bml[10]; + bml_matrix_t *over_bml, *g_bml; + bml_matrix_t *aux_bml, *pcm, *pcm_ref, *inv_bml[10]; + bml_matrix_t *orthox_bml; bml_matrix_type_t matrix_type; bml_matrix_precision_t precision; @@ -36,10 +39,14 @@ main( char dummy[10][21]; // 20 char + null termination '\0' each char sp2conv[10] = "Rel"; - double threshold, gthreshold, idempotency; + double *trace; + double threshold, gthreshold; + double idempotency; double sp2tol, idempotency_tol; double bndfil, tscale, tracelimit, beta; double error_calc, error_tol, errlimit; + double maxeval, mineval, occerrlimit; + double beta0; matrix_type = dense; precision = double_real; @@ -80,10 +87,9 @@ main( prg_build_density_T0(norb, ham, rho, threshold, bndfil, eigenvalues); bml_scale(&scale_factor, rho, rho); - prg_check_idempotency(rho, threshold, idempotency); + prg_check_idempotency(rho, threshold, &idempotency); LOG_INFO("Idempotency for prg_build_density_T0: %.15e\n", idempotency); - if (idempotency > 1.0e-5) { printf("Idempotency is too high %f\n", idempotency); @@ -105,7 +111,7 @@ main( eigenvalues); bml_scale(&scale_factor, rho, rho); - prg_check_idempotency(rho, threshold, idempotency); + prg_check_idempotency(rho, threshold, &idempotency); LOG_INFO("Idempotency for prg_build_density_T: %.15e\n", idempotency); if (idempotency > 1.0e-5) @@ -132,6 +138,8 @@ main( prg_build_density_T_fermi(ham, rho, threshold, kbt, ef, 1, &drho); bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_build_density_T_fermi: %.15e\n", idempotency); if (idempotency > 1.0e-5) { printf("Idempotency is too high %f\n", idempotency); @@ -197,6 +205,76 @@ main( } } + else if (strcmp(test, "prg_implicit_fermi_c") == 0) + { + LOG_INFO + ("Testing the construction of the density matrix at KbT > 0 and at mu = Ef from implicit_fermi_mod \n"); + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + mu = 0.2; + beta = 4.0; + + rho1 = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + + prg_implicit_fermi(ham, rho1, 10, 2, 10.0, mu, beta, 0, 1, 1.0, threshold, 10e-8); + prg_test_density_matrix(ham, rho, beta, mu, 10.0, 1, 1.0, threshold); + + bml_add(rho1, rho, 1.0, -1.0, threshold); + + error_calc = bml_fnorm(rho1); + if (error_calc > 0.1) + { + printf("Error in Implicit Fermi expansion = %f", error_calc); + exit(EXIT_FAILURE); + } + + } + + else if (strcmp(test, "prg_implicit_fermi_save_inverse_c") == 0) + { + int norecs = 10; + int occiter; + int nsiter; + double nocc = 10.0; + + mu = 0.2; + beta = 4.0; + + for (i = 0; i < norecs; i++){ + inv_bml[i] = + bml_identity_matrix(matrix_type, precision, norb, norb, distrib_mode); + } + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + rho1 = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + + prg_implicit_fermi_save_inverse(inv_bml, ham, rho, norecs, nocc, + &mu, beta, 1e-4, threshold, 1e-5, 1, occiter, nsiter); + + prg_test_density_matrix(ham, rho1, beta, mu, nocc, 1, 1.0e-4, threshold); + printf("mu= %f \n", mu); + + bml_scale(&scale_factor, rho, rho); + bml_add(rho1, rho, 1.0, -1.0, threshold); + + error_calc = bml_fnorm(rho1); + if (error_calc > 0.1) + { + printf("Error in Implicit Fermi expansion save inverse= %f \n", error_calc); + exit(EXIT_FAILURE); + } + } + else if (strcmp(test, "prg_sp2_basic_c") == 0) { rho = @@ -207,12 +285,12 @@ main( prg_sp2_basic(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, sp2conv, sp2tol, verbose); bml_scale(&scale_factor, rho, rho); - prg_check_idempotency(rho, threshold, idempotency); + prg_check_idempotency(rho, threshold, &idempotency); LOG_INFO("Idempotency for prg_sp2_basic: %.15e\n", idempotency); if (idempotency > 1.0e-5) { - printf("Idempotency is too high %f\n", idempotency); + printf("Idempotency is too high %f\n", &idempotency); exit(EXIT_FAILURE); } } @@ -227,7 +305,7 @@ main( prg_sp2_alg1(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, sp2conv, sp2tol, verbose); bml_scale(&scale_factor, rho, rho); - prg_check_idempotency(rho, threshold, idempotency); + prg_check_idempotency(rho, threshold, &idempotency); LOG_INFO("Idempotency for prg_sp2_alg1_dense: %.15e\n", idempotency); if (idempotency > 1.0e-5) { @@ -246,7 +324,7 @@ main( prg_sp2_alg2(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, sp2conv, sp2tol, verbose); bml_scale(&scale_factor, rho, rho); - prg_check_idempotency(rho, threshold, idempotency); + prg_check_idempotency(rho, threshold, &idempotency); LOG_INFO("Idempotency for prg_sp2_alg2_dense: %.15e\n", idempotency); if (idempotency > 1.0e-5) { @@ -273,7 +351,7 @@ main( prg_sp2_alg1(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, sp2conv, sp2tol, verbose); bml_scale(&scale_factor, rho, rho); - prg_check_idempotency(rho, threshold, idempotency); + prg_check_idempotency(rho, threshold, &idempotency); LOG_INFO("Idempotency for prg_sp2_alg1_ellpack: %.15e\n", idempotency); if (idempotency > idempotency_tol) @@ -301,7 +379,7 @@ main( prg_sp2_alg2(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, sp2conv, sp2tol, verbose); bml_scale(&scale_factor, rho, rho); - prg_check_idempotency(rho, threshold, idempotency); + prg_check_idempotency(rho, threshold, &idempotency); LOG_INFO("Idempotency for prg_sp2_alg2_ellpack: %.15e\n", idempotency); if (idempotency > idempotency_tol) @@ -310,6 +388,398 @@ main( exit(EXIT_FAILURE); } } + + else if (strcmp(test, "prg_sp2_alg2_ellpack_poly_c") == 0) + { + matrix_type = ellpack; + idempotency_tol = 1.0e-2; + bndfil = 0.5; + norb = 6144; + mdim = 600; + threshold = 1.0e-5; + sp2tol = 1.0e-7; + + rho = + bml_zero_matrix(matrix_type, precision, norb, mdim, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, mdim, distrib_mode); + bml_read_bml_matrix(ham, "poly.512.mtx"); + + prg_sp2_alg2(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, + sp2conv, sp2tol, verbose); + + bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_sp2_alg2_ellpack_poly: %.15e\n", + idempotency); + + if (idempotency > idempotency_tol) + { + printf("Idempotency is too high %f\n", idempotency); + exit(EXIT_FAILURE); + } + } + + else if (strcmp(test, "prg_sp2_alg1_seq_dense_c") == 0) + { + matrix_type = dense; + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + int *pp = malloc(maxsp2iter * sizeof(int)); + double *vv = malloc(maxsp2iter * sizeof(double)); + int icount = 0; + + prg_sp2_alg1_genseq(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, + sp2conv, sp2tol, pp, &icount, vv); + + bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_sp2_alg1_seq_dense: %.15e\n", + idempotency); + if (idempotency > idempotency_tol) + { + printf("Idempotency is too high %f\n", idempotency); + exit(EXIT_FAILURE); + } + } + + else if (strcmp(test, "prg_sp2_alg2_seq_dense_c") == 0) + { + matrix_type = dense; + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + int *pp = malloc(maxsp2iter * sizeof(int)); + double *vv = malloc(maxsp2iter * sizeof(double)); + int icount = 0; + + prg_sp2_alg2_genseq(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, + sp2conv, sp2tol, pp, &icount, vv, verbose); + + bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_sp2_alg2_seq_dense: %.15e\n", + idempotency); + if (idempotency > idempotency_tol) + { + printf("Idempotency is too high %f\n", idempotency); + exit(EXIT_FAILURE); + } + } + + else if (strcmp(test, "prg_sp2_alg1_seq_ellpack_c") == 0) + { + matrix_type = ellpack; + idempotency_tol = 1.0e-6; + bndfil = 0.5; + norb = 6144; + mdim = 600; + threshold = 1.0e-9; + sp2tol = 1.0e-10; + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + int *pp = malloc(maxsp2iter * sizeof(int)); + double *vv = malloc(maxsp2iter * sizeof(double)); + int icount = 0; + + prg_sp2_alg1_genseq(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, + sp2conv, sp2tol, pp, &icount, vv); + + bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_sp2_alg1_seq_ellpack: %.15e\n", + idempotency); + if (idempotency > idempotency_tol) + { + printf("Idempotency is too high %f\n", idempotency); + exit(EXIT_FAILURE); + } + } + + else if (strcmp(test, "prg_sp2_alg2_seq_ellpack_c") == 0) + { + matrix_type = ellpack; + idempotency_tol = 1.0e-6; + bndfil = 0.5; + norb = 6144; + mdim = 600; + threshold = 1.0e-9; + sp2tol = 1.0e-10; + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + int *pp = malloc(maxsp2iter * sizeof(int)); + double *vv = malloc(maxsp2iter * sizeof(double)); + int icount = 0; + + prg_sp2_alg2_genseq(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, + sp2conv, sp2tol, pp, &icount, vv, verbose); + + bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_sp2_alg2_seq_ellpack: %.15e\n", + idempotency); + if (idempotency > idempotency_tol) + { + printf("Idempotency is too high %f\n", idempotency); + exit(EXIT_FAILURE); + } + } + + else if (strcmp(test, "prg_sp2_alg1_seq_inplace_dense_c") == 0) + { + matrix_type = dense; + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + int *pp = malloc(maxsp2iter * sizeof(int)); + double *vv = malloc(maxsp2iter * sizeof(double)); + double* gbnd = NULL; //malloc(2 * sizeof(double)); + int icount = 0; + + prg_sp2_alg1_genseq(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, + sp2conv, sp2tol, pp, &icount, vv); + + bml_copy(ham, rho); + gbnd = bml_gershgorin(rho); + prg_prg_sp2_alg1_seq_inplace(rho, threshold, maxsp2iter, pp, &icount, + vv, gbnd[0], gbnd[1]); + + bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_sp2_alg1_seq_inplace_dense: %.15e\n", + idempotency); + if (idempotency > idempotency_tol) + { + printf("Idempotency is too high %f\n", idempotency); + exit(EXIT_FAILURE); + } + } + + else if (strcmp(test, "prg_sp2_alg2_seq_inplace_dense_c") == 0) + { + matrix_type = dense; + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + int *pp = malloc(maxsp2iter * sizeof(int)); + double *vv = malloc(maxsp2iter * sizeof(double)); + double *gbnd = malloc(2 * sizeof(double)); + //double* gbnd = NULL; //malloc(2 * sizeof(double)); + int icount = 0; + + prg_sp2_alg2_genseq(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, + sp2conv, sp2tol, pp, &icount, vv, verbose); + + bml_copy(ham, rho); + gbnd = bml_gershgorin(rho); + prg_prg_sp2_alg2_seq_inplace(rho, threshold, maxsp2iter, pp, &icount, + vv, gbnd[0], gbnd[1]); + + bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_sp2_alg2_seq_inplace_dense: %.15e\n", + idempotency); + if (idempotency > idempotency_tol) + { + printf("Idempotency is too high %f\n", idempotency); + exit(EXIT_FAILURE); + } + } + + else if (strcmp(test, "prg_sp2_alg1_seq_inplace_ellpack_c") == 0) + { + matrix_type = ellpack; + idempotency_tol = 1.0e-6; + bndfil = 0.5; + norb = 6144; + mdim = 600; + threshold = 1.0e-9; + sp2tol = 1.0e-10; + + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + int *pp = malloc(maxsp2iter * sizeof(int)); + double *vv = malloc(maxsp2iter * sizeof(double)); + double *gbnd = malloc(2 * sizeof(double)); + int icount = 0; + + prg_sp2_alg1_genseq(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, + sp2conv, sp2tol, pp, &icount, vv); + + bml_copy(ham, rho); + gbnd = bml_gershgorin(rho); + prg_prg_sp2_alg1_seq_inplace(rho, threshold, maxsp2iter, pp, &icount, + vv, gbnd[0], gbnd[1]); + + bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_sp2_alg1_seq_inplace_ellpack: %.15e\n", + idempotency); + if (idempotency > idempotency_tol) + { + printf("Idempotency is too high %f\n", idempotency); + exit(EXIT_FAILURE); + } + } + + else if (strcmp(test, "prg_sp2_alg2_seq_inplace_ellpack_c") == 0) + { + matrix_type = ellpack; + idempotency_tol = 1.0e-6; + bndfil = 0.5; + norb = 6144; + mdim = 600; + threshold = 1.0e-9; + sp2tol = 1.0e-10; + + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + int *pp = malloc(maxsp2iter * sizeof(int)); + double *vv = malloc(maxsp2iter * sizeof(double)); + double *gbnd = malloc(2 * sizeof(double)); + int icount = 0; + + prg_sp2_alg2_genseq(ham, rho, threshold, bndfil, minsp2iter, maxsp2iter, + sp2conv, sp2tol, pp, &icount, vv, verbose); + + bml_copy(ham, rho); + gbnd = bml_gershgorin(rho); + prg_prg_sp2_alg2_seq_inplace(rho, threshold, maxsp2iter, pp, &icount, + vv, gbnd[0], gbnd[1]); + + bml_scale(&scale_factor, rho, rho); + prg_check_idempotency(rho, threshold, &idempotency); + LOG_INFO("Idempotency for prg_sp2_alg2_seq_inplace_ellpack: %.15e\n", + idempotency); + if (idempotency > idempotency_tol) + { + printf("Idempotency is too high %f\n", idempotency); + exit(EXIT_FAILURE); + } + } + + else if (strcmp(test, "prg_sp2_fermi_dense_c") == 0) + { + matrix_type = dense; + mdim = -1; + threshold = 1.0e-9; + sp2tol = 1.0e-10; + occerrlimit = 1.0e-9; + tracelimit = 1.0e-12; + mineval = 0.0; + maxeval = 0.0; + beta0 = 20.1; + tscale = 1.0; + + int norecs = 30; + int occsteps = 0; + double drho; + double eps = 1.0e-4; + double nocc = bndfil * norb; + int *signlist = malloc(norecs * sizeof(int)); + memset(signlist, 0, norecs * sizeof(int)); + + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + bml_read_bml_matrix(ham, "hamiltonian_ortho.mtx"); + + orthox_bml = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + + prg_sp2_fermi_init_norecs(ham, &norecs, nocc, tscale, threshold, + occerrlimit, tracelimit, orthox_bml, &mu, + &beta0, &mineval, &maxeval, &signlist[0], 10); + + prg_sp2_fermi_init(ham, norecs, nocc, tscale, threshold, occerrlimit, + tracelimit, orthox_bml, &mu, &beta0, &mineval, &maxeval, + &signlist[0]); + + prg_sp2_fermi(ham, occsteps, norecs, nocc, &mu, &beta0, &mineval, &maxeval, + &signlist[0], threshold, eps, tracelimit, orthox_bml); + + LOG_INFO("BETA %.15e\n", beta0); + LOG_INFO("ETEMP (eV) %.15e\n", 1.0/beta0); + LOG_INFO("NORECS USED %10d\n", norecs); + LOG_INFO("CHEMPOT %.15e\n", mu); + + kbt = 1.0 / beta0; + + prg_build_density_T_fermi(ham, rho, threshold, kbt, mu, 0, &drho); + bml_add(orthox_bml, rho, 2.0, -1.0, threshold); + error_calc = bml_fnorm(orthox_bml); + + if (error_calc > 0.01) + { + printf("Error in sp2 Fermi %f is too high\n", error_calc); + exit(EXIT_FAILURE); + } + + } + + else if(strcmp(test, "prg_pulaycomponent0_c") == 0){ + error_tol = 1.e-9; + matrix_type = dense; + + rho = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + ham = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + pcm = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + pcm_ref = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + zmat_bml = + bml_zero_matrix(matrix_type, precision, norb, norb, distrib_mode); + + bml_read_bml_matrix(zmat_bml,"zmatrix.mtx"); + bml_read_bml_matrix(ham,"hamiltonian.mtx"); + bml_read_bml_matrix(rho,"density.mtx"); + bml_read_bml_matrix(pcm_ref,"pcm.mtx"); + + prg_PulayComponent0(rho, ham, pcm, threshold, mdim, verbose); + bml_add(pcm, pcm_ref, 1.0, -1.0, threshold); + error_calc = bml_fnorm(pcm); + LOG_INFO("Pulay Component error %.12e\n", error_calc); + + if (error_calc > error_tol) + { + printf("Error in PulayComponent0 %f is too high\n", error_calc); + exit(EXIT_FAILURE); + } + + } + else { LOG_INFO("ERROR: unknown test \n");