How to use the pymatsolver.Pardiso function in pymatsolver

def get_prob(mesh, mapping, formulation):
    prb = getattr(EM.TDEM, 'Problem3D_{}'.format(formulation))(
        mesh, sigmaMap=mapping
    )
    prb.timeSteps = [(1e-3, 5), (1e-4, 5), (5e-5, 10), (5e-5, 10), (1e-4, 10)]
    prb.Solver = Solver
    return prb

# Setup the problem object
    sigma1d = M.r(sigBG, 'CC', 'CC', 'M')[0, 0, :]

    if expMap:
        problem = Problem3D_ePrimSec(M, sigmaPrimary=np.log(sigma1d))
        problem.sigmaMap = simpeg.Maps.ExpMap(problem.mesh)
        problem.model = np.log(sig)
    else:
        problem = Problem3D_ePrimSec(M, sigmaPrimary=sigma1d)
        problem.sigmaMap = simpeg.Maps.IdentityMap(problem.mesh)
        problem.model = sig
    problem.pair(survey)
    problem.verbose = False
    try:
        from pymatsolver import Pardiso
        problem.Solver = Pardiso
    except:
        pass

    return (survey, problem)

#     freq=freq
            # ), # less accurate
        ]

        survey = EM.FDEM.Survey(SrcList)

        sig = 1e-1
        sigma = np.ones(mesh.nC)*sig
        sigma[mesh.gridCC[:, 2] > 0] = 1e-8

        prb = EM.FDEM.Problem3D_b(mesh, sigma=sigma)
        prb.pair(survey)

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

        self.prb = prb
        self.mesh = mesh
        self.sig = sig

        print(' starting solve ...')
        u = self.prb.fields()
        print(' ... done')
        self.u = u

orientation=orientation))

        # primary
        self.primaryProblem = FDEM.Problem3D_j(meshp, sigmaMap=primaryMapping)
        self.primaryProblem.solver = Solver
        s_e = np.zeros(meshp.nF)
        inds = meshp.nFx + Utils.closestPoints(meshp, src_loc, gridLoc='Fz')
        s_e[inds] = 1./csz
        primarySrc = FDEM.Src.RawVec_e(
            self.rxlist, freq=freq, s_e=s_e/meshp.area
        )
        self.primarySurvey = FDEM.Survey([primarySrc])

        # Secondary Problem
        self.secondaryProblem = FDEM.Problem3D_e(meshs, sigmaMap=mapping)
        self.secondaryProblem.Solver = Solver
        self.secondarySrc = FDEM.Src.PrimSecMappedSigma(
                self.rxlist, freq, self.primaryProblem,
                self.primarySurvey, primaryMap2Meshs
        )
        self.secondarySurvey = FDEM.Survey([self.secondarySrc])
        self.secondaryProblem.pair(self.secondarySurvey)

        # Full 3D problem to compare with
        self.problem3D = FDEM.Problem3D_e(meshs, sigmaMap=mapping)
        self.problem3D.Solver = Solver
        s_e3D = np.zeros(meshs.nE)
        inds = (meshs.nEx + meshs.nEy +
                Utils.closestPoints(meshs, src_loc, gridLoc='Ez'))
        s_e3D[inds] = [1./(len(inds))] * len(inds)
        self.problem3D.model = model
        src3D = FDEM.Src.RawVec_e(self.rxlist, freq=freq, s_e=s_e3D)

)

    if srctype == "MagDipole":
        src = EM.TDEM.Src.MagDipole(
            [rx], waveform=EM.TDEM.Src.StepOffWaveform(),
            loc=np.array([0., 0., 0.])
        )
    elif srctype == "CircularLoop":
        src = EM.TDEM.Src.CircularLoop(
            [rx], waveform=EM.TDEM.Src.StepOffWaveform(),
            loc=np.array([0., 0., 0.]), radius=0.1
        )

    survey = EM.TDEM.Survey([src])
    prb = EM.TDEM.Problem3D_b(mesh, sigmaMap=mapping)
    prb.Solver = Solver

    prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40),
                     (0.0001, 40), (0.0005, 40)]

    sigma = np.ones(mesh.nCz)*1e-8
    sigma[active] = sig_half
    sigma = np.log(sigma[active])
    prb.pair(survey)
    if srctype == "MagDipole":
        bz_ana = mu_0*EM.Analytics.hzAnalyticDipoleT(rx.locs[0][0]+1e-3,
                                                     rx.times, sig_half)
    elif srctype == "CircularLoop":
        bz_ana = mu_0*EM.Analytics.hzAnalyticDipoleT(13, rx.times, sig_half)

    bz_calc = survey.dpred(sigma)
    ind = np.logical_and(rx.times > bounds[0], rx.times < bounds[1])

rxtimes = np.logspace(-4, -3, 20)
    rx = getattr(EM.TDEM.Rx, 'Point_{}'.format(rxcomp[:-1]))(
        locs=rxlocs, times=rxtimes, orientation=rxcomp[-1]
    )
    Aloc = np.r_[-10., 0., 0.]
    Bloc = np.r_[10., 0., 0.]
    srcloc = np.vstack((Aloc, Bloc))

    src = EM.TDEM.Src.LineCurrent([rx], loc=srcloc, waveform = EM.TDEM.Src.StepOffWaveform())
    survey = EM.TDEM.Survey([src])

    prb = getattr(EM.TDEM, 'Problem3D_{}'.format(prbtype))(mesh, sigmaMap=mapping)

    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]

    prb.Solver = Solver

    m = (
        np.log(1e-1)*np.ones(prb.sigmaMap.nP) +
        1e-3*np.random.randn(prb.sigmaMap.nP)
    )

    prb.pair(survey)
    mesh = mesh

    return prb, m, mesh

def setUpClass(self):

        print('\n------- Testing Primary Secondary Source HJ -> EB --------\n')
        # receivers
        self.rxlist = []
        for rxtype in ['b', 'e']:
            rx = getattr(FDEM.Rx, 'Point_{}'.format(rxtype))
            for orientation in ['x', 'y', 'z']:
                for comp in ['real', 'imag']:
                    self.rxlist.append(rx(rx_locs, component=comp,
                                       orientation=orientation))

        # primary
        self.primaryProblem = FDEM.Problem3D_j(meshp, sigmaMap=primaryMapping)
        self.primaryProblem.solver = Solver
        s_e = np.zeros(meshp.nF)
        inds = meshp.nFx + Utils.closestPoints(meshp, src_loc, gridLoc='Fz')
        s_e[inds] = 1./csz
        primarySrc = FDEM.Src.RawVec_e(
            self.rxlist, freq=freq, s_e=s_e/meshp.area
        )
        self.primarySurvey = FDEM.Survey([primarySrc])

        # Secondary Problem
        self.secondaryProblem = FDEM.Problem3D_e(meshs, sigmaMap=mapping)
        self.secondaryProblem.Solver = Solver
        self.secondarySrc = FDEM.Src.PrimSecMappedSigma(
                self.rxlist, freq, self.primaryProblem,
                self.primarySurvey, primaryMap2Meshs
        )
        self.secondarySurvey = FDEM.Survey([self.secondarySrc])

srcList = [
        EM.TDEM.Src.CircularLoop(
            rxList, loc=srcLoc, radius=radius,
            orientation='z',
            waveform=EM.TDEM.Src.VTEMWaveform(
                    offTime=offTime, peakTime=peakTime, a=3.
                )
        )
    ]
    # solve the problem at these times
    timeSteps = [
        (peakTime/5, 5), ((offTime-peakTime)/5, 5),
        (1e-5, 5), (5e-5, 5), (1e-4, 10), (5e-4, 15)
    ]
    prob = EM.TDEM.Problem3D_e(
        mesh, timeSteps=timeSteps, sigmaMap=mapping, Solver=Solver
    )
    survey = EM.TDEM.Survey(srcList)
    prob.pair(survey)

    src = srcList[0]
    rx = src.rxList[0]
    wave = []
    for time in prob.times:
        wave.append(src.waveform.eval(time))
    wave = np.hstack(wave)
    out = survey.dpred(m0)

    # plot the waveform
    fig = plt.figure(figsize=(5, 3))
    times_off = times-t0
    plt.plot(waveform_skytem[:, 0], waveform_skytem[:, 1], 'k.')

sigma[blk_ind] = sigma_block

    xmin, xmax = -200., 200.
    ymin, ymax = -200., 200.
    x = mesh.vectorCCx[np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCxymin, mesh.vectorCCy

np.array([[rxOffset, 0., src_height_resolve]]),
        orientation='z',
        component='imag'
    )

    # Set Source (In-phase and Quadrature)
    frequency_cp = resolve["frequency_cp"].value
    freqs = frequency_cp.copy()
    srcLoc = np.array([0., 0., src_height_resolve])
    srcList = [EM.FDEM.Src.MagDipole([bzr, bzi], freq, srcLoc, orientation='Z')
               for freq in freqs]

    # Set FDEM survey (In-phase and Quadrature)
    survey = EM.FDEM.Survey(srcList)
    prb = EM.FDEM.Problem3D_b(
        mesh, sigmaMap=mapping, Solver=Solver
    )
    prb.pair(survey)

    # ------------------ RESOLVE Inversion ------------------ #

    # Primary field
    bp = - mu_0/(4*np.pi*rxOffset**3)

    # Observed data
    cpi_inds = [0, 2, 6, 8, 10]
    cpq_inds = [1, 3, 7, 9, 11]
    dobs_re = np.c_[
        resolve["data"][rxind_resolve, :][cpi_inds],
        resolve["data"][rxind_resolve, :][cpq_inds]
    ].flatten() * bp * 1e-6