def save_results(self, eigs, mtx_phi, out=None,
                     mesh_results_name=None, eig_results_name=None):
        mesh_results_name = get_default(mesh_results_name,
                                        self.mesh_results_name)
        eig_results_name = get_default(eig_results_name,
                                       self.eig_results_name)
        pb = self.problem

        save = self.app_options.save_eig_vectors
        n_eigs = self.app_options.n_eigs
        out = get_default(out, {})
        state = pb.create_state()
        aux = {}
        for ii in range(eigs.shape[0]):
            if save is not None:
                if (ii > save[0]) and (ii < (n_eigs - save[1])): continue
            state.set_full(mtx_phi[:,ii])
            aux = state.create_output_dict()
            key = list(aux.keys())[0]
            out[key+'%03d' % ii] = aux[key]

        if aux.get('__mesh__') is not None:
            out['__mesh__'] = aux['__mesh__']

        pb.save_state(mesh_results_name, out=out)

        fd = open(eig_results_name, 'w')
        eigs.tofile(fd, ' ')
        fd.close()

    def get(self, variable, quantity_name, bf=None, integration=None, step=None, time_derivative=None):
        """
        Get the named quantity related to the variable.

        Notes
        -----
        This is a convenience wrapper of Variable.evaluate() that
        initializes the arguments using the term data.
        """
        name = variable.name

        step = get_default(step, self.arg_steps[name])
        time_derivative = get_default(time_derivative, self.arg_derivatives[name])
        integration = get_default(integration, self.geometry_types[name])

        data = variable.evaluate(
            mode=quantity_name,
            region=self.region,
            integral=self.integral,
            integration=integration,
            step=step,
            time_derivative=time_derivative,
            is_trace=self.arg_traces[name],
            bf=bf,
        )
        return data

    def __call__(self, x0, conf=None, obj_fun=None, obj_fun_grad=None,
                 status=None, obj_args=None):
        import inspect

        if conf is not None:
            self.set_method(conf)

        else:
            conf = self.conf

        obj_fun = get_default(obj_fun, self.obj_fun)
        obj_fun_grad = get_default(obj_fun_grad, self.obj_fun_grad)
        status = get_default(status, self.status)
        obj_args = get_default(obj_args, self.obj_args)

        tt = time.clock()

        kwargs = {self._i_max_name[conf.method] : conf.i_max,
                  'args' : obj_args}

        if conf.method in self._has_grad:
            kwargs['fprime'] = obj_fun_grad

        if 'disp' in inspect.getargspec(self.solver)[0]:
            kwargs['disp'] = conf.verbose

        kwargs.update(self.build_solver_kwargs(conf))

        out = self.solver(obj_fun, x0, **kwargs)

        if status is not None:
            status['time_stats'] = time.clock() - tt

        return out

    def __call__(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                 i_max=None, mtx=None, status=None, **kwargs):

        eps_a = get_default(eps_a, self.conf.eps_a)
        eps_r = get_default(eps_r, self.conf.eps_r)
        i_max = get_default(i_max, self.conf.i_max)
        eps_d = self.conf.eps_d

        # There is no use in caching matrix in the solver - always set as new.
        pmtx, psol, prhs = self.set_matrix(mtx)

        ksp = self.ksp
        ksp.setOperators(pmtx)
        ksp.setFromOptions() # PETSc.Options() not used yet...
        ksp.setTolerances(atol=eps_a, rtol=eps_r, divtol=eps_d, max_it=i_max)

        # Set PETSc rhs, solve, get solution from PETSc solution.
        if x0 is not None:
            psol[...] = x0
            ksp.setInitialGuessNonzero(True)
        prhs[...] = rhs
        ksp.solve(prhs, psol)
        sol = psol[...].copy()
        output('%s(%s) convergence: %s (%s)'
               % (self.conf.method, self.conf.precond,
                  ksp.reason, self.converged_reasons[ksp.reason]))

        return sol

    def __call__(self, vec_x0, conf=None, fun=None, fun_grad=None,
                 lin_solver=None, iter_hook=None, status=None):
        if conf is not None:
            self.set_method( conf )
        else:
            conf = self.conf
        fun = get_default( fun, self.fun )
        status = get_default( status, self.status )

        tt = time.clock()

        kwargs = {'iter' : conf.i_max,
                  'alpha' : conf.alpha,
                  'verbose' : conf.verbose}

        if conf.method == 'broyden_generalized':
            kwargs.update( {'M' : conf.M} )

        elif conf.method in ['anderson', 'anderson2']:
            kwargs.update( {'M' : conf.M, 'w0' : conf.w0} )

        vec_x = self.solver( fun, vec_x0, **kwargs )
        vec_x = nm.asarray(vec_x)

        if status is not None:
            status['time_stats'] = time.clock() - tt

        return vec_x

    def solve_pressure_eigenproblem(self, mtx, eig_problem=None,
                                    n_eigs=0, check=False):
        """G = B*AI*BT or B*AI*BT+D"""

        def get_slice(n_eigs, nn):
            if n_eigs > 0:
                ii = slice(0, n_eigs)
            elif n_eigs < 0:
                ii = slice(nn + n_eigs, nn)
            else:
                ii = slice(0, 0)
            return ii

        eig_problem = get_default(eig_problem, self.eig_problem)
        n_eigs = get_default(n_eigs, self.n_eigs)
        check = get_default(check, self.check)

        mtx_c, mtx_b, action_aibt = mtx['C'], mtx['B'], mtx['action_aibt']
        mtx_g = mtx_b * action_aibt.to_array() # mtx_b must be sparse!
        if eig_problem == 'B*AI*BT+D':
            mtx_g += mtx['D'].toarray()

        mtx['G'] = mtx_g
        output(mtx_c.shape, mtx_g.shape)

        eigs, mtx_q = eig(mtx_c.toarray(), mtx_g, method='eig.sgscipy')

        if check:
            ee = nm.diag(sc.dot(mtx_q.T * mtx_c, mtx_q)).squeeze()
            oo = nm.diag(sc.dot(sc.dot(mtx_q.T,  mtx_g), mtx_q)).squeeze()
            try:
                assert_(nm.allclose(ee, eigs))
                assert_(nm.allclose(oo, nm.ones_like(eigs)))
            except ValueError:
                debug()

        nn = mtx_c.shape[0]
        if isinstance(n_eigs, tuple):
            output('required number of eigenvalues: (%d, %d)' % n_eigs)
            if sum(n_eigs) < nn:
                ii0 = get_slice(n_eigs[0], nn)
                ii1 = get_slice(-n_eigs[1], nn)
                eigs = nm.concatenate((eigs[ii0], eigs[ii1]))
                mtx_q = nm.concatenate((mtx_q[:,ii0], mtx_q[:,ii1]), 1) 
        else:
            output('required number of eigenvalues: %d' % n_eigs)
            if (n_eigs != 0) and (abs(n_eigs) < nn):
                ii = get_slice(n_eigs, nn)
                eigs = eigs[ii]
                mtx_q = mtx_q[:,ii]

##         from sfepy.base.plotutils import pylab, iplot
##         pylab.semilogy(eigs)
##         pylab.figure(2)
##         iplot(eigs)
##         pylab.show()
##         debug()

        out = Struct(eigs=eigs, mtx_q=mtx_q)
        return out

    def setup_ic( self, conf_ics = None, functions = None ):
        conf_ics = get_default(conf_ics, self.conf.ics)
        ics = Conditions.from_conf(conf_ics, self.domain.regions)

        functions = get_default(functions, self.functions)

        self.equations.setup_initial_conditions(ics, functions)

    def _standard_call(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                       i_max=None, mtx=None, status=None, context=None,
                       **kwargs):
        tt = time.clock()

        conf = get_default(conf, self.conf)
        mtx = get_default(mtx, self.mtx)
        status = get_default(status, self.status)
        context = get_default(context, self.context)

        assert_(mtx.shape[0] == mtx.shape[1] == rhs.shape[0])
        if x0 is not None:
            assert_(x0.shape[0] == rhs.shape[0])

        result = call(self, rhs, x0, conf, eps_a, eps_r, i_max, mtx, status,
                      context=context, **kwargs)
        if isinstance(result, tuple):
            result, n_iter = result

        else:
            n_iter = -1 # Number of iterations is undefined/unavailable.

        ttt = time.clock() - tt
        if status is not None:
            status['time'] = ttt
            status['n_iter'] = n_iter

        return result

    def _petsc_call(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                    i_max=None, mtx=None, status=None, comm=None,
                    context=None, **kwargs):
        tt = time.clock()

        conf = get_default(conf, self.conf)
        mtx = get_default(mtx, self.mtx)
        status = get_default(status, self.status)
        context = get_default(context, self.context)
        comm = get_default(comm, self.comm)

        mshape = mtx.size if isinstance(mtx, self.petsc.Mat) else mtx.shape
        rshape = [rhs.size] if isinstance(rhs, self.petsc.Vec) else rhs.shape

        assert_(mshape[0] == mshape[1] == rshape[0])
        if x0 is not None:
            xshape = [x0.size] if isinstance(x0, self.petsc.Vec) else x0.shape
            assert_(xshape[0] == rshape[0])

        result = call(self, rhs, x0, conf, eps_a, eps_r, i_max, mtx, status,
                      comm, context=context, **kwargs)

        ttt = time.clock() - tt
        if status is not None:
            status['time'] = ttt
            status['n_iter'] = self.ksp.getIterationNumber()

        return result

    def __call__(self, x0, conf=None, obj_fun=None, obj_fun_grad=None,
                 status=None, obj_args=None):

        if conf is not None:
            self.set_method(conf)

        else:
            conf = self.conf

        obj_fun = get_default(obj_fun, self.obj_fun)
        obj_fun_grad = get_default(obj_fun_grad, self.obj_fun_grad)
        status = get_default(status, self.status)
        obj_args = get_default(obj_args, self.obj_args)

        tt = time.clock()

        kwargs = {self._i_max_name[conf.method] : conf.i_max,
                  'disp' : conf.verbose,
                  'args' : obj_args}

        if conf.method in self._has_grad:
            kwargs['fprime'] = obj_fun_grad

        for key, val in conf.to_dict().iteritems():
            if key not in self._omit:
                kwargs[key] = val

        out = self.solver(obj_fun, x0, **kwargs)

        if status is not None:
            status['time_stats'] = time.clock() - tt

        return out

    def setup_ics(self, ics=None, functions=None):
        """
        Setup the initial conditions for use.
        """
        self.set_ics(get_default(ics, self.ics))

        functions = get_default(functions, self.functions)
        self.equations.setup_initial_conditions(self.ics, functions)

 def time_update(self, ts=None,
                 ebcs=None, epbcs=None, lcbcs=None,
                 functions=None, create_matrix=False):
     self.set_bcs(get_default(ebcs, self.ebcs),
                  get_default(epbcs, self.epbcs),
                  get_default(lcbcs, self.lcbcs))
     self.update_equations(ts, self.ebcs, self.epbcs, self.lcbcs,
                           functions, create_matrix)

    def __init__(self, filename, approx, region_selects, mat_pars, options,
                 evp_options, eigenmomenta_options, band_gaps_options,
                 coefs_save_name='coefs',
                 corrs_save_names=None,
                 incwd=None,
                 output_dir=None, **kwargs):
        Struct.__init__(self, approx=approx, region_selects=region_selects,
                        mat_pars=mat_pars, options=options,
                        evp_options=evp_options,
                        eigenmomenta_options=eigenmomenta_options,
                        band_gaps_options=band_gaps_options,
                        **kwargs)
        self.incwd = get_default(incwd, lambda x: x)

        self.conf = Struct()
        self.conf.filename_mesh = self.incwd(filename)

        output_dir = get_default(output_dir, self.incwd('output'))

        default = {'evp' : 'evp', 'corrs_rs' : 'corrs_rs'}
        self.corrs_save_names = get_default(corrs_save_names,
                                            default)

        io = MeshIO.any_from_filename(self.conf.filename_mesh)
        self.bbox, self.dim = io.read_bounding_box(ret_dim=True)
        rpc_axes = nm.eye(self.dim, dtype=nm.float64) \
                   * (self.bbox[1] - self.bbox[0])

        self.conf.options = options
        self.conf.options.update({
            'output_dir' : output_dir,

            'volume' : {
                'value' : get_lattice_volume(rpc_axes),
            },

            'coefs' : 'coefs',
            'requirements' : 'requirements',

            'coefs_filename' : coefs_save_name,
        })

        self.conf.mat_pars = mat_pars

        self.conf.solvers = self.define_solvers()
        self.conf.regions = self.define_regions()
        self.conf.materials = self.define_materials()
        self.conf.fields = self.define_fields()
        self.conf.variables = self.define_variables()
        (self.conf.ebcs, self.conf.epbcs,
         self.conf.lcbcs, self.all_periodic) = self.define_bcs()
        self.conf.functions = self.define_functions()
        self.conf.integrals = self.define_integrals()

        self.equations, self.expr_coefs = self.define_equations()
        self.conf.coefs = self.define_coefs()
        self.conf.requirements = self.define_requirements()

    def __init__(self, data_names=None, yscales=None,
                 xlabels=None, ylabels=None, is_plot=True, aggregate=200,
                 formats=None, log_filename=None):
        """`data_names` ... tuple of names grouped by subplots:
                            ([name1, name2, ...], [name3, name4, ...], ...)
        where name<n> are strings to display in (sub)plot legends."""
        try:
            import matplotlib as mpl
        except:
            mpl = None

        if (mpl is not None) and mpl.rcParams['backend'] == 'GTKAgg':
            can_live_plot = True
        else:
            can_live_plot = False

        Struct.__init__(self, data_names = {},
                        n_arg = 0, n_gr = 0,
                        data = {}, x_values = {}, n_calls = 0,
                        yscales = {}, xlabels = {}, ylabels = {},
                        plot_pipe = None, formats = {}, output = None)

        if data_names is not None:
            n_gr = len(data_names)
        else:
            n_gr = 0
            data_names = []

        yscales = get_default(yscales, ['linear'] * n_gr)
        xlabels = get_default(xlabels, ['iteration'] * n_gr)
        ylabels = get_default(ylabels, [''] * n_gr )

        if formats is None:
            formats = [None] * n_gr

        for ig, names in enumerate(data_names):
            self.add_group(names, yscales[ig], xlabels[ig], ylabels[ig],
                           formats[ig])

        self.is_plot = get_default( is_plot, True )
        self.aggregate = get_default( aggregate, 100 )

        self.can_plot = (can_live_plot and (mpl is not None)
                         and (Process is not None))

        if log_filename is not None:
            self.output = Output('', filename=log_filename)
            self.output('# started: %s' % time.asctime())
            self.output('# groups: %d' % n_gr)
            for ig, names in enumerate(data_names):
                self.output('#   %d' % ig)
                self.output('#     xlabel: "%s", ylabel: "%s", yscales: "%s"'
                            % (xlabels[ig], ylabels[ig], yscales[ig]))
                self.output('#     names: "%s"' % ', '.join(names))

        if self.is_plot and (not self.can_plot):
            output(_msg_no_live)

    def __call__(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                 i_max=None, mtx=None, status=None, context=None, **kwargs):
        solver_kwargs = self.build_solver_kwargs(conf)

        eps_a = get_default(eps_a, self.conf.eps_a)
        eps_r = get_default(eps_r, self.conf.eps_r)
        i_max = get_default(i_max, self.conf.i_max)

        setup_precond = get_default(kwargs.get('setup_precond', None),
                                    self.conf.setup_precond)
        callback = get_default(kwargs.get('callback', lambda sol: None),
                               self.conf.callback)

        self.iter = 0
        def iter_callback(sol):
            self.iter += 1
            msg = '%s: iteration %d' % (self.conf.name, self.iter)
            if conf.verbose > 2:
                if conf.method not in self._callbacks_res:
                    res = mtx * sol - rhs

                else:
                    res = sol

                rnorm = nm.linalg.norm(res)
                msg += ': |Ax-b| = %e' % rnorm
            output(msg, verbose=conf.verbose > 1)

            # Call an optional user-defined callback.
            callback(sol)

        precond = setup_precond(mtx, context)

        if conf.method == 'qmr':
            prec_args = {'M1' : precond, 'M2' : precond}

        else:
            prec_args = {'M' : precond}

        solver_kwargs.update(prec_args)

        try:
            sol, info = self.solver(mtx, rhs, x0=x0, atol=eps_a, tol=eps_r,
                                    maxiter=i_max, callback=iter_callback,
                                    **solver_kwargs)
        except TypeError:
            sol, info = self.solver(mtx, rhs, x0=x0, tol=eps_r,
                                    maxiter=i_max, callback=iter_callback,
                                    **solver_kwargs)

        output('%s: %s convergence: %s (%s, %d iterations)'
               % (self.conf.name, self.conf.method,
                  info, self.converged_reasons[nm.sign(info)], self.iter),
               verbose=conf.verbose)

        return sol, self.iter

    def __call__(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                 i_max=None, mtx=None, status=None, **kwargs):
        conf = get_default(conf, self.conf)
        mtx = get_default(mtx, self.mtx)
        status = get_default(status, self.status)

        if self.solve is not None:
            # Matrix is already prefactorized.
            return self.solve(rhs)
        else:
            return self.sls.spsolve(mtx, rhs)

    def __call__( self, state0 = None, conf = None,
                  step_fun = None, step_args = None ):

        step_fun = get_default( step_fun, self.step_fun )
        step_args = get_default( step_args, self.step_args )

        for step, time in self.ts:
            output( self.format % (time, step + 1, self.ts.n_step) )

            state = step_fun( self.ts, state0, *step_args )
            state0 = state.copy(deep=True)
            yield self.ts, state

    def __call__(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                 i_max=None, mtx=None, status=None, **kwargs):

        eps_r = get_default(eps_r, self.conf.eps_r)
        i_max = get_default(i_max, self.conf.i_max)

        if (self.mg is None) or (mtx is not self.mtx):
            self.mg = self.solver(mtx)
            self.mtx = mtx

        sol = self.mg.solve(rhs, x0=x0, accel=conf.accel, tol=eps_r,
                            maxiter=i_max)

        return sol

def solve(mtx, rhs, solver_class=None, solver_conf=None):
    """
    Solve the linear system with the matrix `mtx` and the right-hand side
    `rhs`.

    Convenience wrapper around the linear solver classes below.
    """
    solver_class = get_default(solver_class, ScipyDirect)
    solver_conf = get_default(solver_conf, {})

    solver = solver_class(solver_conf, mtx=mtx)
    solution = solver(rhs)

    return solution

    def __call__(self, vec0=None, nls=None, init_fun=None, prestep_fun=None,
                 poststep_fun=None, status=None, **kwargs):
        """
        Solve elastodynamics problems by the generalized :math:`\alpha` method.
        """
        conf = self.conf
        nls = get_default(nls, self.nls)

        rho_inf = conf.rho_inf
        alpha_m = get_default(conf.alpha_m,
                              (2.0 * rho_inf - 1.0) / (rho_inf + 1.0))
        alpha_f = get_default(conf.alpha_f, rho_inf / (rho_inf + 1.0))
        beta = get_default(conf.beta, 0.25 * (1.0 - alpha_m + alpha_f)**2)
        gamma = get_default(conf.gamma, 0.5 - alpha_m + alpha_f)

        output('parameters rho_inf, alpha_m, alpha_f, beta, gamma:',
               verbose=self.verbose)
        output(rho_inf, alpha_m, alpha_f, beta, gamma,
               verbose=self.verbose)

        vec, unpack, pack = self.get_initial_vec(
            nls, vec0, init_fun, prestep_fun, poststep_fun)

        ts = self.ts
        for step, time in ts.iter_from(ts.step):
            output(self.format % (time, step + 1, ts.n_step),
                   verbose=self.verbose)
            dt = ts.dt

            prestep_fun(ts, vec)
            ut, vt, at = unpack(vec)

            nlst = self.create_nlst(nls, dt, alpha_m, alpha_f, gamma, beta,
                                    ut, vt, at)

            ts.set_substep_time((1.0 - alpha_f) * dt)
            am = nlst(at)
            ts.restore_step_time()

            atp = nlst.a1(am)
            vtp = nlst.v1(atp)
            utp = nlst.u1(atp)

            vect = pack(utp, vtp, atp)
            poststep_fun(ts, vect)

            vec = vect

        return vec