feat(debugUtils): added more sparse matrix debug utilities

utils/debugUtils/MFEMAnalysisUtils/MFEMAnalysis-cpp/src/include/mfem_smout.h

@@ -8,6 +8,11 @@
 #include "mfem.hpp"
 #include <iostream>
 #include <fstream>
+#include <vector>
+#include <array>
+#include <iomanip>
+#include <tuple>
+#include <ranges>
 
 /**
  * @brief Saves an mfem::SparseMatrix to a custom compact binary file (.csrbin).
@@ -29,6 +34,58 @@
  * - J array (int64_t * NNZ): CSR Column Indices
  * - Data array (double * NNZ): CSR Non-zero values
  */
+void write_sparse_matrix(const mfem::SparseMatrix &mat, std::ostream &outfile) {
+    // --- Get Data Pointers and Dimensions from MFEM Matrix ---
+    const int* mfem_I = mat.GetI();
+    const int* mfem_J = mat.GetJ();
+    const double* mfem_data = mat.GetData();
+
+    uint64_t height = static_cast<uint64_t>(mat.Height());
+    uint64_t width = static_cast<uint64_t>(mat.Width());
+    uint64_t nnz = static_cast<uint64_t>(mat.NumNonZeroElems());
+    uint64_t i_count = height + 1;
+    uint64_t j_count = nnz;
+    uint64_t data_count = nnz;
+
+
+    // --- Write Header ---
+    const char magic[4] = {'C', 'S', 'R', 'B'};
+    const uint8_t version = 1;
+    const uint8_t int_size = 8;
+    const uint8_t flt_size = 8;
+    const uint8_t reserved = 0;
+
+    outfile.write(magic, 4);
+    outfile.write(reinterpret_cast<const char*>(&version), 1);
+    outfile.write(reinterpret_cast<const char*>(&int_size), 1);
+    outfile.write(reinterpret_cast<const char*>(&flt_size), 1);
+    outfile.write(reinterpret_cast<const char*>(&reserved), 1);
+
+    outfile.write(reinterpret_cast<const char*>(&height), sizeof(height));
+    outfile.write(reinterpret_cast<const char*>(&width), sizeof(width));
+    outfile.write(reinterpret_cast<const char*>(&nnz), sizeof(nnz));
+
+    if (!outfile) throw std::runtime_error("Error writing header.");
+
+    // --- Write Arrays (Converting int to int64_t for I and J) ---
+    std::vector<int64_t> i_buffer(i_count);
+    for (uint64_t idx = 0; idx < i_count; ++idx) {
+        i_buffer[idx] = static_cast<int64_t>(mfem_I[idx]);
+    }
+    outfile.write(reinterpret_cast<const char*>(i_buffer.data()), i_count * sizeof(int64_t));
+    if (!outfile) throw std::runtime_error("Error writing I array.");
+
+    std::vector<int64_t> j_buffer(j_count);
+    for (uint64_t idx = 0; idx < j_count; ++idx) {
+        j_buffer[idx] = static_cast<int64_t>(mfem_J[idx]);
+    }
+    outfile.write(reinterpret_cast<const char*>(j_buffer.data()), j_count * sizeof(int64_t));
+    if (!outfile) throw std::runtime_error("Error writing J array.");
+
+    outfile.write(reinterpret_cast<const char*>(mfem_data), data_count * sizeof(double));
+    if (!outfile) throw std::runtime_error("Error writing Data array.");
+}
+
 bool saveSparseMatrixBinary(const mfem::SparseMatrix& mat, const std::string& filename) {
     std::ofstream outfile(filename, std::ios::binary | std::ios::trunc);
     if (!outfile) {
@@ -37,55 +94,7 @@ bool saveSparseMatrixBinary(const mfem::SparseMatrix& mat, const std::string& fi
     }
 
     try {
-        // --- Get Data Pointers and Dimensions from MFEM Matrix ---
+        write_sparse_matrix(mat, outfile);
-        const int* mfem_I = mat.GetI();
-        const int* mfem_J = mat.GetJ();
-        const double* mfem_data = mat.GetData();
-
-        uint64_t height = static_cast<uint64_t>(mat.Height());
-        uint64_t width = static_cast<uint64_t>(mat.Width());
-        uint64_t nnz = static_cast<uint64_t>(mat.NumNonZeroElems());
-        uint64_t i_count = height + 1;
-        uint64_t j_count = nnz;
-        uint64_t data_count = nnz;
-
-
-        // --- Write Header ---
-        const char magic[4] = {'C', 'S', 'R', 'B'};
-        const uint8_t version = 1;
-        const uint8_t int_size = 8;
-        const uint8_t flt_size = 8;
-        const uint8_t reserved = 0;
-
-        outfile.write(magic, 4);
-        outfile.write(reinterpret_cast<const char*>(&version), 1);
-        outfile.write(reinterpret_cast<const char*>(&int_size), 1);
-        outfile.write(reinterpret_cast<const char*>(&flt_size), 1);
-        outfile.write(reinterpret_cast<const char*>(&reserved), 1);
-
-        outfile.write(reinterpret_cast<const char*>(&height), sizeof(height));
-        outfile.write(reinterpret_cast<const char*>(&width), sizeof(width));
-        outfile.write(reinterpret_cast<const char*>(&nnz), sizeof(nnz));
-
-        if (!outfile) throw std::runtime_error("Error writing header.");
-
-        // --- Write Arrays (Converting int to int64_t for I and J) ---
-        std::vector<int64_t> i_buffer(i_count);
-        for (uint64_t idx = 0; idx < i_count; ++idx) {
-            i_buffer[idx] = static_cast<int64_t>(mfem_I[idx]);
-        }
-        outfile.write(reinterpret_cast<const char*>(i_buffer.data()), i_count * sizeof(int64_t));
-        if (!outfile) throw std::runtime_error("Error writing I array.");
-
-        std::vector<int64_t> j_buffer(j_count);
-        for (uint64_t idx = 0; idx < j_count; ++idx) {
-            j_buffer[idx] = static_cast<int64_t>(mfem_J[idx]);
-        }
-        outfile.write(reinterpret_cast<const char*>(j_buffer.data()), j_count * sizeof(int64_t));
-        if (!outfile) throw std::runtime_error("Error writing J array.");
-
-        outfile.write(reinterpret_cast<const char*>(mfem_data), data_count * sizeof(double));
-        if (!outfile) throw std::runtime_error("Error writing Data array.");
 
 
     } catch (const std::exception& e) {
@@ -163,4 +172,33 @@ void writeDenseMatrixToCSV(const std::string &filename, int precision, const mfe
     writeDenseMatrixToCSV(filename, precision, mat);
 }
 
+void saveBlockFormToBinary(std::vector<mfem::SparseMatrix *> &block_diags, std::vector<std::array<int, 2>> block, std::string filename) {
+    // First write a magic number and version
+
+    // --- Open the file ---
+    std::ofstream outfile(filename, std::ios::binary | std::ios::trunc);
+    if (!outfile) {
+        std::cerr << "Error: Cannot open file for writing: " << filename << std::endl;
+        return;
+    }
+
+    // --- Write Header ---
+    const char magic[4] = {'B', 'L', 'C', 'K'};
+    const char datastart[9] = {'D', 'A', 'T', 'A', 'S', 'T', 'A', 'R', 'T'};
+    const char dataend[7] = {'D', 'A', 'T', 'A', 'E', 'N', 'D'};
+    const uint8_t size = block_diags.size();
+
+    outfile.write(reinterpret_cast<const char*>(&magic), 4);
+    outfile.write(reinterpret_cast<const char*>(&size), sizeof(size));
+
+    for (const auto&& [block_diag, blockIDs] : std::views::zip(block_diags, block)) {
+        // Write the sparse matrix data
+        outfile.write(reinterpret_cast<const char*>(&datastart), 9);
+        outfile.write(reinterpret_cast<const char*>(&blockIDs), sizeof(blockIDs));
+        write_sparse_matrix(*block_diag, outfile);
+        outfile.write(reinterpret_cast<const char*>(&dataend), 7);
+    }
+
+}
+
 #endif //MFEM_SMOUT_H

utils/debugUtils/MFEMAnalysisUtils/SSEDebug/pyproject.toml

@@ -33,4 +33,7 @@ classifiers = [
 package-dir = {"" = "src"}
 
 [tool.setuptools.packages.find]
-where = ["src"]
+where = ["src"]
+
+[project.scripts]
+smanalyze = "SSEDebug.smRead.cli.interface:inspectSMMat"

utils/debugUtils/MFEMAnalysisUtils/SSEDebug/src/SSEDebug/smRead/cli/interface.py

@@ -0,0 +1,31 @@
+import argparse
+
+def inspectSMMat():
+    parser = argparse.ArgumentParser(description="Inspect SM matrix file")
+    parser.add_argument("filename", type=str, help="Path to the SM matrix file")
+    args = parser.parse_args()
+
+    try:
+        with open(args.filename, 'rb') as f:
+            magic = f.read(4)
+            if magic == b'BLCK':
+                print(f"{args.filename} is a valid block form SM matrix file.")
+                from SSEDebug.smRead.smread import loadBlockMatrix as matreader
+            if magic == b"CSRB":
+                print(f"{args.filename} is a valid CSR form SM matrix file.")
+                from SSEDebug.smRead.smread import loadSparseMatrix as matreader
+            else:
+                raise ValueError(f"Unknown file format: {magic}")
+
+        sm = matreader(args.filename)
+        from SSEDebug.smRead import analyze_sparse_matrix
+        analyze_sparse_matrix(sm)
+
+    except ValueError as e:
+        print(f"Invalid file format: {e}")
+    except FileNotFoundError:
+        print(f"File not found: {args.filename}")
+    except Exception as e:
+        print(f"An error occurred: {e}")
+    finally:
+        print("Finished inspecting the SM matrix file.")

utils/debugUtils/MFEMAnalysisUtils/SSEDebug/src/SSEDebug/smRead/smread.py

@@ -7,7 +7,7 @@ import scipy.sparse.linalg as spla # For matrix norm
 import time
 import os
 
-def loadSparseMatrixBinary(filename):
+def loadSparseMatrixBinary(f):
     """
     Loads a sparse matrix from the custom binary format (.csrbin).
 
@@ -27,74 +27,123 @@ def loadSparseMatrixBinary(filename):
     EXPECTED_VERSION = 1
 
     try:
-        with open(filename, 'rb') as f:
+        # --- Read Header ---
-            # --- Read Header ---
+        magic = f.read(4)
-            magic = f.read(4)
+        if magic != EXPECTED_MAGIC:
-            if magic != EXPECTED_MAGIC:
+            raise ValueError(f"Invalid magic number. Expected {EXPECTED_MAGIC}, got {magic}")
-                raise ValueError(f"Invalid magic number. Expected {EXPECTED_MAGIC}, got {magic}")
 
-            version, int_size_file, flt_size_file, reserved = struct.unpack('<BBBB', f.read(4))
+        version, int_size_file, flt_size_file, reserved = struct.unpack('<BBBB', f.read(4))
-            # '<' means little-endian, 'B' means unsigned char (1 byte)
+        # '<' means little-endian, 'B' means unsigned char (1 byte)
 
-            if version != EXPECTED_VERSION:
+        if version != EXPECTED_VERSION:
-                print(f"Warning: File version {version} differs from expected {EXPECTED_VERSION}.")
+            print(f"Warning: File version {version} differs from expected {EXPECTED_VERSION}.")
-            if int_size_file != INT_SIZE:
+        if int_size_file != INT_SIZE:
-                raise ValueError(f"Integer size mismatch. Expected {INT_SIZE}, file has {int_size_file}")
+            raise ValueError(f"Integer size mismatch. Expected {INT_SIZE}, file has {int_size_file}")
-            if flt_size_file != FLT_SIZE:
+        if flt_size_file != FLT_SIZE:
-                raise ValueError(f"Float size mismatch. Expected {FLT_SIZE}, file has {flt_size_file}")
+            raise ValueError(f"Float size mismatch. Expected {FLT_SIZE}, file has {flt_size_file}")
 
-            height, width, nnz = struct.unpack('<QQQ', f.read(24))
+        height, width, nnz = struct.unpack('<QQQ', f.read(24))
-            # '<' means little-endian, 'Q' means unsigned long long (8 bytes)
+        # '<' means little-endian, 'Q' means unsigned long long (8 bytes)
 
-            i_count = height + 1
+        i_count = height + 1
-            j_count = nnz
+        j_count = nnz
-            data_count = nnz
+        data_count = nnz
 
-            if nnz == 0: # Handle empty matrix case
+        if nnz == 0: # Handle empty matrix case
-                print("Warning: Matrix file contains zero non-zero elements.")
+            print("Warning: Matrix file contains zero non-zero elements.")
-                # Return an empty matrix with correct shape
+            # Return an empty matrix with correct shape
-                return sp.csr_matrix((height, width), dtype=np.float64)
+            return sp.csr_matrix((height, width), dtype=np.float64)
 
 
-            # --- Read Arrays ---
+        # --- Read Arrays ---
 
-            # Read I array (Row Pointers)
+        # Read I array (Row Pointers)
-            expected_i_bytes = i_count * INT_SIZE
+        expected_i_bytes = i_count * INT_SIZE
-            I_array = np.fromfile(f, dtype=np.int64, count=i_count) # Read as int64
+        I_array = np.fromfile(f, dtype=np.int64, count=i_count) # Read as int64
-            if I_array.size != i_count:
+        if I_array.size != i_count:
-                raise ValueError(f"Error reading I array. Expected {i_count} elements, read {I_array.size}. File truncated or corrupt?")
+            raise ValueError(f"Error reading I array. Expected {i_count} elements, read {I_array.size}. File truncated or corrupt?")
 
-            # Read J array (Column Indices)
+        # Read J array (Column Indices)
-            expected_j_bytes = j_count * INT_SIZE
+        expected_j_bytes = j_count * INT_SIZE
-            J_array = np.fromfile(f, dtype=np.int64, count=j_count) # Read as int64
+        J_array = np.fromfile(f, dtype=np.int64, count=j_count) # Read as int64
-            if J_array.size != j_count:
+        if J_array.size != j_count:
-                raise ValueError(f"Error reading J array. Expected {j_count} elements, read {J_array.size}. File truncated or corrupt?")
+            raise ValueError(f"Error reading J array. Expected {j_count} elements, read {J_array.size}. File truncated or corrupt?")
 
-            # Read Data array (Values)
+        # Read Data array (Values)
-            expected_data_bytes = data_count * FLT_SIZE
+        expected_data_bytes = data_count * FLT_SIZE
-            Data_array = np.fromfile(f, dtype=np.float64, count=data_count) # Read as float64
+        Data_array = np.fromfile(f, dtype=np.float64, count=data_count) # Read as float64
-            if Data_array.size != data_count:
+        if Data_array.size != data_count:
-                raise ValueError(f"Error reading Data array. Expected {data_count} elements, read {Data_array.size}. File truncated or corrupt?")
+            raise ValueError(f"Error reading Data array. Expected {data_count} elements, read {Data_array.size}. File truncated or corrupt?")
-
-            # --- Check for extra data ---
-            extra_data = f.read()
-            if extra_data:
-                print(f"Warning: {len(extra_data)} extra bytes found at the end of the file.")
 
 
-            # --- Construct SciPy CSR Matrix ---
+        # --- Construct SciPy CSR Matrix ---
-            sparse_matrix = sp.csr_matrix((Data_array, J_array, I_array), shape=(height, width))
+        sparse_matrix = sp.csr_matrix((Data_array, J_array, I_array), shape=(height, width))
 
-            if sparse_matrix.nnz != nnz:
+        if sparse_matrix.nnz != nnz:
-                print(f"Warning: NNZ mismatch after loading. Header NNZ: {nnz}, Scipy NNZ: {sparse_matrix.nnz}")
+            print(f"Warning: NNZ mismatch after loading. Header NNZ: {nnz}, Scipy NNZ: {sparse_matrix.nnz}")
 
 
-            return sparse_matrix
+        return sparse_matrix
 
-    except FileNotFoundError:
-        raise IOError(f"Error: File not found at {filename}")
     except Exception as e:
-        raise RuntimeError(f"An error occurred while reading {filename}: {e}")
+        raise RuntimeError(f"An error occurred while reading: {e}")
+
+def loadSparseMatrix(filename):
+    """
+    Loads a sparse matrix from the custom binary format (.csrbin).
+
+    Args:
+        filename (str): The path to the .csrbin file.
+
+    Returns:
+        scipy.sparse.csr_matrix: The loaded sparse matrix.
+
+    Raises:
+        ValueError: If the file format is incorrect or sizes don't match.
+        IOError: If the file cannot be read.
+    """
+    with open(filename, 'rb') as f:
+        # Check magic number
+        magic = f.read(4)
+        if magic != b'CSRB':
+            raise ValueError(f"Invalid magic number. Expected 'CSRB', got {magic}")
+
+        # Read the rest of the file
+        f.seek(0, 0)
+        sm = loadSparseMatrixBinary(f)
+
+    return sm
+def loadBlockMatrix(filename):
+    smList = list()
+    with open(filename, 'rb') as f:
+        f.seek(0, 2)
+        fileSize = f.tell()
+        f.seek(0, 0)
+        magic = f.read(4)
+        if magic != b'BLCK':
+            raise ValueError(f"Invalid magic number. Expected 'BLCK'. got {magic}")
+        size = struct.unpack('<B', f.read(1))[0]
+        print(f"Size: {size}")
+        while f.tell() < fileSize:
+            dataStartCard = f.read(9)
+            if dataStartCard != b'DATASTART':
+                raise ValueError(f"Invalid data start card. Expected 'DATASTART' Got {dataStartCard}.")
+            blockId = struct.unpack(f'<ii', f.read(8))
+            sm = loadSparseMatrixBinary(f)
+            smList.append((sm, blockId))
+            # unpack 2 ints as the block id
+            dataEndCard = f.read(7)
+            if dataEndCard != b'DATAEND':
+                raise ValueError(f"Invalid data end card. Expected 'DATAEND'. Got {dataEndCard}.")
+    outArray = np.empty(shape=(size, size), dtype=np.object_)
+
+    for sm, blockId in smList:
+        if blockId[0] >= size or blockId[1] >= size:
+            raise ValueError(f"Block ID {blockId} out of range. Size: {size}")
+        outArray[blockId[0], blockId[1]] = sm
+
+    # Check if all blocks are filled
+    return sp.bmat(outArray, format='csr')
 
 
 def analyze_sparse_matrix(sp_mat):
@@ -109,10 +158,6 @@ def analyze_sparse_matrix(sp_mat):
     print("Sparse Matrix Analysis Report")
     print("-" * 50)
 
-    if not isinstance(sp_mat, sp.spmatrix):
-        print("Error: Input is not a SciPy sparse matrix.")
-        return
-
     rows, cols = sp_mat.shape
     print(f"Size (Shape):           {rows} rows x {cols} columns")
 
@@ -129,8 +174,8 @@ def analyze_sparse_matrix(sp_mat):
     else:
         sparsity = 1.0
 
-    print(f"Non-zero elements (NNZ): {nnz}")
+    print(f"Non-zero elements (NNZ): {nnz} (~{nnz*8/(1024**2):.2f} MB)")
-    print(f"Total elements:         {total_elements}")
+    print(f"Total elements:         {total_elements} (~{total_elements*8/(1024**3):.2f} GB)")
     print(f"Sparsity:               {sparsity:.6%} (percentage of zeros)")
 
     if nnz == 0:
@@ -225,6 +270,14 @@ def load_and_analyze_sparse_matrix(filename: str):
     sm = loadSparseMatrixBinary(filename)
     analyze_sparse_matrix(sm)
 
+def compute_frobenius_distance(sparseMat):
+    identityMat = sp.eye(sparseMat.shape[0], sparseMat.shape[1], format='csr')
+    diffMat = sparseMat - identityMat
+    normDistance = np.sqrt(diffMat.data.dot(diffMat.data))
+    frobNormIdentity = np.sqrt(identityMat.shape[0])
+
+    return normDistance/frobNormIdentity
+
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(description="Simple tool to get some statistics about a sparse matrix from mfem")
     parser.add_argument("path", help="path to the output file", type=str)