diff --git a/tests/03_NAO_multik/CASES_GPU.txt b/tests/03_NAO_multik/CASES_GPU.txt index 88d9f4478a..1453c8cf28 100644 --- a/tests/03_NAO_multik/CASES_GPU.txt +++ b/tests/03_NAO_multik/CASES_GPU.txt @@ -33,7 +33,7 @@ scf_out_hsk scf_out_hsk_binary scf_out_hsr scf_out_hsr_binary_spin2 -#scf_out_hsr_spin4 +scf_out_hsr_spin4 scf_out_dh_t scf_out_dos_spin4 scf_out_mul @@ -46,7 +46,7 @@ nscf_out_dos nscf_out_band_pband nscf_out_pot1 nscf_out_mul -#nscf_out_hsr_tr_rr +nscf_out_hsr_tr_rr relax_bfgs2 relax_old_cg relax_cell diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index d37daa7af6..bbe7dc6676 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -445,13 +445,14 @@ fi #----------------------------------- #echo $has_hs2 if ! test -z "$has_hs2" && [ $has_hs2 == 1 ]; then - python3 $COMPARE_SCRIPT hrs1_nao.csr.ref OUT.autotest/hrs1_nao.csr 8 + HSR_CSR_COMPARE="../../integrate/tools/compare_hsr_csr.py" + python3 $HSR_CSR_COMPARE hrs1_nao.csr.ref OUT.autotest/hrs1_nao.csr 8 echo "CompareHR_pass $?" >>$1 if ! test -z "$nspin" && [ "$nspin" -eq 2 ]; then - python3 $COMPARE_SCRIPT hrs2_nao.csr.ref OUT.autotest/hrs2_nao.csr 8 + python3 $HSR_CSR_COMPARE hrs2_nao.csr.ref OUT.autotest/hrs2_nao.csr 8 echo "CompareHR2_pass $?" >>$1 fi - python3 $COMPARE_SCRIPT sr_nao.csr.ref OUT.autotest/sr_nao.csr 8 + python3 $HSR_CSR_COMPARE sr_nao.csr.ref OUT.autotest/sr_nao.csr 8 echo "CompareSR_pass $?" >>$1 elif ! test -z "$has_hs2" && [ "$has_hs2" == 2 ]; then HSR_BINARY_COMPARE="../../integrate/tools/compare_hsr_binary.py" diff --git a/tests/integrate/tools/compare_hsr_csr.py b/tests/integrate/tools/compare_hsr_csr.py new file mode 100755 index 0000000000..6cc510d327 --- /dev/null +++ b/tests/integrate/tools/compare_hsr_csr.py @@ -0,0 +1,222 @@ +#!/usr/bin/env python3 + +import math +import sys + +from compare_hsr_binary import read_text_reference + + +def csr_to_entries(nbasis, values, columns, row_ptr, rvec): + if len(columns) != len(values): + raise ValueError( + "CSR column/value count differs for R = {}: columns={} values={}".format( + rvec, len(columns), len(values) + ) + ) + + if len(row_ptr) != nbasis + 1: + raise ValueError( + "CSR row pointer count differs for R = {}: expected={} actual={}".format( + rvec, nbasis + 1, len(row_ptr) + ) + ) + + if row_ptr[0] != 0 or row_ptr[-1] != len(values): + raise ValueError( + "invalid CSR row pointers for R = {}".format(rvec) + ) + + if any(ptr < 0 or ptr > len(values) for ptr in row_ptr): + raise ValueError( + "CSR row pointer out of range for R = {}".format(rvec) + ) + + if any(lhs > rhs for lhs, rhs in zip(row_ptr, row_ptr[1:])): + raise ValueError( + "non-monotonic CSR row pointers for R = {}".format(rvec) + ) + + entries = {} + + for row in range(nbasis): + begin = row_ptr[row] + end = row_ptr[row + 1] + + for index in range(begin, end): + col = columns[index] + + if col < 0 or col >= nbasis: + raise ValueError( + "CSR column {} out of range for R = {}".format(col, rvec) + ) + + key = (row, col) + + if key in entries: + raise ValueError( + "duplicate CSR entry {} for R = {}".format(key, rvec) + ) + + entries[key] = values[index] + + return entries + + +def blocks_to_map(nbasis, blocks): + result = {} + + for rvec, values, columns, row_ptr in blocks: + if rvec in result: + raise ValueError("duplicate R block: {}".format(rvec)) + + result[rvec] = csr_to_entries( + nbasis, values, columns, row_ptr, rvec + ) + + return result + + +def read_pair(reference_filename, actual_filename): + errors = [] + + # Prefer real because the common nspin=1/2 H(R), S(R) outputs are real. + # Fall back to complex for spinor matrices. + for value_type in ("real", "complex"): + try: + reference_nbasis, reference_blocks = read_text_reference( + reference_filename, value_type + ) + actual_nbasis, actual_blocks = read_text_reference( + actual_filename, value_type + ) + + return ( + value_type, + reference_nbasis, + reference_blocks, + actual_nbasis, + actual_blocks, + ) + except (OSError, ValueError) as error: + errors.append("{}: {}".format(value_type, error)) + + raise ValueError( + "failed to parse files as real or complex CSR: {}".format( + "; ".join(errors) + ) + ) + + +def compare(reference_filename, actual_filename, tolerance): + ( + value_type, + reference_nbasis, + reference_blocks, + actual_nbasis, + actual_blocks, + ) = read_pair(reference_filename, actual_filename) + + if reference_nbasis != actual_nbasis: + raise ValueError( + "matrix dimension differs: reference={} actual={}".format( + reference_nbasis, actual_nbasis + ) + ) + + reference = blocks_to_map(reference_nbasis, reference_blocks) + actual = blocks_to_map(actual_nbasis, actual_blocks) + + max_difference = 0.0 + max_location = None + + # Missing R blocks and missing CSR entries are mathematically zero. + for rvec in sorted(set(reference) | set(actual)): + reference_entries = reference.get(rvec, {}) + actual_entries = actual.get(rvec, {}) + + for row_col in sorted(set(reference_entries) | set(actual_entries)): + expected = reference_entries.get(row_col, 0.0) + calculated = actual_entries.get(row_col, 0.0) + + if not math.isfinite(abs(expected)) or not math.isfinite(abs(calculated)): + row, col = row_col + raise ValueError( + "non-finite matrix value at R={}, row={}, col={}: " + "reference={} actual={}".format( + rvec, row, col, expected, calculated + ) + ) + + difference = abs(calculated - expected) + + if difference > max_difference: + max_difference = difference + max_location = (rvec, row_col, expected, calculated) + + if difference > tolerance: + row, col = row_col + raise ValueError( + "matrix value differs at R={}, row={}, col={}: " + "reference={} actual={} difference={} tolerance={}".format( + rvec, + row, + col, + expected, + calculated, + difference, + tolerance, + ) + ) + + if max_location is None: + print( + "H(R)/S(R) CSR comparison passed: " + "type={} max_abs_difference=0".format(value_type) + ) + else: + rvec, (row, col), expected, calculated = max_location + print( + "H(R)/S(R) CSR comparison passed: " + "type={} max_abs_difference={} " + "at R={}, row={}, col={}, reference={}, actual={}".format( + value_type, + max_difference, + rvec, + row, + col, + expected, + calculated, + ) + ) + + +def main(): + if len(sys.argv) != 4: + print( + "usage: compare_hsr_csr.py " + "REFERENCE_CSR ACTUAL_CSR ACCURACY" + ) + return 2 + + reference_filename = sys.argv[1] + actual_filename = sys.argv[2] + + try: + tolerance = 10.0 ** (-int(sys.argv[3])) + + compare( + reference_filename, + actual_filename, + tolerance, + ) + except (OSError, ValueError) as error: + print( + "failed to compare H(R)/S(R) CSR output: {}".format(error) + ) + return 1 + + return 0 + + +if __name__ == "__main__": + sys.exit(main())