DRAFT: Implementation of Galerkin Quadrature methods

.gitignore

@@ -76,4 +76,8 @@ doc/source/tutorials/**/*.pvtu
 doc/source/tutorials/**/*.vtu
 doc/source/tutorials/**/*.py
 
-doc/source/tutorials/**/.ipynb_checkpoints/
+doc/source/tutorials/**/.ipynb_checkpoints/
+
+pyOpenSn_env/
+pyopensn.egg-info/
+devFolder/

framework/math/math.cc

@@ -4,6 +4,10 @@
 #include "framework/math/math.h"
 #include "framework/runtime.h"
 #include <cassert>
+#include <petscmat.h>
+#include <petscksp.h>
+#include <iostream>
+#include <iomanip>
 
 namespace opensn
 {
@@ -49,6 +53,247 @@ OmegaToPhiThetaSafe(const Vector3& omega)
     return {2.0 * M_PI - acos(cos_phi), theta};
 }
 
+std::vector<std::vector<double>>
+Transpose(const std::vector<std::vector<double>>& matrix)
+{
+  // Check if matrix is empty
+  if (matrix.empty())
+  {
+    throw std::runtime_error("Cannot transpose empty matrix");
+  }
+
+  const size_t m = matrix.size();    // number of rows
+  const size_t n = matrix[0].size(); // number of columns
+
+  // Verify consistent row dimensions
+  for (size_t i = 1; i < m; ++i)
+  {
+    if (matrix[i].size() != n)
+    {
+      throw std::runtime_error("Matrix must have consistent row dimensions");
+    }
+  }
+
+  // Create transposed matrix with dimensions n x m
+  std::vector<std::vector<double>> transposed(n, std::vector<double>(m));
+
+  // Perform transpose: transposed[j][i] = matrix[i][j]
+  for (size_t i = 0; i < m; ++i)
+  {
+    for (size_t j = 0; j < n; ++j)
+    {
+      transposed[j][i] = matrix[i][j];
+    }
+  }
+
+  return transposed;
+}
+
+std::vector<std::vector<double>>
+InvertMatrix(const std::vector<std::vector<double>>& matrix)
+{
+  const PetscInt n = static_cast<PetscInt>(matrix.size());
+
+  // Check if matrix is square
+  if (n == 0)
+  {
+    throw std::runtime_error("Matrix must not be empty for inversion");
+  }
+
+  for (PetscInt i = 0; i < n; ++i)
+  {
+    if (static_cast<PetscInt>(matrix[i].size()) != n)
+    {
+      throw std::runtime_error("Matrix must be square for inversion");
+    }
+  }
+
+  Mat A, B, X;
+  MatFactorInfo info;
+
+  // Create dense matrix A directly
+  MatCreateSeqDense(PETSC_COMM_SELF, n, n, NULL, &A);
+
+  // Fill matrix A using array access (column-major storage)
+  PetscScalar* aArray;
+  MatDenseGetArrayWrite(A, &aArray);
+
+  for (PetscInt j = 0; j < n; ++j) // columns
+  {
+    for (PetscInt i = 0; i < n; ++i) // rows
+    {
+      aArray[j * n + i] = matrix[i][j]; // Column-major: col j, row i
+    }
+  }
+
+  MatDenseRestoreArrayWrite(A, &aArray);
+  MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
+  MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
+
+  // Create identity matrix B (right-hand side)
+  MatCreateSeqDense(PETSC_COMM_SELF, n, n, NULL, &B);
+
+  // Fill B with identity matrix using array access
+  PetscScalar* bArray;
+  MatDenseGetArrayWrite(B, &bArray);
+
+  // Initialize to zero first
+  for (PetscInt idx = 0; idx < n * n; ++idx)
+  {
+    bArray[idx] = 0.0;
+  }
+  // Set diagonal to 1
+  for (PetscInt i = 0; i < n; ++i)
+  {
+    bArray[i * n + i] = 1.0; // Column-major: diagonal elements
+  }
+
+  MatDenseRestoreArrayWrite(B, &bArray);
+  MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
+  MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
+
+  // Create solution matrix X
+  MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &X);
+
+  // Perform LU factorization directly
+  IS rowPerm, colPerm;
+  Mat F; // Factored matrix
+
+  MatGetOrdering(A, MATORDERINGND, &rowPerm, &colPerm);
+  MatFactorInfoInitialize(&info);
+  MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F);
+  MatLUFactorSymbolic(F, A, rowPerm, colPerm, &info);
+
+  // This is the critical operation that can fail for singular matrices
+  PetscErrorCode ierr = MatLUFactorNumeric(F, A, &info);
+  if (ierr != 0)
+  {
+    // Clean up before throwing
+    MatDestroy(&A);
+    MatDestroy(&B);
+    MatDestroy(&X);
+    MatDestroy(&F);
+    ISDestroy(&rowPerm);
+    ISDestroy(&colPerm);
+    throw std::runtime_error("LU factorization failed - matrix may be singular");
+  }
+
+  // Check if factorization succeeded (check for zero pivot)
+  PetscBool zeropivot;
+  MatFactorGetError(F, (MatFactorError*)&zeropivot);
+
+  if (zeropivot)
+  {
+    // Clean up before throwing
+    MatDestroy(&A);
+    MatDestroy(&B);
+    MatDestroy(&X);
+    MatDestroy(&F);
+    ISDestroy(&rowPerm);
+    ISDestroy(&colPerm);
+    throw std::runtime_error(
+      "Matrix is singular or nearly singular (zero pivot in LU factorization)");
+  }
+
+  // Solve AX = B (where B is identity, so X will be A^-1)
+  MatMatSolve(F, B, X);
+
+  // Extract solution back to std::vector format
+  std::vector<std::vector<double>> inverse(n, std::vector<double>(n));
+
+  // Get array from PETSc dense matrix (read-only)
+  const PetscScalar* xArray;
+  MatDenseGetArrayRead(X, &xArray);
+
+  // Copy data (PETSc uses column-major storage for dense matrices)
+  for (PetscInt i = 0; i < n; ++i)
+  {
+    for (PetscInt j = 0; j < n; ++j)
+    {
+      inverse[i][j] = PetscRealPart(xArray[j * n + i]); // Transpose due to column-major
+    }
+  }
+
+  MatDenseRestoreArrayRead(X, &xArray);
+
+  // Clean up PETSc objects
+  MatDestroy(&A);
+  MatDestroy(&B);
+  MatDestroy(&X);
+  MatDestroy(&F);
+  ISDestroy(&rowPerm);
+  ISDestroy(&colPerm);
+
+  return inverse;
+}
+
+std::vector<std::vector<double>>
+OrthogonalizeHouseholder(const std::vector<std::vector<double>>& matrix,
+                         const std::vector<double>& weights)
+{
+  // Check an empty matrix
+  if (matrix.empty())
+  {
+    throw std::runtime_error("Cannot orthogonalize empty matrix");
+  }
+
+  const size_t m = matrix.size();
+  const size_t n = matrix[0].size();
+
+  // Verify rectangular matrix
+  for (size_t i = 1; i < m; ++i)
+  {
+    if (matrix[i].size() != n)
+    {
+      throw std::runtime_error("Matrix must have consistent column dimensions");
+    }
+  }
+
+  // Orthogonalize columns using Modified Gram-Schmidt with weighted inner product
+  // This keeps each column as close to its original form as possible while ensuring orthogonality
+
+  // Create working copy of the matrix - we will orthogonalize its columns in place
+  std::vector<std::vector<double>> Q = matrix;
+
+  // Orthogonalize columns using Modified Gram-Schmidt
+  for (size_t j = 0; j < n; ++j)
+  {
+    // Orthogonalize column j against all previous columns
+    for (size_t i = 0; i < j; ++i)
+    {
+      // Compute weighted inner product <Q[:,i], Q[:,j]>
+      double dot_product = 0.0;
+      for (size_t row = 0; row < m; ++row)
+      {
+        dot_product += weights[row] * Q[row][i] * Q[row][j];
+      }
+
+      // Subtract projection: Q[:,j] -= dot_product * Q[:,i]
+      for (size_t row = 0; row < m; ++row)
+      {
+        Q[row][j] -= dot_product * Q[row][i];
+      }
+    }
+
+    // Normalize column j under weighted inner product
+    double norm_squared = 0.0;
+    for (size_t row = 0; row < m; ++row)
+    {
+      norm_squared += weights[row] * Q[row][j] * Q[row][j];
+    }
+    double norm = std::sqrt(norm_squared);
+
+    if (norm > 1e-14)
+    {
+      for (size_t row = 0; row < m; ++row)
+      {
+        Q[row][j] /= norm;
+      }
+    }
+  }
+  return Q;
+}
+
 void
 PrintVector(const std::vector<double>& x)
 {

framework/math/math.h

@@ -78,6 +78,45 @@ enum class NormType : int
  *  refines the view until the final linear search can be performed.*/
 int SampleCDF(double x, std::vector<double> cdf_bin);
 
+/**
+ * Transposes a matrix.
+ *
+ * Takes a matrix and returns its transpose. For a matrix A with dimensions m x n,
+ * the transpose A^T has dimensions n x m where A^T[j][i] = A[i][j].
+ *
+ * @param matrix Input matrix (m x n) represented as vector<vector<double>> where
+ *               matrix[i][j] represents row i, column j.
+ * @return Transposed matrix (n x m)
+ * @throws std::runtime_error if matrix is empty or has inconsistent row dimensions
+ */
+std::vector<std::vector<double>> Transpose(const std::vector<std::vector<double>>& matrix);
+
+/**
+ * Inverts a square matrix using PETSc LU decomposition.
+ *
+ * @param matrix The square matrix to invert
+ * @return The inverted matrix
+ * @throws std::runtime_error if the matrix is not square or inversion fails
+ */
+std::vector<std::vector<double>> InvertMatrix(const std::vector<std::vector<double>>& matrix);
+
+/**
+ * Orthogonalizes a matrix of column vectors using Householder reflections.
+ *
+ * Takes a matrix where each column represents a vector and returns a matrix
+ * with orthonormalized columns that span the same space. Uses Householder
+ * QR decomposition for numerical stability.
+ *
+ * @param matrix Input matrix (m x n) where each column is a vector to orthogonalize.
+ *               Matrix is represented as vector<vector<double>> where matrix[i][j]
+ *               represents row i, column j.
+ * @return Orthogonalized matrix (m x n) with orthonormal columns
+ * @throws std::runtime_error if matrix is empty or has inconsistent dimensions
+ */
+std::vector<std::vector<double>>
+OrthogonalizeHouseholder(const std::vector<std::vector<double>>& matrix,
+                         const std::vector<double>& weights);
+
 /// Computes the factorial of an integer.
 double Factorial(int x);