How to use the pymatsolver.PardisoSolver function in pymatsolver

srcs = []
    for ii in range(nSrc):
        rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
        srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))]

    survey = EM.TDEM.SurveyTDEM(srcs)

    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
    # prb.timeSteps = [1e-5]
    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
    # prb.timeSteps = [(1e-05, 100)]

    try:
        from pymatsolver import PardisoSolver
        prb.Solver = PardisoSolver
    except ImportError:
        prb.Solver  = SolverLU

    sigma = np.ones(mesh.nCz)*1e-8
    sigma[mesh.vectorCCz<0] = 1e-1
    sigma = np.log(sigma[active])

    prb.pair(survey)
    return prb, mesh, sigma

[np.r_[A0loc, 0.], np.r_[A1loc, 0.]],
                    rxloc, sighalf, earth_type="halfspace")

        rx = DC.Rx.Pole_ky(M)
        src0 = DC.Src.Dipole([rx], A0loc, A1loc)
        survey = DC.Survey_ky([src0])

        self.survey = survey
        self.mesh = mesh
        self.sigma = sigma
        self.data_ana = data_ana
        self.plotIt = False

        try:
            from pymatsolver import PardisoSolver
            self.Solver = PardisoSolver
        except ImportError:
            self.Solver = SolverLU

data_ana = EM.Analytics.DCAnalytic_Pole_Pole(
                    np.r_[A0loc, 0.],
                    rxloc, sighalf, earth_type="halfspace")

        rx = DC.Rx.Pole_ky(M)
        src0 = DC.Src.Pole([rx], A0loc)
        survey = DC.Survey_ky([src0])

        self.survey = survey
        self.mesh = mesh
        self.sigma = sigma
        self.data_ana = data_ana

        try:
            from pymatsolver import PardisoSolver
            self.Solver = PardisoSolver
        except ImportError:
            self.Solver = SolverLU

def setUp(self):

        mesh = Mesh.TensorMesh([20, 20, 20], "CCN")
        sigma = np.ones(mesh.nC)*1./100.
        actind = mesh.gridCC[:, 2] < -0.2
        # actMap = Maps.InjectActiveCells(mesh, actind, 0.)

        xyzM = Utils.ndgrid(np.ones_like(mesh.vectorCCx[:-1])*-0.4, np.ones_like(mesh.vectorCCy)*-0.4, np.r_[-0.3])
        xyzN = Utils.ndgrid(mesh.vectorCCx[1:], mesh.vectorCCy, np.r_[-0.3])

        problem = SP.Problem_CC(mesh, sigma=sigma, qMap=Maps.IdentityMap(mesh), Solver=PardisoSolver)
        rx = SP.Rx.Dipole(xyzN, xyzM)
        src = SP.Src.StreamingCurrents([rx], L=np.ones(mesh.nC), mesh=mesh,
                                       modelType="CurrentSource")
        survey = SP.Survey([src])
        survey.pair(problem)

        q = np.zeros(mesh.nC)
        inda = Utils.closestPoints(mesh, np.r_[-0.5, 0., -0.8])
        indb = Utils.closestPoints(mesh, np.r_[0.5, 0., -0.8])
        q[inda] = 1.
        q[indb] = -1.

        mSynth = q.copy()
        survey.makeSyntheticData(mSynth)

        # Now set up the problem to do some minimization

endl = np.array([[xmin, ymin, zmin], [xmax, ymax, zmax]])
survey = DCUtils.gen_DCIPsurvey(
    endl, survey_type="dipole-dipole", dim=mesh.dim,
    a=1, b=1, n=10, d2flag='2D'
)

# Setup Problem with exponential mapping and Active cells only in the core mesh
expmap = Maps.ExpMap(mesh)
mapactive = Maps.InjectActiveCells(
    mesh=mesh,  indActive=actind,
    valInactive=-5.
)
mapping = expmap * mapactive
problem = DC.Problem3D_CC(mesh,  sigmaMap=mapping, storeJ=True)
problem.pair(survey)
problem.Solver = PardisoSolver

survey.dpred(mtrue[actind])
survey.makeSyntheticData(mtrue[actind], std=0.01, force=True)

#####################
# Tikhonov Inversion#
#####################
m0 = np.median(ln_sigback) * np.ones(mapping.nP)
dmis = DataMisfit.l2_DataMisfit(survey)
regT = Regularization.Simple(mesh, indActive=actind)

# Personal preference for this solver with a Jacobi preconditioner
opt = Optimization.ProjectedGNCG(maxIter=20, lower=-10, upper=10,
                                 maxIterLS=20, maxIterCG=30, tolCG=1e-4)

opt.remember('xc')

ymin, ymax = 0., 0.
zmin, zmax = 0, 0
endl = np.array([[xmin, ymin, zmin], [xmax, ymax, zmax]])
survey = DCUtils.gen_DCIPsurvey(
    endl, survey_type="dipole-dipole", dim=mesh.dim,
    a=1, b=1, n=10, d2flag='2D'
)

# Setup Problem with exponential mapping and Active cells only in the core mesh
expmap = Maps.ExpMap(mesh)
mapactive = Maps.InjectActiveCells(mesh=mesh, indActive=actind,
                                   valInactive=-5.)
mapping = expmap * mapactive
problem = DC.Problem3D_CC(mesh, sigmaMap=mapping)
problem.pair(survey)
problem.Solver = PardisoSolver

survey.dpred(mtrue[actind])
survey.makeSyntheticData(mtrue[actind], std=0.05, force=True)


#####################
# Tikhonov Inversion#
#####################
m0 = np.median(ln_sigback) * np.ones(mapping.nP)
dmis = DataMisfit.l2_DataMisfit(survey)
regT = Regularization.Simple(
    mesh,
    alpha_s=1.,
    alpha_x=1.,
    alpha_y=1., indActive=actind
)

bzi = EM.FDEM.Rx.Point_b(
        np.array([[rxOffset, 0., src_height]]),
        orientation='z',
        component='imag'
    )

    # source location
    srcLoc = np.array([0., 0., src_height])
    srcList = [
        EM.FDEM.Src.MagDipole([bzr, bzi], freq, srcLoc, orientation='Z')
        for freq in freqs
    ]

    # construct a forward simulation
    survey = EM.FDEM.Survey(srcList)
    prb = EM.FDEM.Problem3D_b(mesh, sigmaMap=mapping, Solver=PardisoSolver)
    prb.pair(survey)

    # ------------------- Inversion ------------------- #
    # data misfit term
    survey.dobs = dobs
    dmisfit = DataMisfit.l2_DataMisfit(survey)
    uncert = abs(dobs) * std + floor
    dmisfit.W = 1./uncert

    # regularization
    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
    reg = Regularization.Simple(regMesh)
    reg.mref = mref

    # optimization
    opt = Optimization.InexactGaussNewton(maxIter=10)

sigma[active] = sig_half
sigma[layer] = sig_layer
#sigma[Rlayer] = sig_Rlayer
mtrue = np.log(sigma[active])

rxOffset = 1e-3
rx = EM.TDEM.Rx.Point_b(
np.array([[rxOffset, 0., 30]]),
np.logspace(-5, -3, 31),
'z'
)
src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 80]))
survey = EM.TDEM.Survey([src])
prb = EM.TDEM.Problem3D_b(mesh, sigmaMap=mapping)

prb.Solver = PardisoSolver
prb.timeSteps = [(1e-06, 20), (1e-05, 20), (0.0001, 20)]
prb.pair(survey)

# create observed data
std = 0.01

survey.dobs = survey.makeSyntheticData(mtrue, std)
survey.std = std
survey.eps = 1e-5*np.linalg.norm(survey.dobs)

dmisfit = DataMisfit.l2_DataMisfit(survey)
regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
reg = Regularization.Tikhonov(regMesh, alpha_s=1e-2, alpha_x=1.)
opt = Optimization.InexactGaussNewton(maxIter=10, LSshorten=0.5)
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)

endl = np.array([[xmin, ymin, zmin], [xmax, ymax, zmax]])
survey = DCUtils.gen_DCIPsurvey(
    endl, survey_type="dipole-dipole", dim=mesh.dim,
    a=1, b=1, n=10, d2flag='2D'
)

# Setup Problem with exponential mapping and Active cells only in the core mesh
expmap = Maps.ExpMap(mesh)
mapactive = Maps.InjectActiveCells(
    mesh=mesh,  indActive=actind,
    valInactive=-5.
)
mapping = expmap * mapactive
problem = DC.Problem3D_CC(mesh,  sigmaMap=mapping, storeJ=True)
problem.pair(survey)
problem.Solver = PardisoSolver

survey.dpred(mtrue[actind])
survey.makeSyntheticData(mtrue[actind], std=0.01, force=True)

#####################
# Tikhonov Inversion#
#####################
m0 = np.median(ln_sigback) * np.ones(mapping.nP)
dmis = DataMisfit.l2_DataMisfit(survey)
regT = Regularization.Simple(mesh, indActive=actind)

# Personal preference for this solver with a Jacobi preconditioner
opt = Optimization.ProjectedGNCG(maxIter=20, lower=-10, upper=10,
                                 maxIterLS=20, maxIterCG=30, tolCG=1e-4)

opt.remember('xc')

mesh = survey.setMesh(
    sigma=0.01, # approximate conductivity of the background
    max_depth_core=max_depth_core,  # extent of the core region of the mesh
    ncell_per_skind=10,  # number of cells per the smallest skin depth
    n_skind=2,  # number of skin depths that the mesh should extend to ensure the lowest-frequency fields have decayed
    core_meshType = "log",   # cell spacings in the core region of the mesh ("linear" or "log")
    max_hz_core=1000.  # If using a logarithmic core mesh, what is the maximum cell size?
 )

M = mesh

prob = MT1DProblem(
    mesh,  # The mesh contains the geometry, grids, etc necessary for constructing the discrete PDE system
    sigmaMap=Maps.ExpMap(mesh),  # in the inversion, we want to invert for log-conductivity (enforces postivity, electrical conductivity tends to vary logarithmically)
    verbose=False,  # print information as we are setting up and solving
    Solver=Solver  # solver to employ for solving Ax = b
)

# tell the problem and survey about each other so we can construct our matrix system
# and right hand-side
prob.pair(survey)

# start with nans so we can do a check to make sure all
# layer conductivities have been properly assigned
rho = np.ones(mesh.nC) * np.nan

# loop over each layer in the model and assign to mesh
for layer_top, rho_layer in zip(layer_tops, rho_layers):
    inds = mesh.vectorCCx < layer_top
    rho[inds] = rho_layer

sigma = 1./rho