diff --git a/src/core/model/matrix-array.cc b/src/core/model/matrix-array.cc
index 72ec3f99a..854266dee 100644
--- a/src/core/model/matrix-array.cc
+++ b/src/core/model/matrix-array.cc
@@ -156,6 +156,70 @@ MatrixArray<T>::Transpose() const
     return res;
 }
 
+template <class T>
+MatrixArray<T>
+MatrixArray<T>::Determinant() const
+{
+    MatrixArray<T> res{1, 1, m_numPages};
+    NS_ASSERT_MSG(m_numRows == m_numCols, "Matrix is not square");
+    // In case of small matrices, we use a fast path
+    if (m_numRows == 1)
+    {
+        return *this;
+    }
+    // Calculate determinant for each matrix
+    for (size_t page = 0; page < m_numPages; ++page)
+    {
+        res(0, 0, page) = 0;
+        auto pageValues = GetPagePtr(page);
+
+        // Fast path for 2x2 matrices
+        if (m_numRows == 2)
+        {
+            res(0, 0, page) = pageValues[0] * pageValues[3] - pageValues[1] * pageValues[2];
+            continue;
+        }
+        for (size_t detN = 0; detN < m_numRows; detN++)
+        {
+            auto partDetP = T{0} + 1.0;
+            auto partDetN = T{0} + 1.0;
+            for (size_t row = 0; row < m_numRows; row++)
+            {
+                // Wraparound not to have to extend the matrix
+                // Positive determinant
+                size_t col = (row + detN) % m_numCols;
+                partDetP *= pageValues[row * m_numCols + col];
+
+                // Negative determinant
+                col = m_numCols - 1 - (row + detN) % m_numCols;
+                partDetN *= pageValues[row * m_numCols + col];
+            }
+            res(0, 0, page) += partDetP - partDetN;
+        }
+    }
+    return res;
+}
+
+template <class T>
+MatrixArray<T>
+MatrixArray<T>::FrobeniusNorm() const
+{
+    MatrixArray<T> res{1, 1, m_numPages};
+    for (size_t page = 0; page < m_numPages; ++page)
+    {
+        // Calculate the sum of squared absolute values of each matrix page
+        res[page] = 0;
+        auto pagePtr = this->GetPagePtr(page);
+        for (size_t i = 0; i < m_numRows * m_numCols; i++)
+        {
+            auto absVal = std::abs(pagePtr[i]);
+            res[page] += absVal * absVal;
+        }
+        res[page] = sqrt(res[page]);
+    }
+    return res;
+}
+
 template <class T>
 MatrixArray<T>
 MatrixArray<T>::MultiplyByLeftAndRightMatrix(const MatrixArray<T>& lMatrix,
@@ -229,6 +293,84 @@ MatrixArray<T>::HermitianTranspose() const
     return retMatrix;
 }
 
+template <class T>
+MatrixArray<T>
+MatrixArray<T>::MakeNCopies(size_t nCopies) const
+{
+    NS_ASSERT_MSG(m_numPages == 1, "The MatrixArray should have only one page to be copied.");
+    auto copiedMatrix = MatrixArray<T>{m_numRows, m_numCols, nCopies};
+    for (size_t copy = 0; copy < nCopies; copy++)
+    {
+        for (size_t i = 0; i < m_numRows * m_numCols; i++)
+        {
+            copiedMatrix.GetPagePtr(copy)[i] = m_values[i];
+        }
+    }
+    return copiedMatrix;
+}
+
+template <class T>
+MatrixArray<T>
+MatrixArray<T>::ExtractPage(size_t page) const
+{
+    NS_ASSERT_MSG(page < m_numPages, "The page to extract from the MatrixArray is out of bounds.");
+    auto extractedPage = MatrixArray<T>{m_numRows, m_numCols, 1};
+
+    for (size_t i = 0; i < m_numRows * m_numCols; ++i)
+    {
+        extractedPage.m_values[i] = GetPagePtr(page)[i];
+    }
+    return extractedPage;
+}
+
+template <class T>
+MatrixArray<T>
+MatrixArray<T>::JoinPages(const std::vector<MatrixArray<T>>& pages)
+{
+    auto jointMatrix =
+        MatrixArray<T>{pages.front().GetNumRows(), pages.front().GetNumCols(), pages.size()};
+    for (size_t page = 0; page < jointMatrix.GetNumPages(); page++)
+    {
+        NS_ASSERT_MSG(pages[page].GetNumRows() == jointMatrix.GetNumRows(),
+                      "All page matrices should have the same number of rows");
+        NS_ASSERT_MSG(pages[page].GetNumCols() == jointMatrix.GetNumCols(),
+                      "All page matrices should have the same number of columns");
+        NS_ASSERT_MSG(pages[page].GetNumPages() == 1,
+                      "All page matrices should have a single page");
+
+        size_t i = 0;
+        for (auto a : pages[page].GetValues())
+        {
+            jointMatrix.GetPagePtr(page)[i] = a;
+            i++;
+        }
+    }
+    return jointMatrix;
+}
+
+template <class T>
+MatrixArray<T>
+MatrixArray<T>::IdentityMatrix(const size_t size, const size_t pages)
+{
+    auto identityMatrix = MatrixArray<T>{size, size, pages};
+    for (std::size_t page = 0; page < pages; page++)
+    {
+        for (std::size_t i = 0; i < size; i++)
+        {
+            identityMatrix(i, i, page) = 1.0;
+        }
+    }
+    return identityMatrix;
+}
+
+template <class T>
+MatrixArray<T>
+MatrixArray<T>::IdentityMatrix(const MatrixArray& likeme)
+{
+    NS_ASSERT_MSG(likeme.GetNumRows() == likeme.GetNumCols(), "Template array is not square.");
+    return IdentityMatrix(likeme.GetNumRows(), likeme.GetNumPages());
+}
+
 template MatrixArray<std::complex<double>> MatrixArray<std::complex<double>>::HermitianTranspose()
     const;
 template class MatrixArray<std::complex<double>>;
diff --git a/src/core/model/matrix-array.h b/src/core/model/matrix-array.h
index 34305e214..3e3bf6ca0 100644
--- a/src/core/model/matrix-array.h
+++ b/src/core/model/matrix-array.h
@@ -189,6 +189,16 @@ class MatrixArray : public ValArray<T>
      * \return The resulting MatrixArray composed of the array of transposed matrices.
      */
     MatrixArray Transpose() const;
+    /**
+     * \brief This operator calculates a vector o determinants, one for each page
+     * \return The resulting MatrixArray with the determinants for each page.
+     */
+    MatrixArray Determinant() const;
+    /**
+     * \brief This operator calculates a vector of Frobenius norm, one for each page
+     * \return The resulting MatrixArray with the Frobenius norm for each page.
+     */
+    MatrixArray FrobeniusNorm() const;
     /**
      *\brief Multiply each matrix in the array by the left and the right matrix.
      * For each page of this MatrixArray the operation performed is
@@ -229,6 +239,40 @@ class MatrixArray : public ValArray<T>
               typename = std::enable_if_t<(std::is_same_v<T, std::complex<double>> && EnableBool)>>
     MatrixArray<T> HermitianTranspose() const;
 
+    /**
+     * \brief Function that copies the current 1-page matrix into a new matrix with n copies of the
+     * original matrix
+     * \param nCopies Number of copies to make of the input 1-page MatrixArray
+     * \return Returns a new matrix with n page copies
+     */
+    MatrixArray<T> MakeNCopies(size_t nCopies) const;
+    /**
+     * \brief Function extracts a page from a MatrixArray
+     * \param page Index of the page to be extracted from the n-page MatrixArray
+     * \return Returns an extracted page of the original MatrixArray
+     */
+    MatrixArray<T> ExtractPage(size_t page) const;
+    /**
+     * \brief Function joins multiple pages into a single MatrixArray
+     * \param pages Vector containing n 1-page MatrixArrays to be merged into a n-page MatrixArray
+     * \return Returns a joint MatrixArray
+     */
+    static MatrixArray<T> JoinPages(const std::vector<MatrixArray<T>>& pages);
+    /**
+     * \brief Function produces an identity MatrixArray with the specified size
+     * \param size Size of the output identity matrix
+     * \param pages Number of copies of the identity matrix
+     * \return Returns an identity MatrixArray
+     */
+    static MatrixArray<T> IdentityMatrix(const size_t size, const size_t pages = 1);
+    /**
+     * \brief Function produces an identity MatrixArray with the same size
+     * as the input MatrixArray
+     * \param likeme Input matrix used to determine the size of the identity matrix
+     * \return Returns an identity MatrixArray
+     */
+    static MatrixArray<T> IdentityMatrix(const MatrixArray& likeme);
+
   protected:
     // To simplify functions in MatrixArray that are using members from the template base class
     using ValArray<T>::m_numRows;
diff --git a/src/core/test/matrix-array-test-suite.cc b/src/core/test/matrix-array-test-suite.cc
index 044106675..3dc0e1ee0 100644
--- a/src/core/test/matrix-array-test-suite.cc
+++ b/src/core/test/matrix-array-test-suite.cc
@@ -40,6 +40,27 @@ namespace tests
 
 NS_LOG_COMPONENT_DEFINE("MatrixArrayTest");
 
+/**
+ * \brief Function casts an input valArray "in" (type IN) to an output valArray "out" (type T)
+ * \param in Input valarray to be casted
+ * \param out Output valarray to receive casted values
+ */
+template <typename IN, typename T>
+void
+CastStdValarray(const std::valarray<IN>& in, std::valarray<T>& out)
+{
+    // Ensure output valarray is the right size
+    if (out.size() != in.size())
+    {
+        out.resize(in.size());
+    }
+
+    // Perform the cast operation
+    std::transform(std::begin(in), std::end(in), std::begin(out), [](IN i) {
+        return static_cast<T>(i);
+    });
+}
+
 /**
  * \ingroup matrixArray-tests
  *  MatrixArray test case for testing constructors, operators and other functions
@@ -443,6 +464,154 @@ MatrixArrayTestCase<T>::DoRun()
     MatrixArray<T> m27 = MatrixArray<T>(std::move(jCasted));
     NS_LOG_INFO("m27.GetSize ()" << m27.GetSize());
     NS_TEST_ASSERT_MSG_EQ(jCastedSize, m27.GetSize(), "The number of elements are not equal.");
+
+    // test determinant
+    {
+        std::vector<std::pair<std::valarray<int>, T>> detTestCases{
+            // right-wraparound
+            {{1, 0, 7, 4, 2, 0, 6, 5, 3}, 62},
+            // left-wraparound
+            {{1, 4, 6, 0, 2, 5, 7, 0, 3}, 62},
+            // identity rank 3
+            {{1, 0, 0, 0, 1, 0, 0, 0, 1}, 1},
+            // identity rank 4
+            {{1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}, 1},
+            // permutation matrix rank 4
+            {{0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0}, -1},
+            // positive det matrix rank 2
+            {{36, -5, -5, 43}, 1523},
+            // single value matrix rank 1
+            {{1}, 1},
+        };
+        for (const auto& [detVal, detRef] : detTestCases)
+        {
+            std::valarray<T> detCast(detVal.size());
+            CastStdValarray(detVal, detCast);
+            auto side = sqrt(detVal.size());
+            MatrixArray<T> detMat = MatrixArray<T>(side, side, std::move(detCast));
+            NS_TEST_ASSERT_MSG_EQ(detMat.Determinant()(0),
+                                  static_cast<T>(detRef),
+                                  "The determinants are not equal.");
+        }
+    }
+    // clang-format off
+    std::valarray<int> multiPageMatrixValues{
+            // page 0: identity matrix
+            1, 0, 0,
+            0, 1, 0,
+            0, 0, 1,
+            // page 1: permutation matrix
+            0, 0, 1,
+            0, 1, 0,
+            1, 0, 0,
+            // page 2: random matrix
+            1, 4, 6,
+            0, 2, 5,
+            7, 0, 3,
+            // page 3: upper triangular
+            5, -9, 2,
+            0, -4, -10,
+            0, 0, -3,
+            // page 4: lower triangular
+            -7, 0, 0,
+            -1, -9, 0,
+            -5, 6, -5,
+            // page 5: zero diagonal
+            0, 1, 2,
+            1, 0, 1,
+            2, 1, 0,
+            // page 6: zero bidiagonal
+            0, 1, 0,
+            1, 0, 1,
+            0, 1, 0
+    };
+    // clang-format on
+    std::valarray<T> castMultiPageMatrixValues(multiPageMatrixValues.size());
+    CastStdValarray(multiPageMatrixValues, castMultiPageMatrixValues);
+    MatrixArray<T> multiPageMatrix = MatrixArray<T>(3,
+                                                    3,
+                                                    multiPageMatrixValues.size() / 9,
+                                                    std::move(castMultiPageMatrixValues));
+    // test the determinant in a multi-page MatrixArray
+    std::vector<int> determinants = {1, -1, 62, 60, -315, 4, 0};
+    auto det = multiPageMatrix.Determinant();
+    for (size_t page = 0; page < multiPageMatrix.GetNumPages(); page++)
+    {
+        NS_TEST_ASSERT_MSG_EQ(det(page),
+                              static_cast<T>(determinants[page]),
+                              "The determinants from the page " << std::to_string(page)
+                                                                << " are not equal.");
+    }
+
+    // test Frobenius norm in a multi-page MatrixArray
+    std::vector<double> fnorms = {sqrt(3), sqrt(3), 11.8322, 15.3297, 14.7309, 3.4641, 2};
+    auto frob = multiPageMatrix.FrobeniusNorm();
+    for (size_t page = 0; page < multiPageMatrix.GetNumPages(); page++)
+    {
+        NS_TEST_ASSERT_MSG_EQ_TOL(std::abs(frob(page)),
+                                  std::abs(static_cast<T>(fnorms[page])),
+                                  0.0001,
+                                  "The Frobenius norm from the page " << std::to_string(page)
+                                                                      << " are not equal.");
+    }
+
+    // test page copying
+    for (size_t noOfCopies = 1; noOfCopies < 4; noOfCopies++)
+    {
+        auto copies = m27.MakeNCopies(noOfCopies);
+        NS_TEST_ASSERT_MSG_EQ(copies.GetNumPages(),
+                              noOfCopies,
+                              "Creating " << std::to_string(noOfCopies) << " copies failed.");
+        NS_TEST_ASSERT_MSG_EQ(copies.GetNumRows(),
+                              m27.GetNumRows(),
+                              "The copy doesn't have the same number of rows as the original.");
+        NS_TEST_ASSERT_MSG_EQ(copies.GetNumCols(),
+                              m27.GetNumCols(),
+                              "The copy doesn't have the same number of columns as the original.");
+        for (size_t page = 0; page < copies.GetNumPages(); page++)
+        {
+            T diff{};
+            for (size_t row = 0; row < copies.GetNumRows(); row++)
+            {
+                for (size_t col = 0; col < copies.GetNumCols(); col++)
+                {
+                    diff += m27(row, col, 0) - copies(row, col, page);
+                }
+            }
+            NS_TEST_ASSERT_MSG_EQ(diff, T{}, "Mismatch in copied values.");
+        }
+    }
+
+    // test page 1 and 0 extraction
+    std::vector<MatrixArray<T>> pages{multiPageMatrix.ExtractPage(1),
+                                      multiPageMatrix.ExtractPage(0)};
+
+    // test page 1 and 0 joining
+    auto jointPagesMatrix = MatrixArray<T>::JoinPages(pages);
+    NS_TEST_ASSERT_MSG_EQ(jointPagesMatrix.GetNumPages(), 2, "Mismatch in number of join pages.");
+    for (size_t page = 0; page < jointPagesMatrix.GetNumPages(); page++)
+    {
+        T diff{};
+        for (size_t row = 0; row < jointPagesMatrix.GetNumRows(); row++)
+        {
+            for (size_t col = 0; col < jointPagesMatrix.GetNumCols(); col++)
+            {
+                diff += multiPageMatrix(row, col, 1 - page) - jointPagesMatrix(row, col, page);
+            }
+        }
+        NS_TEST_ASSERT_MSG_EQ(diff, T{}, "Mismatching pages.");
+    }
+
+    // test identity matrix
+    auto identityRank3Reference = multiPageMatrix.ExtractPage(0);
+    auto identityRank3 = MatrixArray<T>::IdentityMatrix(3);
+    NS_TEST_ASSERT_MSG_EQ(identityRank3, identityRank3Reference, "Mismatch in identity matrices.");
+
+    identityRank3 = MatrixArray<T>::IdentityMatrix(3, 10).ExtractPage(9);
+    NS_TEST_ASSERT_MSG_EQ(identityRank3, identityRank3Reference, "Mismatch in identity matrices.");
+
+    identityRank3 = MatrixArray<T>::IdentityMatrix(identityRank3Reference);
+    NS_TEST_ASSERT_MSG_EQ(identityRank3, identityRank3Reference, "Mismatch in identity matrices.");
 }
 
 /**