Source code for Muscat.LinAlg.MumpsSolver

import numpy as np
from Muscat.Types import MuscatIndex, MuscatFloat




class MumpsSolver():
    """MumpsSolver class to keep the mumps functionnality within Muscat,
    However this relies on the python-mumps package
    """


    def __init__(self):
        from mumps import Context
        self.ctx = Context()



    def SetICNTL(self, index: int, var:int):
        """Set a ICNTL option """
        self.ctx.mumps_instance.icntl[index] = var




    def GetICNTL(self, index: int) -> int :
        """Get a ICNTL option """
        return self.ctx.mumps_instance.icntl[index]




    def SetCNTL(self, index: int , var):
        self.ctx.mumps_instance.cntl[index] = var




    def GetCNTL(self, index:int) -> float:
        return self.ctx.mumps_instance.cntl[index]




    def SetOp(self, matrix):
        self.ctx.set_matrix(matrix)
        self.DoAnalysisStep()
        self.DoFactorisationStep(reuse_analysis=True)





    def DoAnalysisStep(self):
        self.ctx.analyze()




    def DoFactorisationStep(self,reuse_analysis=False,ordering = "auto" ):
        self.ctx.factor(reuse_analysis=reuse_analysis,ordering = ordering)




    def Solve(self, rhs, sol ) :
        sol[:] = self.ctx.solve(rhs)
        return sol


    def __dealloc__(self):
        del self.ctx




def CheckIntegrity():
    # test data
    try:
        from mumps import Context
    except ImportError:
        return "skip"
    from scipy.sparse import coo_matrix

    # op
    data = np.array([1., 2.], dtype=MuscatFloat)
    i = np.array([0,1], dtype=MuscatIndex)
    j = np.array([0,1], dtype=MuscatIndex)
    m = coo_matrix((data,(i,j)), shape=(2,2) )

    # rsh
    rhs = np.zeros(2, dtype=MuscatFloat)
    rhs[0] = 1.0
    rhs[1] = 4.0

    # solution
    sol = np.zeros_like(rhs)

    solver = MumpsSolver()

    solver.SetOp(m)
    # test using int and
    solver.SetICNTL(1, 0)
    solver.SetICNTL(1,solver.GetICNTL(1))


    solver.SetCNTL(1 ,solver.GetCNTL(1) )



    solver.Solve(rhs, sol)

    error = sol - [1, 2]
    if np.max(np.abs(error)) > 1.e-16 :
        print(("sol",sol))
        print(("error",error))
        raise Exception("Error in the test case of Mumps")

    del solver

    return "ok"


if __name__ == '__main__':
    print(CheckIntegrity()) #pragma: no cover