Manually setting matrix entries

[mkofler96]

Hello,

I was wondering if it is possible to manually set entries of the matrix in the linear system.

If I look at the example given at the main page, I would like to manually change the sparse matrix AA.

# Form the linear system of equations (AX=B)
A = mfem.OperatorPtr()
B = mfem.Vector()
X = mfem.Vector()
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B)
print("Size of linear system: " + str(A.Height()))

# Solve the linear system using PCG and store the solution in x
AA = mfem.OperatorHandle2SparseMatrix(A)

# works
AA.Set(0, 0, 1)

# does not work
AA.Set(0, 10, 1)

While changinging the entry (0,0) works, because it is already present, adding a new entry at (0,10) gives the following error:

MFEM abort: Could not find entry for row = 0, col = 10
 ... in function: mfem::real_t& mfem::SparseMatrix::SearchRow(int, int)
 ... in file: /__w/PyMFEM/PyMFEM/PyMFEM/external/mfem/linalg/sparsemat.hpp:930

Traceback (most recent call last):
  File "/usr2/mkofler/PyMFEM/examplefrom_git.py", line 43, in <module>
    AA.Set(0, 10, 1)
RuntimeError: PyMFEM error (mfem::ErrorException): MFEM abort: Could not find entry for row = 0, col = 10
 ... in function: mfem::real_t& mfem::SparseMatrix::SearchRow(int, int)
 ... in file: /__w/PyMFEM/PyMFEM/PyMFEM/external/mfem/linalg/sparsemat.hpp:930

Is there another possibility of how I could do this? I was also looking for a way to convert the sparse matrix to a dense matrix and changing the entries there, but I have not found a way of doing this.

I would be happy if someone could help me out here and thanks for the great work.

[sshiraiwa]

Hi @mkofler96 I guess that the error indicates that you are trying to set the matrix element which is not occupied (non-zero). A sparsematrix in MFEM is a CSR matrix and it is expensive to change the matrix filling pattern.

My favorite recipe in this situation is to convert the MFEM matrix to scipy's sparse matrix, where you can manipulate more efficiently using other data format ( including the conversion to full dense matrix.) Then, I put it back to MFEM SparseMatrix.

Converting it to Scipy sparse matrix can be done like

###  this is found in mfem.common.sparse_utils
def sparsemat_to_scipycsr(mat, dtype):
     w, h = mat.Width(), mat.Height()
     I = mat.GetIArray()
     J = mat.GetJArray()
     data = mat.GetDataArray()
     mat.LoseData()
     m1 =csr_matrix((data, J, I), shape = (h, w),
                    dtype = dtype)
     return m1

An example of creating SparseMatrix from scratch is

def run_test():
    print("Test sparsemat module")
    indptr = np.array([0, 2, 3, 6], dtype=np.int32)
    indices = np.array([0, 2, 2, 0, 1, 2], dtype=np.int32)
    data = np.array([1, 2, 3, 4, 5, 6], dtype=float)
    
    mat = mfem.SparseMatrix([indptr, indices, data, 3,3])
    mat.Print()
    mat.PrintInfo()

@v-dobrev, @tzanio : if there is a different way, which can be done by just calling MFEM's API, please let us know.

[mkofler96]

Thank you @sshiraiwa vey much for the quick answer, this was exactly what I was looking for.

I also found that there is a method ToDenseMatrix on the mfem._ser.sparsemat.SparseMatrix class. But of course working directly on sparse matrices is way better.

[tzanio]

The reason you can't add entries in the SparseMatrix is that it is finalized in FormLinearSystem.

The code for FormLinearSystem in the serial conforming case (no AMR) and no Hybridization or static condensation is essentially

void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
                                    Vector &b, OperatorHandle &A, Vector &X,
                                    Vector &B, int copy_interior)
{
   EliminateVDofs(ess_tdof_list, diag_policy);
   Finalize(0);
   A.Reset(mat, false);
   EliminateVDofsInRHS(ess_tdof_list, x, b);
   X.MakeRef(x, 0, x.Size());
   B.MakeRef(b, 0, b.Size());
   if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
}

you can modify the entries in mat any time before the Finalize call