#! /usr/bin/env python3
#-------------------------------*- python -*-----------------------------------
# ==  == ====== ====   ====    |
#                   \\     ||  | Multiphase Code Repository by HZDR
# ======   //   ||  || ===//   | Website: https://doi.org/10.14278/rodare.767
# ||  ||  //    ||  // || \\   | License: GPL-3.0-or-later
# ==  == ====== ====   ==  ==  |
#------------------------------------------------------------------------------
# License
#     This file is part of the Multiphase Code Repository by HZDR.
#
#     Copyright (C) 2025 by Helmholtz-Zentrum Dresden-Rossendorf e.V. (HZDR),
#     Website: https://hzdr.de
#
#     Multiphase Code Repository by HZDR is based on the free software for
#     computational fluid dynamics (CFD) from the OpenFOAM Foundation.
#     Copyright (C) 2025 by OpenFOAM Foundation, Website: https://openfoam.org
#
#     If you are interested in which files are original OpenFOAM Foundation
#     files, which OpenFOAM Foundation files were modified, and which files were
#     newly created, see FILES.md.
#
#     Multiphase Code Repository by HZDR is free software: you can redistribute
#     it and/or modify it under the terms of the GNU General Public License as
#     published by the Free Software Foundation, either version 3 of the
#     License, or (at your option) any later version.
#
#     Multiphase Code Repository by HZDR is distributed in the hope that it will
#     be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
#     Public License for more details.
#
#     You should have received a copy of the GNU General Public License along
#     with the Multiphase Code Repository by HZDR. If not, see
#     <http://www.gnu.org/licenses/>.
#
#------------------------------------------------------------------------------

"""Calculates the condition of a matrix. Checks if it is positive (negative)
definite and diagonal dominant or not."""

import argparse
import numpy as np

# Parse arguments
parser = argparse.ArgumentParser(description='Calculate matrix properties, \
e.g., definite, dominance and condition.')

# Add positional arguments
parser.add_argument('file', help='File with matrix coefficients')

args = parser.parse_args()

# Load matrix data
print("Loading matrix.")
a_mat = np.loadtxt(args.file)

# Eigenvalues
print("Check if matrix is definite.")
eigenvalues = np.linalg.eigvalsh(a_mat)
pos_def = np.all(eigenvalues > 0)
pos_semidef = np.all(eigenvalues >= 0)
neg_semidef = np.all(eigenvalues <= 0)
neg_def = np.all(eigenvalues < 0)

# Condition
print("Calculation of matrix condition.")
if pos_def or neg_def:
    print("  Using eigenvalues to calculate condition (matrix definite).")
    cond = np.amax(eigenvalues)/np.amin(eigenvalues)
else:
    print("  Using matrix norm to calculate condition.")
    cond = np.linalg.cond(a_mat)

# Diagonal dominance
print("Check diagonal dominance.")
diag_sum = np.diag(np.abs(a_mat))
row_sum = np.sum(np.abs(a_mat), axis=1)
weak_dominant = np.all(2.0*diag_sum >= row_sum)
strong_dominant = np.all(2.0*diag_sum > row_sum)
irreducible_dominant = np.any(2.0*diag_sum > row_sum)

print("Matrix properties are:")
print("  Condition                    : ", cond)
print("  Positive definite            : ", pos_def)
print("  Positive semidefinite        : ", pos_semidef)
print("  Negative definite            : ", neg_def)
print("  Negative semidefinite        : ", neg_semidef)
print("  Irreducible diagonal dominant: ", irreducible_dominant)
print("  Weak diagonal dominant       : ", weak_dominant)
print("  Strong diagonal dominant     : ", strong_dominant)
