Loading doc/src/config.dox +2 −1 Original line number Diff line number Diff line Loading @@ -986,7 +986,8 @@ INPUT_FILE_ENCODING = FILE_PATTERNS = *.cpp \ *.f \ *.h \ *.md *.md \ *.py # The RECURSIVE tag can be used to specify whether or not subdirectories should # be searched for input files as well. Loading src/scripts/pycompare→src/scripts/pycompare.py +47 −21 Original line number Diff line number Diff line #!/bin/python ## @package pycompare # Script to perform output consistency tests # # Comparing the numeric output can be rendered hard by the amount of information # contained in a typical output file and the necessity to determine whether a # difference is actually significant or just caused by numeric noise hitting # negligible values. The task of this script is to compare two output files, in # the assumption that they were written by the FORTRAN and the C++ versions of # the code and to flag all the possible inconsistencies according to various # severity levels (namely: NOISE, WARNING, and ERROR). import re from math import log10 from sys import argv ## \cond number_reg = re.compile(r'-?[0-9]\.[0-9]+E[-+][0-9]{2,2}') ## \endcond ## \brief Main execution code # # `main()` is the function that handles the creation of the script configuration # and the execution of the comparison. It returns an integer value corresponding # to the number of detected error-level inconsistencies. # # \returns errors: `int` Number of detected error-level inconsistencies. def main(): config = parse_arguments() errors, warnings, noisy = (0, 0, 0) Loading @@ -32,6 +52,7 @@ def main(): )) return errors ## \brief Perform the comparison of two files. def compare_files(config): mismatch_count = { 'errors': 0, Loading Loading @@ -78,8 +99,14 @@ def compare_files(config): l_file.write(" </body>\n") l_file.write("</html>\n") l_file.close() else: mismatch_count['errors'] = len(c_lines) print("ERROR: {0:s} and {1:s} have different numbers of lines!".format( config['fortran_file_name'], config['c_file_name'] )) return mismatch_count ## \brief Perform the comparison of two file lines. def compare_lines(f_line, c_line, config, line_num=0, num_len=1, log_file=None): errors = 0 warnings = 0 Loading Loading @@ -156,26 +183,21 @@ def compare_lines(f_line, c_line, config, line_num=0, num_len=1, log_file=None): log_file.write(log_line + "</code></pre></div>\n") return (errors, warnings, noisy) ## \brief Determine the severity of a numerical mismatch. # # The severity scale is currently designed with the following integer codes: # 0 - the values are equal # 1 - the values are subject to suspect numerical noise (green fonts) # 2 - the values are different but below error threshold (blue fonts) # 3 - the values differ more than error threshold (red fonts) # # \param str_f_values: `array(string)` The strings representing the numeric # values read from the FORTRAN output file. # \param str_c_values: `array(string)` The strings representing the numeric # values read from the C++ output file. # \param config: `dict` A dictionary containing the configuration options from # which to read the warning and the error threshold. def mismatch_severities(str_f_values, str_c_values, config): """Determine the severity of a numerical mismatch. The severiti scale is currently designed with the following integer codes: 0 - the values are equal 1 - the values are subject to suspect numerical noise (green fonts) 2 - the values are different but below error threshold (blue fonts) 3 - the values differ more than error threshold (red fonts) ----------- Parameters: str_f_values: `array(string)` The strings representing the numeric values read from the FORTRAN output file. str_c_values: `array(string)` The strings representing the numeric values read from the C++ output file. config: `dict` A dictionary containing the configuration options from which to read the warning and the error threshold. """ result = [0 for ri in range(len(str_f_values))] for i in range(len(str_f_values)): if (str_f_values[i] != str_c_values[i]): Loading Loading @@ -204,6 +226,7 @@ def mismatch_severities(str_f_values, str_c_values, config): else: result[i] = 3 return result ## \brief Parse the command line arguments. def parse_arguments(): config = { 'fortran_file_name': '', Loading Loading @@ -234,13 +257,14 @@ def parse_arguments(): raise Exception("Unrecognized argument \'{0:s}\'".format(arg)) return config ## \brief Print a command-line help summary. def print_help(): print(" ") print("*** PYCOMPARE ***") print(" ") print("Compare the output of C++ and FORTRAN codes.") print(" ") print("Usage: \"./pycompare OPTIONS\" ") print("Usage: \"./pycompare.py OPTIONS\" ") print(" ") print("Valid options are: ") print("--ffile=FORTRAN_OUTPUT File containing the output of the FORTRAN code (mandatory).") Loading @@ -253,7 +277,9 @@ def print_help(): print(" ") ### PROGRAM EXECUTION ### # ### PROGRAM EXECUTION ### ## \cond res = main() ## \endcond if (res > 0): exit(1) exit(0) src/sphere/edfb.cppdeleted 100644 → 0 +0 −495 Original line number Diff line number Diff line /*! \file edfb.cpp */ #include <cstdio> #include <cmath> #include <complex> #include <cstring> #include <iostream> #include <fstream> #include "../include/file_io.h" #include "../include/List.h" using namespace std; /*! \brief Load a text file as a sequence of strings in memory. * * The configuration of the field expansion code in FORTRAN uses * shared memory access and file I/O operations managed by different * functions. Although this approach could be theoretically replicated, * it is more convenient to handle input and output to distinct files * using specific functions. load_file() helps in the task of handling * input such as configuration files or text data structures that need * to be loaded entirely. The function performs a line-by line scan of * the input file and returns an array of strings that can be later * parsed and ingested by the concerned code blocks. An optional pointer * to integer allows the function to keep track of the number of file * lines that were read, if needed. * * \param file_name: `string` The path of the file to be read. * \param count: `int*` Pointer to an integer recording the number of * read lines [OPTIONAL, default=NULL]. * \return array_lines `string*` An array of strings, one for each input * file line. */ string *load_file(string file_name, int *count); /*! \brief C++ implementation of EDFB * * This code aims at replicating the original work-flow in C++. */ int main(int argc, char **argv) { // Common variables set complex<double> *dc0, ***dc0m; double *ros, **rcf; int *iog, *nshl; double *xiv, *wns, *wls, *pus, *evs, *vss; string vns[5]; int ici; // Input file reading section int num_lines = 0; int last_read_line = 0; // Keep track of where the input stream was left string *file_lines = load_file("../../test_data/sphere/DEDFB", &num_lines); // Configuration code int nsph, ies; sscanf(file_lines[last_read_line].c_str(), " %d %d", &nsph, &ies); if (ies != 0) ies = 1; double exdc, wp, xip; int exdc_exp, wp_exp, xip_exp; int idfc, nxi, instpc, insn; int nsh; sscanf( file_lines[++last_read_line].c_str(), " %9lf D%d %9lf D%d %8lf D%d %d %d %d %d", &exdc, &exdc_exp, &wp, &wp_exp, &xip, &xip_exp, &idfc, &nxi, &instpc, &insn ); exdc *= pow(10.0, exdc_exp); wp *= pow(10.0, wp_exp); xip *= pow(10.0, xip_exp); FILE *output = fopen("c_OEDFB", "w"); // FORTRAN starts subroutine INXI at this point const double pigt = acos(0.0) * 4.0; const double evc = 6.5821188e-16; if (idfc >= 0) { // Not walked by default input data // This part of the code in not tested vss = new double[nxi]; xiv = new double[nxi]; pus = new double[nxi]; evs = new double[nxi]; wns = new double[nxi]; wls = new double[nxi]; if (instpc == 0) { // The variable vector is explicitly defined double vs; int vs_exp; for (int jxi_r = 0; jxi_r < nxi; jxi_r++) { sscanf(file_lines[++last_read_line].c_str(), " %lf D%d", &vs, &vs_exp); vs *= pow(10.0, vs_exp); vss[jxi_r] = vs; } switch (insn) { case 1: //xi vector definition vns[insn - 1] = "XIV"; fprintf(output, " JXI XIV WNS WLS PUS EVS\n"); for (int jxi210w = 0; jxi210w < nxi; jxi210w++) { xiv[jxi210w] = vss[jxi210w]; pus[jxi210w] = xiv[jxi210w] * wp; evs[jxi210w] = pus[jxi210w] * evc; wns[jxi210w] = pus[jxi210w] / 3.0e8; wls[jxi210w] = pigt / wns[jxi210w]; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi210w + 1), xiv[jxi210w], wns[jxi210w], wls[jxi210w], pus[jxi210w], evs[jxi210w] ); } break; case 2: //wave number vector definition vns[insn - 1] = "WNS"; fprintf(output, " JXI WNS WLS PUS EVS XIV\n"); for (int jxi230w = 0; jxi230w < nxi; jxi230w++) { wns[jxi230w] = vss[jxi230w]; wls[jxi230w] = pigt / wns[jxi230w]; xiv[jxi230w] = 3.0e8 * wns[jxi230w] / wp; pus[jxi230w] = xiv[jxi230w] * wp; evs[jxi230w] = pus[jxi230w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi230w + 1), wns[jxi230w], wls[jxi230w], pus[jxi230w], evs[jxi230w], xiv[jxi230w] ); } break; case 3: //wavelength vector definition vns[insn - 1] = "WLS"; fprintf(output, " JXI WLS WNS PUS EVS XIV\n"); for (int jxi250w = 0; jxi250w < nxi; jxi250w++) { wls[jxi250w] = vss[jxi250w]; wns[jxi250w] = pigt / wls[jxi250w]; xiv[jxi250w] = 3.0e8 * wns[jxi250w] / wp; pus[jxi250w] = xiv[jxi250w] * wp; evs[jxi250w] = pus[jxi250w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi250w + 1), wls[jxi250w], wns[jxi250w], pus[jxi250w], evs[jxi250w], xiv[jxi250w] ); } break; case 4: //pu vector definition vns[insn - 1] = "PUS"; fprintf(output, " JXI PUS WNS WLS EVS XIV\n"); for (int jxi270w = 0; jxi270w < nxi; jxi270w++) { pus[jxi270w] = vss[jxi270w]; xiv[jxi270w] = pus[jxi270w] / wp; wns[jxi270w] = pus[jxi270w] / 3.0e8; wls[jxi270w] = pigt / wns[jxi270w]; evs[jxi270w] = pus[jxi270w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi270w + 1), pus[jxi270w], wns[jxi270w], wls[jxi270w], evs[jxi270w], xiv[jxi270w] ); } break; case 5: //eV vector definition vns[insn - 1] = "EVS"; fprintf(output, " JXI EVS WNS WLS PUS XIV\n"); for (int jxi290w = 0; jxi290w < nxi; jxi290w++) { evs[jxi290w] = vss[jxi290w]; pus[jxi290w] = evs[jxi290w] / evc; xiv[jxi290w] = pus[jxi290w] / wp; wns[jxi290w] = pus[jxi290w] / 3.0e8; wls[jxi290w] = pigt / wns[jxi290w]; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi290w + 1), evs[jxi290w], wns[jxi290w], wls[jxi290w], pus[jxi290w], xiv[jxi290w] ); } break; } } else { // The variable vector needs to be computed in steps double vs, vs_step; int vs_exp, vs_step_exp; sscanf(file_lines[++last_read_line].c_str(), " %lf D%d %lf D%d", &vs, &vs_exp, &vs_step, &vs_step_exp); vs *= pow(10.0, vs_exp); vs_step *= pow(10.0, vs_step_exp); switch (insn) { case 1: //xi vector definition vns[insn - 1] = "XIV"; fprintf(output, " JXI XIV WNS WLS PUS EVS\n"); for (int jxi110w = 0; jxi110w < nxi; jxi110w++) { vss[jxi110w] = vs; xiv[jxi110w] = vss[jxi110w]; pus[jxi110w] = xiv[jxi110w] * wp; wns[jxi110w] = pus[jxi110w] / 3.0e8; evs[jxi110w] = pus[jxi110w] * evc; wls[jxi110w] = pigt / wns[jxi110w]; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi110w + 1), xiv[jxi110w], wns[jxi110w], wls[jxi110w], pus[jxi110w], evs[jxi110w] ); vs += vs_step; } break; case 2: //wave number vector definition vns[insn - 1] = "WNS"; fprintf(output, " JXI WNS WLS PUS EVS XIV\n"); for (int jxi130w = 0; jxi130w < nxi; jxi130w++) { vss[jxi130w] = vs; wns[jxi130w] = vss[jxi130w]; xiv[jxi130w] = 3.0e8 * wns[jxi130w] / wp; pus[jxi130w] = xiv[jxi130w] * wp; wls[jxi130w] = pigt / wns[jxi130w]; evs[jxi130w] = pus[jxi130w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi130w + 1), wns[jxi130w], wls[jxi130w], pus[jxi130w], evs[jxi130w], xiv[jxi130w] ); vs += vs_step; } break; case 3: //wavelength vector definition vns[insn - 1] = "WLS"; fprintf(output, " JXI WLS WNS PUS EVS XIV\n"); for (int jxi150w = 0; jxi150w < nxi; jxi150w++) { vss[jxi150w] = vs; wls[jxi150w] = vss[jxi150w]; wns[jxi150w] = pigt / wls[jxi150w]; xiv[jxi150w] = 3.0e8 * wns[jxi150w] / wp; pus[jxi150w] = xiv[jxi150w] * wp; evs[jxi150w] = pus[jxi150w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi150w + 1), wls[jxi150w], wns[jxi150w], pus[jxi150w], evs[jxi150w], xiv[jxi150w] ); vs += vs_step; } break; case 4: //pu vector definition vns[insn - 1] = "PUS"; fprintf(output, " JXI PUS WNS WLS EVS XIV\n"); for (int jxi170w = 0; jxi170w < nxi; jxi170w++) { vss[jxi170w] = vs; pus[jxi170w] = vss[jxi170w]; xiv[jxi170w] = pus[jxi170w] / wp; wns[jxi170w] = pus[jxi170w] / 3.0e8; wls[jxi170w] = pigt / wns[jxi170w]; evs[jxi170w] = pus[jxi170w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi170w + 1), pus[jxi170w], wns[jxi170w], wls[jxi170w], evs[jxi170w], xiv[jxi170w] ); vs += vs_step; } break; case 5: //eV vector definition vns[insn - 1] = "EVS"; fprintf(output, " JXI EVS WNS WLS PUS XIV\n"); for (int jxi190w = 0; jxi190w < nxi; jxi190w++) { vss[jxi190w] = vs; evs[jxi190w] = vss[jxi190w]; pus[jxi190w] = evs[jxi190w] / evc; xiv[jxi190w] = pus[jxi190w] / wp; wns[jxi190w] = pus[jxi190w] / 3.0e8; wls[jxi190w] = pigt / wns[jxi190w]; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi190w + 1), evs[jxi190w], wns[jxi190w], wls[jxi190w], pus[jxi190w], xiv[jxi190w] ); vs += vs_step; } break; } } // End of the untested code section. } else { if (instpc < 1) { // In this case the XI vector is explicitly defined. // Test input comes this way. double xi, pu, wn; int xi_exp; vns[insn - 1] = "XIV"; List<double> xi_vector; sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp); xi *= pow(10.0, xi_exp); xi_vector.set(0, xi); for (int jxi310 = 1; jxi310 < nxi; jxi310++) { sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp); xi *= pow(10.0, xi_exp); xi_vector.append(xi); } vss = xi_vector.to_array(); xiv = xi_vector.to_array(); pu = xip * wp; wn = pu / 3.0e8; fprintf(output, " XIP WN WL PU EV\n"); fprintf(output, " %13.4lE", xip); fprintf(output, "%13.4lE", wn); fprintf(output, "%13.4lE", pigt / wn); fprintf(output, "%13.4lE", pu); fprintf(output, "%13.4lE\n", pu * evc); fprintf(output, " SCALE FACTORS XI\n"); for (int jxi6612 = 1; jxi6612 <= nxi; jxi6612++) fprintf(output, "%5d%13.4lE\n", jxi6612, xiv[jxi6612 - 1]); //INXI branch ends here. } } last_read_line++; iog = new int[nsph]; for (int i = 0; i < nsph; i++) { string read_format = ""; for (int j = 0; j < i; j++) read_format += " %*d"; read_format += " %d"; sscanf(file_lines[last_read_line].c_str(), read_format.c_str(), (iog + i)); } nshl = new int[nsph]; ros = new double[nsph]; rcf = new double*[nsph]; for (int i113 = 1; i113 <= nsph; i113++) { int i_val; double ros_val; int ros_val_exp; if (iog[i113 - 1] < i113) continue; sscanf(file_lines[++last_read_line].c_str(), " %d %9lf D%d", &i_val, &ros_val, &ros_val_exp); nshl[i113 - 1] = i_val; ros[i113 - 1] = ros_val * pow(10.0, ros_val_exp); nsh = nshl[i113 -1]; if (i113 == 1) nsh += ies; rcf[i113 - 1] = new double[nsh]; for (int ns = 0; ns < nsh; ns++) { double ns_rcf; int ns_rcf_exp; sscanf(file_lines[++last_read_line].c_str(), " %8lf D%d", &ns_rcf, &ns_rcf_exp); rcf[i113 -1][ns] = ns_rcf * pow(10.0, ns_rcf_exp); } } // The FORTRAN code writes an auxiliary file in binary format. This should // be avoided or possibly replaced with the use of standard file formats for // scientific use (e.g. FITS). int uid = 27; string bin_file_name = "c_TEDF"; string status = "UNKNOWN"; string mode = "UNFORMATTED"; open_file_(&uid, bin_file_name.c_str(), status.c_str(), mode.c_str()); write_int_(&uid, &nsph); for (int iogi = 0; iogi < nsph; iogi++) write_int_(&uid, (iog + iogi)); write_double_(&uid, &exdc); write_double_(&uid, &wp); write_double_(&uid, &xip); write_int_(&uid, &idfc); write_int_(&uid, &nxi); for (int xivi = 0; xivi < nxi; xivi++) write_double_(&uid, (xiv + xivi)); for (int i115 = 1; i115 <= nsph; i115++) { if (iog[i115 - 1] < i115) continue; write_int_(&uid, (nshl + i115 -1)); write_double_(&uid, (ros + i115 - 1)); nsh = nshl[i115 - 1]; if (i115 == 1) nsh += ies; for (int ins = 0; ins < nsh; ins++) write_double_(&uid, (rcf[i115 - 1] + ins)); } // Remake the dc0m matrix. dc0m = new complex<double>**[nsph]; for (int dim1 = 0; dim1 < nsph; dim1++) { dc0m[dim1] = new complex<double>*[nsph]; for (int dim2 = 0; dim2 < nxi; dim2++) { dc0m[dim1][dim2] = new complex<double>[nxi]; } } for (int jxi468 = 1; jxi468 <= nxi; jxi468++) { if (idfc != 0 && jxi468 > 1) continue; for (int i162 = 1; i162 <= nsph; i162++) { if (iog[i162 - 1] < i162) continue; nsh = nshl[i162 - 1]; ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here? if (i162 == 1) ici = ici + ies; for (int i157 = 0; i157 < ici; i157++) { double dc0_real, dc0_img; int dc0_real_exp, dc0_img_exp; sscanf(file_lines[++last_read_line].c_str(), " (%8lf D%d, %8lf D%d)", &dc0_real, &dc0_real_exp, &dc0_img, &dc0_img_exp); dc0_real *= pow(10.0, dc0_real_exp); dc0_img *= pow(10.0, dc0_img_exp); dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img; // The FORTRAN code writes the complex numbers as a 16-byte long binary stream. // Here we assume that the 16 bytes are equally split in 8 bytes to represent the // real part and 8 bytes to represent the imaginary one. write_complex_(&uid, &dc0_real, &dc0_img); } } } close_file_(&uid); if (idfc != 0) { fprintf(output, " DIELECTRIC CONSTANTS\n"); for (int i473 = 1; i473 <= nsph; i473++) { if (iog[i473 - 1] != i473) continue; ici = (nshl[i473 - 1] + 1) / 2; if (i473 == 1) ici += ies; fprintf(output, " SPHERE N. %4d\n", i473); for (int ic472 = 0; ic472 < ici; ic472++) { double dc0_real = dc0m[ic472][i473 - 1][0].real(), dc0_img = dc0m[ic472][i473 - 1][0].imag(); fprintf(output, "%5d %12.4lE%12.4lE\n", (ic472 + 1), dc0_real, dc0_img); } } } else { fprintf(output, " DIELECTRIC FUNCTIONS\n"); for (int i478 = 1; i478 <= nsph; i478++) { if (iog[i478 - 1] != i478) continue; ici = (nshl[i478 - 1] + 1) / 2; if (i478 == 1) ici += ies; fprintf(output, " SPHERE N. %4d\n", i478); for (int ic477 = 1; ic477 <= ici; ic477++) { fprintf(output, " NONTRANSITION LAYER N. %2d , SCALE = %3c\n", ic477, vns[insn - 1].c_str()); for (int jxi476 = 0; jxi476 < nxi; jxi476++) { double dc0_real = dc0m[ic477 - 1][i478 - 1][jxi476].real(); double dc0_img = dc0m[ic477 - 1][i478 - 1][jxi476].imag(); fprintf(output, "%5d (%12.4lE,%12.4lE)\n", (jxi476 + 1), dc0_real, dc0_img); } } } } fclose(output); return 0; } string *load_file(string file_name, int *count = 0) { fstream input_file(file_name.c_str(), ios::in); List<string> file_lines = List<string>(); string line; if (input_file.is_open()) { getline(input_file, line); file_lines.set(0, line); while (getline(input_file, line)) { file_lines.append(line); } input_file.close(); } string *array_lines = file_lines.to_array(); if (count != 0) *count = file_lines.length(); return array_lines; } Loading
doc/src/config.dox +2 −1 Original line number Diff line number Diff line Loading @@ -986,7 +986,8 @@ INPUT_FILE_ENCODING = FILE_PATTERNS = *.cpp \ *.f \ *.h \ *.md *.md \ *.py # The RECURSIVE tag can be used to specify whether or not subdirectories should # be searched for input files as well. Loading
src/scripts/pycompare→src/scripts/pycompare.py +47 −21 Original line number Diff line number Diff line #!/bin/python ## @package pycompare # Script to perform output consistency tests # # Comparing the numeric output can be rendered hard by the amount of information # contained in a typical output file and the necessity to determine whether a # difference is actually significant or just caused by numeric noise hitting # negligible values. The task of this script is to compare two output files, in # the assumption that they were written by the FORTRAN and the C++ versions of # the code and to flag all the possible inconsistencies according to various # severity levels (namely: NOISE, WARNING, and ERROR). import re from math import log10 from sys import argv ## \cond number_reg = re.compile(r'-?[0-9]\.[0-9]+E[-+][0-9]{2,2}') ## \endcond ## \brief Main execution code # # `main()` is the function that handles the creation of the script configuration # and the execution of the comparison. It returns an integer value corresponding # to the number of detected error-level inconsistencies. # # \returns errors: `int` Number of detected error-level inconsistencies. def main(): config = parse_arguments() errors, warnings, noisy = (0, 0, 0) Loading @@ -32,6 +52,7 @@ def main(): )) return errors ## \brief Perform the comparison of two files. def compare_files(config): mismatch_count = { 'errors': 0, Loading Loading @@ -78,8 +99,14 @@ def compare_files(config): l_file.write(" </body>\n") l_file.write("</html>\n") l_file.close() else: mismatch_count['errors'] = len(c_lines) print("ERROR: {0:s} and {1:s} have different numbers of lines!".format( config['fortran_file_name'], config['c_file_name'] )) return mismatch_count ## \brief Perform the comparison of two file lines. def compare_lines(f_line, c_line, config, line_num=0, num_len=1, log_file=None): errors = 0 warnings = 0 Loading Loading @@ -156,26 +183,21 @@ def compare_lines(f_line, c_line, config, line_num=0, num_len=1, log_file=None): log_file.write(log_line + "</code></pre></div>\n") return (errors, warnings, noisy) ## \brief Determine the severity of a numerical mismatch. # # The severity scale is currently designed with the following integer codes: # 0 - the values are equal # 1 - the values are subject to suspect numerical noise (green fonts) # 2 - the values are different but below error threshold (blue fonts) # 3 - the values differ more than error threshold (red fonts) # # \param str_f_values: `array(string)` The strings representing the numeric # values read from the FORTRAN output file. # \param str_c_values: `array(string)` The strings representing the numeric # values read from the C++ output file. # \param config: `dict` A dictionary containing the configuration options from # which to read the warning and the error threshold. def mismatch_severities(str_f_values, str_c_values, config): """Determine the severity of a numerical mismatch. The severiti scale is currently designed with the following integer codes: 0 - the values are equal 1 - the values are subject to suspect numerical noise (green fonts) 2 - the values are different but below error threshold (blue fonts) 3 - the values differ more than error threshold (red fonts) ----------- Parameters: str_f_values: `array(string)` The strings representing the numeric values read from the FORTRAN output file. str_c_values: `array(string)` The strings representing the numeric values read from the C++ output file. config: `dict` A dictionary containing the configuration options from which to read the warning and the error threshold. """ result = [0 for ri in range(len(str_f_values))] for i in range(len(str_f_values)): if (str_f_values[i] != str_c_values[i]): Loading Loading @@ -204,6 +226,7 @@ def mismatch_severities(str_f_values, str_c_values, config): else: result[i] = 3 return result ## \brief Parse the command line arguments. def parse_arguments(): config = { 'fortran_file_name': '', Loading Loading @@ -234,13 +257,14 @@ def parse_arguments(): raise Exception("Unrecognized argument \'{0:s}\'".format(arg)) return config ## \brief Print a command-line help summary. def print_help(): print(" ") print("*** PYCOMPARE ***") print(" ") print("Compare the output of C++ and FORTRAN codes.") print(" ") print("Usage: \"./pycompare OPTIONS\" ") print("Usage: \"./pycompare.py OPTIONS\" ") print(" ") print("Valid options are: ") print("--ffile=FORTRAN_OUTPUT File containing the output of the FORTRAN code (mandatory).") Loading @@ -253,7 +277,9 @@ def print_help(): print(" ") ### PROGRAM EXECUTION ### # ### PROGRAM EXECUTION ### ## \cond res = main() ## \endcond if (res > 0): exit(1) exit(0)
src/sphere/edfb.cppdeleted 100644 → 0 +0 −495 Original line number Diff line number Diff line /*! \file edfb.cpp */ #include <cstdio> #include <cmath> #include <complex> #include <cstring> #include <iostream> #include <fstream> #include "../include/file_io.h" #include "../include/List.h" using namespace std; /*! \brief Load a text file as a sequence of strings in memory. * * The configuration of the field expansion code in FORTRAN uses * shared memory access and file I/O operations managed by different * functions. Although this approach could be theoretically replicated, * it is more convenient to handle input and output to distinct files * using specific functions. load_file() helps in the task of handling * input such as configuration files or text data structures that need * to be loaded entirely. The function performs a line-by line scan of * the input file and returns an array of strings that can be later * parsed and ingested by the concerned code blocks. An optional pointer * to integer allows the function to keep track of the number of file * lines that were read, if needed. * * \param file_name: `string` The path of the file to be read. * \param count: `int*` Pointer to an integer recording the number of * read lines [OPTIONAL, default=NULL]. * \return array_lines `string*` An array of strings, one for each input * file line. */ string *load_file(string file_name, int *count); /*! \brief C++ implementation of EDFB * * This code aims at replicating the original work-flow in C++. */ int main(int argc, char **argv) { // Common variables set complex<double> *dc0, ***dc0m; double *ros, **rcf; int *iog, *nshl; double *xiv, *wns, *wls, *pus, *evs, *vss; string vns[5]; int ici; // Input file reading section int num_lines = 0; int last_read_line = 0; // Keep track of where the input stream was left string *file_lines = load_file("../../test_data/sphere/DEDFB", &num_lines); // Configuration code int nsph, ies; sscanf(file_lines[last_read_line].c_str(), " %d %d", &nsph, &ies); if (ies != 0) ies = 1; double exdc, wp, xip; int exdc_exp, wp_exp, xip_exp; int idfc, nxi, instpc, insn; int nsh; sscanf( file_lines[++last_read_line].c_str(), " %9lf D%d %9lf D%d %8lf D%d %d %d %d %d", &exdc, &exdc_exp, &wp, &wp_exp, &xip, &xip_exp, &idfc, &nxi, &instpc, &insn ); exdc *= pow(10.0, exdc_exp); wp *= pow(10.0, wp_exp); xip *= pow(10.0, xip_exp); FILE *output = fopen("c_OEDFB", "w"); // FORTRAN starts subroutine INXI at this point const double pigt = acos(0.0) * 4.0; const double evc = 6.5821188e-16; if (idfc >= 0) { // Not walked by default input data // This part of the code in not tested vss = new double[nxi]; xiv = new double[nxi]; pus = new double[nxi]; evs = new double[nxi]; wns = new double[nxi]; wls = new double[nxi]; if (instpc == 0) { // The variable vector is explicitly defined double vs; int vs_exp; for (int jxi_r = 0; jxi_r < nxi; jxi_r++) { sscanf(file_lines[++last_read_line].c_str(), " %lf D%d", &vs, &vs_exp); vs *= pow(10.0, vs_exp); vss[jxi_r] = vs; } switch (insn) { case 1: //xi vector definition vns[insn - 1] = "XIV"; fprintf(output, " JXI XIV WNS WLS PUS EVS\n"); for (int jxi210w = 0; jxi210w < nxi; jxi210w++) { xiv[jxi210w] = vss[jxi210w]; pus[jxi210w] = xiv[jxi210w] * wp; evs[jxi210w] = pus[jxi210w] * evc; wns[jxi210w] = pus[jxi210w] / 3.0e8; wls[jxi210w] = pigt / wns[jxi210w]; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi210w + 1), xiv[jxi210w], wns[jxi210w], wls[jxi210w], pus[jxi210w], evs[jxi210w] ); } break; case 2: //wave number vector definition vns[insn - 1] = "WNS"; fprintf(output, " JXI WNS WLS PUS EVS XIV\n"); for (int jxi230w = 0; jxi230w < nxi; jxi230w++) { wns[jxi230w] = vss[jxi230w]; wls[jxi230w] = pigt / wns[jxi230w]; xiv[jxi230w] = 3.0e8 * wns[jxi230w] / wp; pus[jxi230w] = xiv[jxi230w] * wp; evs[jxi230w] = pus[jxi230w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi230w + 1), wns[jxi230w], wls[jxi230w], pus[jxi230w], evs[jxi230w], xiv[jxi230w] ); } break; case 3: //wavelength vector definition vns[insn - 1] = "WLS"; fprintf(output, " JXI WLS WNS PUS EVS XIV\n"); for (int jxi250w = 0; jxi250w < nxi; jxi250w++) { wls[jxi250w] = vss[jxi250w]; wns[jxi250w] = pigt / wls[jxi250w]; xiv[jxi250w] = 3.0e8 * wns[jxi250w] / wp; pus[jxi250w] = xiv[jxi250w] * wp; evs[jxi250w] = pus[jxi250w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi250w + 1), wls[jxi250w], wns[jxi250w], pus[jxi250w], evs[jxi250w], xiv[jxi250w] ); } break; case 4: //pu vector definition vns[insn - 1] = "PUS"; fprintf(output, " JXI PUS WNS WLS EVS XIV\n"); for (int jxi270w = 0; jxi270w < nxi; jxi270w++) { pus[jxi270w] = vss[jxi270w]; xiv[jxi270w] = pus[jxi270w] / wp; wns[jxi270w] = pus[jxi270w] / 3.0e8; wls[jxi270w] = pigt / wns[jxi270w]; evs[jxi270w] = pus[jxi270w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi270w + 1), pus[jxi270w], wns[jxi270w], wls[jxi270w], evs[jxi270w], xiv[jxi270w] ); } break; case 5: //eV vector definition vns[insn - 1] = "EVS"; fprintf(output, " JXI EVS WNS WLS PUS XIV\n"); for (int jxi290w = 0; jxi290w < nxi; jxi290w++) { evs[jxi290w] = vss[jxi290w]; pus[jxi290w] = evs[jxi290w] / evc; xiv[jxi290w] = pus[jxi290w] / wp; wns[jxi290w] = pus[jxi290w] / 3.0e8; wls[jxi290w] = pigt / wns[jxi290w]; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi290w + 1), evs[jxi290w], wns[jxi290w], wls[jxi290w], pus[jxi290w], xiv[jxi290w] ); } break; } } else { // The variable vector needs to be computed in steps double vs, vs_step; int vs_exp, vs_step_exp; sscanf(file_lines[++last_read_line].c_str(), " %lf D%d %lf D%d", &vs, &vs_exp, &vs_step, &vs_step_exp); vs *= pow(10.0, vs_exp); vs_step *= pow(10.0, vs_step_exp); switch (insn) { case 1: //xi vector definition vns[insn - 1] = "XIV"; fprintf(output, " JXI XIV WNS WLS PUS EVS\n"); for (int jxi110w = 0; jxi110w < nxi; jxi110w++) { vss[jxi110w] = vs; xiv[jxi110w] = vss[jxi110w]; pus[jxi110w] = xiv[jxi110w] * wp; wns[jxi110w] = pus[jxi110w] / 3.0e8; evs[jxi110w] = pus[jxi110w] * evc; wls[jxi110w] = pigt / wns[jxi110w]; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi110w + 1), xiv[jxi110w], wns[jxi110w], wls[jxi110w], pus[jxi110w], evs[jxi110w] ); vs += vs_step; } break; case 2: //wave number vector definition vns[insn - 1] = "WNS"; fprintf(output, " JXI WNS WLS PUS EVS XIV\n"); for (int jxi130w = 0; jxi130w < nxi; jxi130w++) { vss[jxi130w] = vs; wns[jxi130w] = vss[jxi130w]; xiv[jxi130w] = 3.0e8 * wns[jxi130w] / wp; pus[jxi130w] = xiv[jxi130w] * wp; wls[jxi130w] = pigt / wns[jxi130w]; evs[jxi130w] = pus[jxi130w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi130w + 1), wns[jxi130w], wls[jxi130w], pus[jxi130w], evs[jxi130w], xiv[jxi130w] ); vs += vs_step; } break; case 3: //wavelength vector definition vns[insn - 1] = "WLS"; fprintf(output, " JXI WLS WNS PUS EVS XIV\n"); for (int jxi150w = 0; jxi150w < nxi; jxi150w++) { vss[jxi150w] = vs; wls[jxi150w] = vss[jxi150w]; wns[jxi150w] = pigt / wls[jxi150w]; xiv[jxi150w] = 3.0e8 * wns[jxi150w] / wp; pus[jxi150w] = xiv[jxi150w] * wp; evs[jxi150w] = pus[jxi150w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi150w + 1), wls[jxi150w], wns[jxi150w], pus[jxi150w], evs[jxi150w], xiv[jxi150w] ); vs += vs_step; } break; case 4: //pu vector definition vns[insn - 1] = "PUS"; fprintf(output, " JXI PUS WNS WLS EVS XIV\n"); for (int jxi170w = 0; jxi170w < nxi; jxi170w++) { vss[jxi170w] = vs; pus[jxi170w] = vss[jxi170w]; xiv[jxi170w] = pus[jxi170w] / wp; wns[jxi170w] = pus[jxi170w] / 3.0e8; wls[jxi170w] = pigt / wns[jxi170w]; evs[jxi170w] = pus[jxi170w] * evc; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi170w + 1), pus[jxi170w], wns[jxi170w], wls[jxi170w], evs[jxi170w], xiv[jxi170w] ); vs += vs_step; } break; case 5: //eV vector definition vns[insn - 1] = "EVS"; fprintf(output, " JXI EVS WNS WLS PUS XIV\n"); for (int jxi190w = 0; jxi190w < nxi; jxi190w++) { vss[jxi190w] = vs; evs[jxi190w] = vss[jxi190w]; pus[jxi190w] = evs[jxi190w] / evc; xiv[jxi190w] = pus[jxi190w] / wp; wns[jxi190w] = pus[jxi190w] / 3.0e8; wls[jxi190w] = pigt / wns[jxi190w]; fprintf( output, "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n", (jxi190w + 1), evs[jxi190w], wns[jxi190w], wls[jxi190w], pus[jxi190w], xiv[jxi190w] ); vs += vs_step; } break; } } // End of the untested code section. } else { if (instpc < 1) { // In this case the XI vector is explicitly defined. // Test input comes this way. double xi, pu, wn; int xi_exp; vns[insn - 1] = "XIV"; List<double> xi_vector; sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp); xi *= pow(10.0, xi_exp); xi_vector.set(0, xi); for (int jxi310 = 1; jxi310 < nxi; jxi310++) { sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp); xi *= pow(10.0, xi_exp); xi_vector.append(xi); } vss = xi_vector.to_array(); xiv = xi_vector.to_array(); pu = xip * wp; wn = pu / 3.0e8; fprintf(output, " XIP WN WL PU EV\n"); fprintf(output, " %13.4lE", xip); fprintf(output, "%13.4lE", wn); fprintf(output, "%13.4lE", pigt / wn); fprintf(output, "%13.4lE", pu); fprintf(output, "%13.4lE\n", pu * evc); fprintf(output, " SCALE FACTORS XI\n"); for (int jxi6612 = 1; jxi6612 <= nxi; jxi6612++) fprintf(output, "%5d%13.4lE\n", jxi6612, xiv[jxi6612 - 1]); //INXI branch ends here. } } last_read_line++; iog = new int[nsph]; for (int i = 0; i < nsph; i++) { string read_format = ""; for (int j = 0; j < i; j++) read_format += " %*d"; read_format += " %d"; sscanf(file_lines[last_read_line].c_str(), read_format.c_str(), (iog + i)); } nshl = new int[nsph]; ros = new double[nsph]; rcf = new double*[nsph]; for (int i113 = 1; i113 <= nsph; i113++) { int i_val; double ros_val; int ros_val_exp; if (iog[i113 - 1] < i113) continue; sscanf(file_lines[++last_read_line].c_str(), " %d %9lf D%d", &i_val, &ros_val, &ros_val_exp); nshl[i113 - 1] = i_val; ros[i113 - 1] = ros_val * pow(10.0, ros_val_exp); nsh = nshl[i113 -1]; if (i113 == 1) nsh += ies; rcf[i113 - 1] = new double[nsh]; for (int ns = 0; ns < nsh; ns++) { double ns_rcf; int ns_rcf_exp; sscanf(file_lines[++last_read_line].c_str(), " %8lf D%d", &ns_rcf, &ns_rcf_exp); rcf[i113 -1][ns] = ns_rcf * pow(10.0, ns_rcf_exp); } } // The FORTRAN code writes an auxiliary file in binary format. This should // be avoided or possibly replaced with the use of standard file formats for // scientific use (e.g. FITS). int uid = 27; string bin_file_name = "c_TEDF"; string status = "UNKNOWN"; string mode = "UNFORMATTED"; open_file_(&uid, bin_file_name.c_str(), status.c_str(), mode.c_str()); write_int_(&uid, &nsph); for (int iogi = 0; iogi < nsph; iogi++) write_int_(&uid, (iog + iogi)); write_double_(&uid, &exdc); write_double_(&uid, &wp); write_double_(&uid, &xip); write_int_(&uid, &idfc); write_int_(&uid, &nxi); for (int xivi = 0; xivi < nxi; xivi++) write_double_(&uid, (xiv + xivi)); for (int i115 = 1; i115 <= nsph; i115++) { if (iog[i115 - 1] < i115) continue; write_int_(&uid, (nshl + i115 -1)); write_double_(&uid, (ros + i115 - 1)); nsh = nshl[i115 - 1]; if (i115 == 1) nsh += ies; for (int ins = 0; ins < nsh; ins++) write_double_(&uid, (rcf[i115 - 1] + ins)); } // Remake the dc0m matrix. dc0m = new complex<double>**[nsph]; for (int dim1 = 0; dim1 < nsph; dim1++) { dc0m[dim1] = new complex<double>*[nsph]; for (int dim2 = 0; dim2 < nxi; dim2++) { dc0m[dim1][dim2] = new complex<double>[nxi]; } } for (int jxi468 = 1; jxi468 <= nxi; jxi468++) { if (idfc != 0 && jxi468 > 1) continue; for (int i162 = 1; i162 <= nsph; i162++) { if (iog[i162 - 1] < i162) continue; nsh = nshl[i162 - 1]; ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here? if (i162 == 1) ici = ici + ies; for (int i157 = 0; i157 < ici; i157++) { double dc0_real, dc0_img; int dc0_real_exp, dc0_img_exp; sscanf(file_lines[++last_read_line].c_str(), " (%8lf D%d, %8lf D%d)", &dc0_real, &dc0_real_exp, &dc0_img, &dc0_img_exp); dc0_real *= pow(10.0, dc0_real_exp); dc0_img *= pow(10.0, dc0_img_exp); dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img; // The FORTRAN code writes the complex numbers as a 16-byte long binary stream. // Here we assume that the 16 bytes are equally split in 8 bytes to represent the // real part and 8 bytes to represent the imaginary one. write_complex_(&uid, &dc0_real, &dc0_img); } } } close_file_(&uid); if (idfc != 0) { fprintf(output, " DIELECTRIC CONSTANTS\n"); for (int i473 = 1; i473 <= nsph; i473++) { if (iog[i473 - 1] != i473) continue; ici = (nshl[i473 - 1] + 1) / 2; if (i473 == 1) ici += ies; fprintf(output, " SPHERE N. %4d\n", i473); for (int ic472 = 0; ic472 < ici; ic472++) { double dc0_real = dc0m[ic472][i473 - 1][0].real(), dc0_img = dc0m[ic472][i473 - 1][0].imag(); fprintf(output, "%5d %12.4lE%12.4lE\n", (ic472 + 1), dc0_real, dc0_img); } } } else { fprintf(output, " DIELECTRIC FUNCTIONS\n"); for (int i478 = 1; i478 <= nsph; i478++) { if (iog[i478 - 1] != i478) continue; ici = (nshl[i478 - 1] + 1) / 2; if (i478 == 1) ici += ies; fprintf(output, " SPHERE N. %4d\n", i478); for (int ic477 = 1; ic477 <= ici; ic477++) { fprintf(output, " NONTRANSITION LAYER N. %2d , SCALE = %3c\n", ic477, vns[insn - 1].c_str()); for (int jxi476 = 0; jxi476 < nxi; jxi476++) { double dc0_real = dc0m[ic477 - 1][i478 - 1][jxi476].real(); double dc0_img = dc0m[ic477 - 1][i478 - 1][jxi476].imag(); fprintf(output, "%5d (%12.4lE,%12.4lE)\n", (jxi476 + 1), dc0_real, dc0_img); } } } } fclose(output); return 0; } string *load_file(string file_name, int *count = 0) { fstream input_file(file_name.c_str(), ios::in); List<string> file_lines = List<string>(); string line; if (input_file.is_open()) { getline(input_file, line); file_lines.set(0, line); while (getline(input_file, line)) { file_lines.append(line); } input_file.close(); } string *array_lines = file_lines.to_array(); if (count != 0) *count = file_lines.length(); return array_lines; }