def init_solvers(self, nls_status=None, ls_conf=None, nls_conf=None, mtx=None, **kwargs):
        ls_conf = get_default(ls_conf, self.ls_conf, "you must set linear solver!")

        nls_conf = get_default(nls_conf, self.nls_conf, "you must set nonlinear solver!")

        ev = self.get_evaluator(**kwargs)
        ls = Solver.any_from_conf(ls_conf, mtx=mtx)
        nls = Solver.any_from_conf(nls_conf, evaluator=ev, lin_solver=ls, status=nls_status)
        self.solvers = Struct(name="solvers", ls=ls, nls=nls)

    def test_eigenvalue_solvers(self):
        eig_confs = self._list_eigenvalue_solvers(self.conf.solvers)

        n_eigs = 5

        ok = True
        tt = []
        for eig_conf in eig_confs:
            self.report(eig_conf.name)

            eig_solver = Solver.any_from_conf(eig_conf)

            t0 = time.clock()
            eigs, vecs = eig_solver(self.mtx, n_eigs=n_eigs, eigenvectors=True)
            tt.append((' '.join((eig_conf.name, eig_conf.kind)),
                       time.clock() - t0))

            self.report(eigs)

            _ok = nm.allclose(eigs.real, eigs_expected, rtol=0.0, atol=1e-8)
            ok = ok and (_ok or (eig_conf.kind in self.can_fail))

        tt.sort(key=lambda x: x[1])
        self.report('solution times:')
        for row in tt:
            self.report('%.2f [s]' % row[1], ':', row[0])

        return ok

    def presolve(self, mtx):
        """Prepare A^{-1} B^T for the Schur complement."""

        mtx_a = mtx['A']
        mtx_bt = mtx['BT']
        output('full A size: %.3f MB' % (8.0 * nm.prod(mtx_a.shape) / 1e6))
        output('full B size: %.3f MB' % (8.0 * nm.prod(mtx_bt.shape) / 1e6))

        ls = Solver.any_from_conf(self.problem.ls_conf,
                                   presolve=True, mtx=mtx_a)
        if self.mode == 'explicit':
            tt = time.clock()
            mtx_aibt = nm.zeros(mtx_bt.shape, dtype=mtx_bt.dtype)
            for ic in xrange(mtx_bt.shape[1]):
                mtx_aibt[:,ic] = ls(mtx_bt[:,ic].toarray().squeeze())
            output('mtx_aibt: %.2f s' % (time.clock() - tt))
            action_aibt = MatrixAction.from_array(mtx_aibt)
        else:
            ##
            # c: 30.08.2007, r: 13.02.2008
            def fun_aibt(vec):
                # Fix me for sparse mtx_bt...
                rhs = sc.dot(mtx_bt, vec)
                out = ls(rhs)
                return out
            action_aibt = MatrixAction.from_function(fun_aibt,
                                                    (mtx_a.shape[0],
                                                     mtx_bt.shape[1]),
                                                    nm.float64)
        mtx['action_aibt'] = action_aibt

def solve_optimize(conf, options):
    opts = conf.options
    trunk = io.get_trunk(conf.filename_mesh)
    data = {}

    dpb = ProblemDefinition.from_conf(conf, init_equations=False)
    equations = getattr(conf, "_".join(("equations_direct", opts.problem)))

    dpb.set_equations(equations)

    dpb.name = "direct"
    dpb.time_update(None)

    apb = dpb.copy("adjoint")
    equations = getattr(conf, "_".join(("equations_adjoint", opts.problem, opts.objective_function)))

    apb.set_equations(equations)
    apb.time_update(None)
    apb.ebcs.zero_dofs()
    apb.update_equations(None, ebcs=apb.ebcs)

    ls_conf = dpb.get_solver_conf(opts.ls)
    dnls_conf = dpb.get_solver_conf(opts.nls_direct)
    anls_conf = dpb.get_solver_conf(opts.nls_adjoint)
    opt_conf = dpb.get_solver_conf(opts.optimizer)

    dpb.init_solvers(ls_conf=ls_conf, nls_conf=dnls_conf)

    apb.init_solvers(ls_conf=ls_conf, nls_conf=anls_conf)

    shape_opt = so.ShapeOptimFlowCase.from_conf(conf, dpb, apb)
    design0 = shape_opt.dsg_vars.val
    shape_opt.cache = Struct(design=design0 + 100, state=None, i_mesh=-1)

    opt_status = IndexedStruct()
    optimizer = Solver.any_from_conf(
        opt_conf, obj_fun=so.obj_fun, obj_fun_grad=so.obj_fun_grad, status=opt_status, obj_args=(shape_opt, opts)
    )

    ##
    # State problem solution for the initial design.
    vec_dp0 = so.solve_problem_for_design(dpb, design0, shape_opt, opts)

    dpb.save_state(trunk + "_direct_initial.vtk", vec_dp0)

    ##
    # Optimize.
    des = optimizer(design0)
    print opt_status

    ##
    # Save final state (for "optimal" design).
    dpb.domain.mesh.write(trunk + "_opt.mesh", io="auto")
    dpb.save_state(trunk + "_direct_current.vtk", shape_opt.cache.state)

    print des

    def __init__(self, problem, options):
        self.mtx_mass = problem.evaluate(options.mass, mode='weak',
                                         auto_init=True, dw_mode='matrix')

        if options.lumped:
            raise NotImplementedError

        else:
            # Initialize solvers (and possibly presolve the matrix).
            self.ls = Solver.any_from_conf(problem.ls_conf, mtx=self.mtx_mass,
                                           presolve=True)

    def init_solvers(self, nls_status=None, ls_conf=None, nls_conf=None,
                     force=False):
        """
        Create and initialize solver instances.

        Parameters
        ----------
        nls_status : dict-like, IndexedStruct, optional
            The user-supplied object to hold nonlinear solver convergence
            statistics.
        ls_conf : Struct, optional
            The linear solver options.
        nls_conf : Struct, optional
            The nonlinear solver options.
        force : bool
            If True, re-create the solver instances even if they already exist
            in `self.nls` attribute.
        """
        if (self.nls is None) or force:
            ls_conf = get_default(ls_conf, self.ls_conf,
                                  'you must set linear solver!')

            nls_conf = get_default(nls_conf, self.nls_conf,
                                   'you must set nonlinear solver!')

            ls = Solver.any_from_conf(ls_conf, problem=self)

            ev = self.get_evaluator()

            if self.conf.options.get('ulf', False):
                self.nls_iter_hook = ev.new_ulf_iteration

            nls = Solver.any_from_conf(nls_conf, fun=ev.eval_residual,
                                       fun_grad=ev.eval_tangent_matrix,
                                       lin_solver=ls,
                                       iter_hook=self.nls_iter_hook,
                                       status=nls_status, problem=self)

            self.set_solver(nls)

def eig( mtx_a, mtx_b = None, num = None, eigenvectors = True,
         return_time = None, method = 'eig.scipy', **ckwargs ):

    kwargs = {'name' : 'aux', 'kind' : method}
    kwargs.update( ckwargs )
    conf = Struct( **kwargs )
    solver = Solver.any_from_conf( conf )

    status = {}
    out = solver( mtx_a, mtx_b, num, eigenvectors, status )
    if return_time is not None:
        return_time[0] = status['time']
        
    return out

    def init_solvers(self, nls_status=None, ls_conf=None, nls_conf=None,
                     mtx=None, presolve=False):
        """Create and initialize solvers."""
        ls_conf = get_default( ls_conf, self.ls_conf,
                               'you must set linear solver!' )

        nls_conf = get_default( nls_conf, self.nls_conf,
                              'you must set nonlinear solver!' )

        if presolve:
            tt = time.clock()
        if get_default_attr(ls_conf, 'needs_problem_instance', False):
            extra_args = {'problem' : self}
        else:
            extra_args = {}

        ls = Solver.any_from_conf(ls_conf, mtx=mtx, presolve=presolve,
                                  **extra_args)
        if presolve:
            tt = time.clock() - tt
            output('presolve: %.2f [s]' % tt)

        if get_default_attr(nls_conf, 'needs_problem_instance', False):
            extra_args = {'problem' : self}
        else:
            extra_args = {}
        ev = self.get_evaluator()

        if get_default_attr(self.conf.options, 'ulf', False):
            self.nls_iter_hook = ev.new_ulf_iteration

        nls = Solver.any_from_conf(nls_conf, fun=ev.eval_residual,
                                   fun_grad=ev.eval_tangent_matrix,
                                   lin_solver=ls, iter_hook=self.nls_iter_hook,
                                   status=nls_status, **extra_args)

        self.solvers = Struct( name = 'solvers', ls = ls, nls = nls )

    def prepare_matrices(self, problem):
        """
        A = K + B^T D^{-1} B
        """
        equations = problem.equations
        mtx = equations.eval_tangent_matrices(None, problem.mtx_a,
                                              by_blocks=True)

        ls = Solver.any_from_conf(problem.ls_conf,
                                  presolve=True, mtx=mtx['D'])

        mtx_b, mtx_m = mtx['B'], mtx['M']
        mtx_dib = nm.empty(mtx_b.shape, dtype=mtx_b.dtype)
        for ic in xrange(mtx_b.shape[1]):
            mtx_dib[:,ic] = ls(mtx_b[:,ic].toarray().squeeze())
        mtx_a = mtx['K'] + mtx_b.T * mtx_dib

        return mtx_a, mtx_m, mtx_dib

    def solve_eigen_problem(self):
        opts = self.app_options
        pb = self.problem

        pb.set_equations(pb.conf.equations)
        pb.time_update()

        output('assembling lhs...')
        tt = time.clock()
        mtx_a = pb.evaluate(pb.conf.equations['lhs'], mode='weak',
                            auto_init=True, dw_mode='matrix')
        output('...done in %.2f s' % (time.clock() - tt))

        if 'rhs' in pb.conf.equations:
            output('assembling rhs...')
            tt = time.clock()
            mtx_b = pb.evaluate(pb.conf.equations['rhs'], mode='weak',
                                dw_mode='matrix')
            output('...done in %.2f s' % (time.clock() - tt))

        else:
            mtx_b = None

        _n_eigs = get_default(opts.n_eigs, mtx_a.shape[0])

        output('solving eigenvalue problem for {} values...'.format(_n_eigs))
        eig = Solver.any_from_conf(pb.get_solver_conf(opts.evps))
        if opts.eigs_only:
            eigs = eig(mtx_a, mtx_b, opts.n_eigs, eigenvectors=False)
            svecs = None

        else:
            eigs, svecs = eig(mtx_a, mtx_b, opts.n_eigs, eigenvectors=True)

        output('...done')

        vecs = self.make_full(svecs)
        self.save_results(eigs, vecs)

        return Struct(pb=pb, eigs=eigs, vecs=vecs)

    def test_eigenvalue_solvers(self):
        eig_confs = self._list_eigenvalue_solvers(self.conf.solvers)

        n_eigs = 5

        ok = True
        tt = []
        for eig_conf in eig_confs:
            self.report(eig_conf.name)

            try:
                eig_solver = Solver.any_from_conf(eig_conf)

            except (ValueError, ImportError):
                if eig_conf.kind in self.can_fail:
                    continue

                else:
                    raise

            t0 = time.clock()
            eigs, vecs = eig_solver(self.mtx, n_eigs=n_eigs, eigenvectors=True)
            tt.append([' '.join((eig_conf.name, eig_conf.kind)),
                       time.clock() - t0])

            self.report(eigs)

            _ok = nm.allclose(eigs.real, eigs_expected, rtol=0.0, atol=1e-8)
            tt[-1].append(_ok)

            ok = ok and (_ok or (eig_conf.kind in self.can_fail)
                         or (eig_conf.name in self.can_miss))

        tt.sort(key=lambda x: x[1])
        self.report('solution times:')
        for row in tt:
            self.report('%.2f [s] : %s (ok: %s)' % (row[1], row[0], row[2]))

        return ok

    def solve_eigen_problem_n( self ):
        opts = self.app_options
        pb = self.problem

        dim = pb.domain.mesh.dim

        pb.set_equations( pb.conf.equations )
        pb.select_bcs( ebc_names = ['ZeroSurface'] )

        output( 'assembling rhs...' )
        tt = time.clock()
        mtx_b = pb.evaluate(pb.conf.equations['rhs'], mode='weak',
                            auto_init=True, dw_mode='matrix')
        output( '...done in %.2f s' % (time.clock() - tt) )
        assert_( nm.alltrue( nm.isfinite( mtx_b.data ) ) )

        ## mtx_b.save( 'b.txt', format='%d %d %.12f\n' )

        aux = pb.create_evaluable(pb.conf.equations['lhs'], mode='weak',
                                  dw_mode='matrix')
        mtx_a_equations, mtx_a_variables = aux

        if self.options.plot:
            log_conf = {
                'is_plot' : True,
                'aggregate' : 1,
                'yscales' : ['linear', 'log'],
            }
        else:
            log_conf = {
                'is_plot' : False,
            }
        log =  Log.from_conf( log_conf, ([r'$|F(x)|$'], [r'$|F(x)-x|$']) )

        file_output = Output('', opts.log_filename, combined = True)

        eig_conf = pb.get_solver_conf( opts.eigen_solver )
        eig_solver = Solver.any_from_conf( eig_conf )

        # Just to get the shape. Assumes one element group only!!!
        v_hxc_qp = pb.evaluate('dq_state_in_volume_qp.i1.Omega(Psi)')
        v_hxc_qp.fill(0.0)
        self.qp_shape = v_hxc_qp.shape
        vec_v_hxc = self._interp_to_nodes(v_hxc_qp)

        self.norm_v_hxc0 = nla.norm(vec_v_hxc)
        self.itercount = 0
        aux = wrap_function(self.iterate,
                            (eig_solver,
                             mtx_a_equations, mtx_a_variables,
                             mtx_b, log, file_output))
        ncalls, times, nonlin_v, results = aux

        # Create and call the DFT solver.
        dft_conf = pb.get_solver_conf(opts.dft_solver)
        dft_status = {}
        dft_solver = Solver.any_from_conf(dft_conf,
                                          fun = nonlin_v,
                                          status = dft_status)
        v_hxc_qp = dft_solver(v_hxc_qp.ravel())

        v_hxc_qp = nm.array(v_hxc_qp, dtype=nm.float64)
        v_hxc_qp.shape = self.qp_shape
        eigs, mtx_s_phi, vec_n, vec_v_h, v_ion_qp, v_xc_qp, v_hxc_qp = results
        output( 'DFT iteration time [s]:', dft_status['time_stats'] )

        fun = pb.materials['mat_v'].function
        variable = self.problem.create_variables(['scalar'])['scalar']
        vec_v_ion = fun(None, variable.field.get_coor(),
                        mode='qp')['V_ion'].squeeze()

        vec_v_xc = self._interp_to_nodes(v_xc_qp)
        vec_v_hxc = self._interp_to_nodes(v_hxc_qp)
        vec_v_sum = self._interp_to_nodes(v_hxc_qp + v_ion_qp)

        coor = pb.domain.get_mesh_coors()
        r2 = norm_l2_along_axis(coor, squared=True)
        vec_nr2 = vec_n * r2

        pb.select_bcs( ebc_names = ['ZeroSurface'] )
        mtx_phi = self.make_full( mtx_s_phi )

        out = {}
        update_state_to_output(out, pb, vec_n, 'n')
        update_state_to_output(out, pb, vec_nr2, 'nr2')
        update_state_to_output(out, pb, vec_v_h, 'V_h')
        update_state_to_output(out, pb, vec_v_xc, 'V_xc')
        update_state_to_output(out, pb, vec_v_ion, 'V_ion')
        update_state_to_output(out, pb, vec_v_hxc, 'V_hxc')
        update_state_to_output(out, pb, vec_v_sum, 'V_sum')
        self.save_results(eigs, mtx_phi, out=out)

        if self.options.plot:
            log( save_figure = opts.iter_fig_name )
            pause()
            log(finished=True)

        return Struct( pb = pb, eigs = eigs, mtx_phi = mtx_phi,
                       vec_n = vec_n, vec_nr2 = vec_nr2,
                       vec_v_h = vec_v_h, vec_v_xc = vec_v_xc )

    def test_semismooth_newton(self):
        import numpy as nm
        from sfepy.solvers import Solver

        ns = [0, 6, 2, 2]

        offsets = nm.cumsum(ns)
        nx = offsets[-1]

        iw = slice(offsets[0], offsets[1])
        ig = slice(offsets[1], offsets[2])
        il = slice(offsets[2], offsets[3])

        def fun_smooth(vec_x):
            xw = vec_x[iw]
            xg = vec_x[ig]
            xl = vec_x[il]

            m = self.ms
            rw = m['A'] * xw - m['Bb'].T * xg - self.m['fs'] 
            rg = m['Bb'] * xw + xl * xg

            rwg = nm.r_[rw, rg]

            return rwg

        def fun_smooth_grad(vec_x):
            xw = vec_x[iw]
            xl = vec_x[il]
            xg = vec_x[ig]

            m = self.m

            mzl = nm.zeros((6, 2), dtype=nm.float64)

            mw = nm.c_[m['A'], -m['Bb'].T, mzl]
            mg = nm.c_[m['Bb'], nm.diag(xl), nm.diag(xg)]

            mx = nm.r_[mw, mg]

            mx = sps.csr_matrix(mx)
            return mx

        def fun_a(vec_x):
            xw = vec_x[iw]
            xg = vec_x[ig]

            subsd = {}
            for ii, key in enumerate(self.w_names):
                subsd[key] = xw[ii]

            sn = eval_matrix(self.m['sn'], **subsd).squeeze()

            ra = nm.abs(xg) - fc * nm.abs(sn)

            return -ra

        def fun_a_grad(vec_x):
            xw = vec_x[iw]
            xg = vec_x[ig]
            xl = vec_x[il]

            subsd = {}
            for ii, key in enumerate(self.w_names):
                subsd[key] = xw[ii]

            md = eval_matrix(self.m['D'], **subsd)
            sn = eval_matrix(self.m['sn'], **subsd).squeeze()

            ma = nm.zeros((xl.shape[0], nx), dtype=nm.float64)
            ma[:,iw] = - fc * nm.sign(sn)[:,None] * md
            ma[:,ig] = nm.sign(xg)[:,None] * self.m['C']

            return -sps.csr_matrix(ma)

        def fun_b(vec_x):
            xl = vec_x[il]

            return xl

        def fun_b_grad(vec_x):
            xl = vec_x[il]

            mb = nm.zeros((xl.shape[0], nx), dtype=nm.float64)
            mb[:,il] = self.m['C']

            return sps.csr_matrix(mb)

        vec_x0 = 0.1 * nm.ones((nx,), dtype=nm.float64)

        lin_solver = Solver.any_from_conf(dict_to_struct(ls_conf))
        status = {}
        solver = Solver.any_from_conf(dict_to_struct(conf),
                                      fun_smooth=fun_smooth,
                                      fun_smooth_grad=fun_smooth_grad,
                                      fun_a=fun_a,
                                      fun_a_grad=fun_a_grad,
                                      fun_b=fun_b,
                                      fun_b_grad=fun_b_grad,
                                      lin_solver=lin_solver,
#.........这里部分代码省略.........

def setup_precond(mtx, problem):
    """
    Setup a preconditioner for `mtx`.
    """
    import scipy.sparse.linalg as spla
    from sfepy.solvers import Solver

    # Get active DOF indices for u, p.
    adi = problem.get_variables().adi
    iu = adi.indx['u']
    ip = adi.indx['p']

    # Get the diagonal blocks of the linear system matrix.
    K = mtx[iu, iu]
    M = mtx[ip, ip]

    # Create solvers for K, M blocks to be used in matvec_bj(). A different
    # solver for each block could be used.
    conf = problem.solver_confs['direct']
    # conf = problem.solver_confs['cg-s']
    # conf = problem.solver_confs['cg-p']
    # conf = problem.solver_confs['pyamg']
    ls1 = Solver.any_from_conf(conf, mtx=K, context=problem)
    ls2 = Solver.any_from_conf(conf, mtx=M, context=problem)

    def matvec_bj(vec):
        """
        The application of the Block Jacobi preconditioner.

        The exact version (as with the `'direct'` solver) can be obtained also
        by using the following PETSs command-line options, together with the
        `'iterative-p'` solver::

          -ksp_monitor -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_u_ksp_type preonly -fieldsplit_u_pc_type lu -fieldsplit_p_ksp_type preonly -fieldsplit_p_pc_type lu

        The inexact version (20 iterations of a CG solver for each block, as
        with the `'cg-s'` or `'cg-p'` solvers) can be obtained also by using
        the following PETSs command-line options, together with the
        `'iterative-p'` solver::

          -ksp_monitor -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_u_ksp_type cg -fieldsplit_u_pc_type none -fieldsplit_p_ksp_type cg -fieldsplit_p_pc_type none -fieldsplit_u_ksp_max_it 20 -fieldsplit_p_ksp_max_it 20
        """
        vu = ls1(vec[iu])
        vp = ls2(vec[ip])

        return nm.r_[vu, vp]

    def matvec_j(vec):
        """
        The application of the Jacobi (diagonal) preconditioner.

        The same effect can be obtained also by using the following PETSs
        command-line options, together with the `'iterative-p'` solver::

          -ksp_monitor -pc_type jacobi
        """
        D = mtx.diagonal()

        return vec / D

    # Create the preconditioner, using one of matvec_bj() or matvec_j().
    precond = Struct(name='precond', shape=mtx.shape, matvec=matvec_bj)
    precond = spla.aslinearoperator(precond)

    return precond

def solve_navier_stokes(conf, options):
    opts = conf.options

    dpb = ProblemDefinition.from_conf(conf, init_equations=False)
    equations = getattr(conf, "_".join(("equations_direct", opts.problem)))
    dpb.set_equations(equations)

    ls_conf = dpb.get_solver_conf(opts.ls)
    nls_conf = dpb.get_solver_conf(opts.nls_direct)

    method = opts.direct_method
    if method == "stationary":
        data = {}
        dpb.time_update(None)
        state_dp = dpb.solve(nls_conf=nls_conf)

    elif method == "transient":
        ls = Solver.any_from_conf(ls_conf)
        ts_conf = dpb.get_solver_conf(opts.ts_direct)

        data = {"ts": Struct(dt=ts_conf.dt)}

        # Plug in mass term.
        mequations = {}
        for key, eq in equations.iteritems():
            if "dw_div_grad" in eq:
                eq = "+".join((ts_conf.mass_term, eq)).replace("++", "+")
            mequations[key] = eq

        if ts_conf.stokes_init:
            state_dp0 = solve_stokes(dpb, conf.equations_direct_stokes, nls_conf)
            dpb.set_equations(mequations)
        else:
            dpb.set_equations(mequations)
            state_dp0 = dpb.create_state()
            dpb.time_update(None)
            state_dp0.apply_ebc()

        from sfepy.base.log import Log

        log = Log.from_conf(Struct(is_plot=True), ([r"$||u||$"], [r"$||p||$"]))

        output("Navier-Stokes...")
        ev = BasicEvaluator(dpb, ts=Struct(dt=ts_conf.dt))
        nls = Solver.any_from_conf(nls_conf, evaluator=ev, lin_solver=ls)

        n_step = ts_conf.n_step
        step = 0
        while 1:
            for ii in xrange(n_step):
                output(step)

                vec_u = state_dp0("w")
                vec_p = state_dp0("r")
                log(nm.linalg.norm(vec_u), nm.linalg.norm(vec_p))

                dpb.variables.set_data_from_state("w_0", state_dp0(), "w")
                vec_dp = nls(state_dp0())

                step += 1
                state_dp = state_dp0.copy()
                state_dp.set_reduced(vec_dp)

                state_dp0 = state_dp

            if ts_conf.interactive:
                try:
                    n_step = int(raw_input("continue: "))
                    if n_step <= 0:
                        break
                except:
                    break

        vec_u = state_dp("w")
        vec_p = state_dp("r")
        log(nm.linalg.norm(vec_u), nm.linalg.norm(vec_p), finished=True)

    else:
        raise "unknown Navier-Stokes solution method (%s)!" % method

    return dpb, state_dp, data

    def get_time_solver(self, ts_conf=None, **kwargs):
        ts_conf = get_default(ts_conf, self.ts_conf, "you must set time-stepping solver!")

        return Solver.any_from_conf(ts_conf, **kwargs)

    def test_ls_reuse(self):
        import numpy as nm
        from sfepy.solvers import Solver
        from sfepy.discrete.state import State

        self.problem.init_solvers(ls_conf=self.problem.solver_confs['d00'])
        nls = self.problem.get_nls()

        state0 = State(self.problem.equations.variables)
        state0.apply_ebc()
        vec0 = state0.get_reduced()

        self.problem.update_materials()

        rhs = nls.fun(vec0)
        mtx = nls.fun_grad(vec0)

        ok = True
        for name in ['i12', 'i01']:
            solver_conf = self.problem.solver_confs[name]
            method = solver_conf.get('method', '')
            precond = solver_conf.get('precond', '')
            name = ' '.join((solver_conf.name, solver_conf.kind,
                             method, precond)).rstrip()
            self.report(name)
            try:
                ls = Solver.any_from_conf(solver_conf)

            except:
                self.report('skipped!')
                continue

            conf = ls.conf.copy()
            conf.force_reuse = True

            sol00 = ls(rhs, mtx=mtx, conf=conf)
            digest00 = ls.mtx_digest

            sol0 = ls(rhs, mtx=mtx)
            digest0 = ls.mtx_digest

            sol1 = ls(rhs, mtx=2*mtx, conf=conf)
            digest1 = ls.mtx_digest

            sol2 = ls(rhs, mtx=2*mtx)
            digest2 = ls.mtx_digest
            ls(rhs, mtx=2*mtx)
            digest3 = ls.mtx_digest

            _ok = digest00 != digest0
            self.report(digest00, '!=', digest0, ':', _ok); ok = ok and _ok
            _ok = digest0 == digest1
            self.report(digest0, '==', digest1, ':', _ok); ok = ok and _ok
            _ok = digest1 != digest2
            self.report(digest1, '!=', digest2, ':', _ok); ok = ok and _ok
            _ok = digest2[1] == digest3[1]
            self.report(digest2[1], '==', digest3[1], ':', _ok); ok = ok and _ok
            _ok = nm.allclose(sol00, sol0, atol=1e-12, rtol=0.0)
            self.report('sol00 == sol0:', _ok); ok = ok and _ok
            _ok = nm.allclose(sol0, sol1, atol=1e-12, rtol=0.0)
            self.report('sol0 == sol1:', _ok); ok = ok and _ok
            _ok = nm.allclose(sol0, 2 * sol2, atol=1e-12, rtol=0.0)
            self.report('sol0 == 2 * sol2:', _ok); ok = ok and _ok

        return ok

    def solve_eigen_problem( self, ofn_trunk = None, post_process_hook = None ):

        if self.cached_evp is not None:
            return self.cached_evp

        problem = self.problem
        ofn_trunk = get_default( ofn_trunk, problem.ofn_trunk,
                                 'output file name trunk missing!' )
        post_process_hook = get_default( post_process_hook,
                                         self.post_process_hook )
 
        conf = self.conf
        
        eig_problem = self.app_options.eig_problem
        if eig_problem in ['simple', 'simple_liquid']:
            problem.set_equations( conf.equations )
            problem.time_update()

            mtx_a = problem.evaluate(conf.equations['lhs'], mode='weak',
                                     auto_init=True, dw_mode='matrix')

            mtx_m = problem.evaluate(conf.equations['rhs'], mode='weak',
                                     dw_mode='matrix')

        elif eig_problem == 'schur':
            # A = K + B^T D^{-1} B.
            mtx = assemble_by_blocks( conf.equations, self.problem,
                                      ebcs = conf.ebcs,
                                      epbcs = conf.epbcs )
            problem.set_equations( conf.equations )
            problem.time_update()

            ls = Solver.any_from_conf( problem.ls_conf,
                                       presolve = True, mtx = mtx['D'] )

            mtx_b, mtx_m = mtx['B'], mtx['M']
            mtx_dib = nm.empty( mtx_b.shape, dtype = mtx_b.dtype )
            for ic in xrange( mtx_b.shape[1] ):
                mtx_dib[:,ic] = ls( mtx_b[:,ic].toarray().squeeze() )
            mtx_a = mtx['K'] + mtx_b.T * mtx_dib

        else:
            raise NotImplementedError

##     from sfepy.base.plotutils import spy, plt
##     spy( mtx_b, eps = 1e-12 )
##     plt.show()
##     mtx_a.save( 'a.txt', format='%d %d %.12f\n' )
##     mtx_b.save( 'b.txt', format='%d %d %.12f\n' )
##     pause()

        output( 'computing resonance frequencies...' )
        tt = [0]

        if isinstance( mtx_a, sc.sparse.spmatrix ):
            mtx_a = mtx_a.toarray()
        if isinstance( mtx_m, sc.sparse.spmatrix ):
            mtx_m = mtx_m.toarray()

        eigs, mtx_s_phi = eig(mtx_a, mtx_m, return_time=tt,
                              method=self.app_options.eigensolver)
        eigs[eigs<0.0] = 0.0
        output( '...done in %.2f s' % tt[0] )
        output( 'original eigenfrequencies:' )        
        output( eigs )
        opts = self.app_options
        epsilon2 = opts.scale_epsilon * opts.scale_epsilon
        eigs_rescaled = (opts.elasticity_contrast / epsilon2)  * eigs
        output( 'rescaled eigenfrequencies:' )        
        output( eigs_rescaled )
        output( 'number of eigenfrequencies: %d' % eigs.shape[0] )

        try:
            assert_( nm.isfinite( eigs ).all() )
        except ValueError:
            debug()

        # B-orthogonality check.
##         print nm.dot( mtx_s_phi[:,5], nm.dot( mtx_m, mtx_s_phi[:,5] ) )
##         print nm.dot( mtx_s_phi[:,5], nm.dot( mtx_m, mtx_s_phi[:,0] ) )
##         debug()

        n_eigs = eigs.shape[0]

        variables = problem.get_variables()

        mtx_phi = nm.empty( (variables.di.ptr[-1], mtx_s_phi.shape[1]),
                           dtype = nm.float64 )

        make_full = variables.make_full_vec
        if eig_problem in ['simple', 'simple_liquid']:
            for ii in xrange( n_eigs ):
                mtx_phi[:,ii] = make_full( mtx_s_phi[:,ii] )
            eig_vectors = mtx_phi
            
        elif eig_problem == 'schur':
            # Update also eliminated variables.
            schur = self.app_options.schur
            primary_var = schur['primary_var']
            eliminated_var = schur['eliminated_var']
#.........这里部分代码省略.........

def solve_navier_stokes(conf, options):
    opts = conf.options

    dpb = ProblemDefinition.from_conf(conf, init_equations=False)
    equations = getattr(conf, '_'.join(('equations_direct', opts.problem)))
    dpb.set_equations(equations)

    ls_conf = dpb.get_solver_conf( opts.ls )
    nls_conf = dpb.get_solver_conf(opts.nls_direct)

    method = opts.direct_method
    if method == 'stationary':
        data = {}
        dpb.time_update(None)
        vec_dp = dpb.solve(nls_conf=nls_conf)

    elif method == 'transient':
        ls = Solver.any_from_conf( ls_conf )
        ts_conf = dpb.get_solver_conf( opts.ts_direct )

        data = {'ts' : Struct( dt = ts_conf.dt )}

        # Plug in mass term.
        mequations = {}
        for key, eq in equations.iteritems():
            if 'dw_div_grad' in eq:
                eq = '+'.join( (ts_conf.mass_term, eq) ).replace( '++', '+')
            mequations[key] = eq

        if ts_conf.stokes_init:
            vec_dp0 = solve_stokes( dpb, conf.equations_direct_stokes, nls_conf )
            dpb.set_equations( mequations )
        else:
            dpb.set_equations( mequations )
            vec_dp0 = dpb.create_state_vector()
            dpb.time_update( None )
            dpb.apply_ebc( vec_dp0 )

        from sfepy.base.log import Log

        log = Log.from_conf( Struct( is_plot = True ),
                            ([r'$||u||$'], [r'$||p||$']) )

        output( 'Navier-Stokes...' )
        ev = BasicEvaluator( dpb, ts = Struct( dt = ts_conf.dt ) )
        nls = Solver.any_from_conf( nls_conf, evaluator = ev, lin_solver = ls )

        n_step = ts_conf.n_step
        step = 0
        while 1:
            for ii in xrange( n_step ):
                output( step )

                vec_u = dpb.variables.get_state_part_view( vec_dp0, 'w' )
                vec_p = dpb.variables.get_state_part_view( vec_dp0, 'r' )
                log( nm.linalg.norm( vec_u ), nm.linalg.norm( vec_p ) )

                dpb.variables.non_state_data_from_state( 'w_0', vec_dp0, 'w' )
                vec_dp = nls( vec_dp0 )

                step += 1
                vec_dp0 = vec_dp.copy()

            if ts_conf.interactive:
                try:
                    n_step = int( raw_input( 'continue: ' ) )
                    if n_step <= 0: break
                except:
                    break

        vec_u = dpb.variables.get_state_part_view( vec_dp, 'w' )
        vec_p = dpb.variables.get_state_part_view( vec_dp, 'r' )
        log( nm.linalg.norm( vec_u ), nm.linalg.norm( vec_p ), finished = True )

    else:
        raise 'unknown Navier-Stokes solution method (%s)!'  % method
    
    return dpb, vec_dp, data

#.........这里部分代码省略.........

    pars = apply_units(options.pars, options.unit_multipliers)
    output('material parameters with applied unit multipliers:')
    output(pars)

    if options.mode == 'omega':
        rng = copy(options.range)
        rng[:2] = apply_unit_multipliers(options.range[:2],
                                         ['wave_number', 'wave_number'],
                                         options.unit_multipliers)
        output('wave number range with applied unit multipliers:', rng)

    else:
        rng = copy(options.range)
        rng[:2] = apply_unit_multipliers(options.range[:2],
                                         ['frequency', 'frequency'],
                                         options.unit_multipliers)
        output('frequency range with applied unit multipliers:', rng)

    pb, wdir, bzone, mtxs = assemble_matrices(define, mod, pars, set_wave_dir,
                                              options)
    dim = pb.domain.shape.dim

    if dim != 2:
        options.plane = 'strain'

    if options.save_regions:
        pb.save_regions_as_groups(os.path.join(output_dir, 'regions'))

    if options.save_materials:
        save_materials(output_dir, pb, options)

    conf = pb.solver_confs['eig']
    eig_solver = Solver.any_from_conf(conf)

    n_eigs, options.n_eigs = setup_n_eigs(options, pb, mtxs)

    get_color = lambda ii: plt.cm.viridis((float(ii) / (options.n_eigs - 1)))
    plot_kwargs = [{'color' : get_color(ii), 'ls' : '', 'marker' : 'o'}
                  for ii in range(options.n_eigs)]

    log_names = []
    log_plot_kwargs = []
    if options.log_std_waves:
        std_wave_fun, log_names, log_plot_kwargs = get_std_wave_fun(
            pb, options)

    else:
        std_wave_fun = None

    stepper = get_stepper(rng, pb, options)

    if options.mode == 'omega':
        eigenshapes_filename = os.path.join(output_dir,
                                            'frequency-eigenshapes-%s.vtk'
                                            % stepper.suffix)

        log = Log([[r'$\lambda_{%d}$' % ii for ii in range(options.n_eigs)],
                   [r'$\omega_{%d}$'
                    % ii for ii in range(options.n_eigs)] + log_names],
                  plot_kwargs=[plot_kwargs, plot_kwargs + log_plot_kwargs],
                  formats=[['{:.5e}'] * options.n_eigs,
                           ['{:.5e}'] * (options.n_eigs + len(log_names))],
                  yscales=['linear', 'linear'],
                  xlabels=[r'$\kappa$', r'$\kappa$'],
                  ylabels=[r'eigenvalues $\lambda_i$',