Python solvers.eig函数代码示例

    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 __call__(self, problem=None, data=None):
        problem = get_default(problem, self.problem)
        opts = self.app_options

        problem.set_equations(self.equations)
        problem.select_bcs(ebc_names=self.ebcs, epbc_names=self.epbcs,
                           lcbc_names=self.get('lcbcs', []))
        problem.update_materials(problem.ts)

        self.init_solvers(problem)

        mtx_a, mtx_m, data = self.prepare_matrices(problem)

        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=opts.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:
            from sfepy.base.base import debug; debug()

        mtx_phi, eig_vectors = self.post_process(eigs, mtx_s_phi, data,
                                                 problem)

        self.save(eigs, mtx_phi, problem)

        evp = Struct(name='evp', eigs=eigs, eigs_rescaled=eigs_rescaled,
                     eig_vectors=eig_vectors)

        return evp

    def compute_phase_velocity( self ):
        from sfepy.homogenization.phono import compute_density_volume_info
        opts = self.app_options
        dim = self.problem.domain.mesh.dim

        christoffel = self.compute_cat()

        self.problem.update_materials()
        dv_info = compute_density_volume_info( self.problem, opts.volume,
                                               opts.region_to_material )
        output( 'average density:', dv_info.average_density )

        eye = nm.eye( dim, dim, dtype = nm.float64 )
        mtx_mass = eye * dv_info.average_density

        meigs, mvecs = eig( mtx_mass, mtx_b = christoffel,
                            eigenvectors = True, method = opts.eigensolver )
        phase_velocity = 1.0 / nm.sqrt( meigs )

        return phase_velocity

def main():
    parser = OptionParser(usage=usage, version="%prog")
    parser.add_option(
        "-b", "--basis", metavar="name", action="store", dest="basis", default="lagrange", help=help["basis"]
    )
    parser.add_option(
        "-n",
        "--max-order",
        metavar="order",
        type=int,
        action="store",
        dest="max_order",
        default=10,
        help=help["max_order"],
    )
    parser.add_option(
        "-m",
        "--matrix",
        metavar="type",
        action="store",
        dest="matrix_type",
        default="laplace",
        help=help["matrix_type"],
    )
    parser.add_option(
        "-g", "--geometry", metavar="name", action="store", dest="geometry", default="2_4", help=help["geometry"]
    )
    options, args = parser.parse_args()

    dim, n_ep = int(options.geometry[0]), int(options.geometry[2])
    output("reference element geometry:")
    output("  dimension: %d, vertices: %d" % (dim, n_ep))

    n_c = {"laplace": 1, "elasticity": dim}[options.matrix_type]

    output("matrix type:", options.matrix_type)
    output("number of variable components:", n_c)

    output("polynomial space:", options.basis)

    output("max. order:", options.max_order)

    mesh = Mesh.from_file(data_dir + "/meshes/elements/%s_1.mesh" % options.geometry)
    domain = Domain("domain", mesh)
    omega = domain.create_region("Omega", "all")

    orders = nm.arange(1, options.max_order + 1, dtype=nm.int)
    conds = []

    order_fix = 0 if options.geometry in ["2_4", "3_8"] else 1

    for order in orders:
        output("order:", order, "...")

        field = Field.from_args(
            "fu", nm.float64, n_c, omega, approx_order=order, space="H1", poly_space_base=options.basis
        )

        to = field.approx_order
        quad_order = 2 * (max(to - order_fix, 0))
        output("quadrature order:", quad_order)

        integral = Integral("i", order=quad_order)
        qp, _ = integral.get_qp(options.geometry)
        output("number of quadrature points:", qp.shape[0])

        u = FieldVariable("u", "unknown", field, n_c)
        v = FieldVariable("v", "test", field, n_c, primary_var_name="u")

        m = Material("m", lam=1.0, mu=1.0)

        if options.matrix_type == "laplace":
            term = Term.new("dw_laplace(m.mu, v, u)", integral, omega, m=m, v=v, u=u)
            n_zero = 1

        else:
            assert_(options.matrix_type == "elasticity")
            term = Term.new("dw_lin_elastic_iso(m.lam, m.mu, v, u)", integral, omega, m=m, v=v, u=u)
            n_zero = (dim + 1) * dim / 2

        term.setup()

        output("assembling...")
        tt = time.clock()
        mtx, iels = term.evaluate(mode="weak", diff_var="u")
        output("...done in %.2f s" % (time.clock() - tt))
        mtx = mtx[0][0, 0]

        try:
            assert_(nm.max(nm.abs(mtx - mtx.T)) < 1e-10)

        except:
            from sfepy.base.base import debug

            debug()

        output("matrix shape:", mtx.shape)

        eigs = eig(mtx, method="eig.sgscipy", eigenvectors=False)
        eigs.sort()
#.........这里部分代码省略.........

def main():
    parser = OptionParser(usage=usage, version='%prog')
    parser.add_option('-b', '--basis', metavar='name',
                      action='store', dest='basis',
                      default='lagrange', help=help['basis'])
    parser.add_option('-n', '--max-order', metavar='order', type=int,
                      action='store', dest='max_order',
                      default=10, help=help['max_order'])
    parser.add_option('-m', '--matrix', metavar='type',
                      action='store', dest='matrix_type',
                      default='laplace', help=help['matrix_type'])
    parser.add_option('-g', '--geometry', metavar='name',
                      action='store', dest='geometry',
                      default='2_4', help=help['geometry'])
    options, args = parser.parse_args()

    dim, n_ep = int(options.geometry[0]), int(options.geometry[2])
    output('reference element geometry:')
    output('  dimension: %d, vertices: %d' % (dim, n_ep))

    n_c = {'laplace' : 1, 'elasticity' : dim}[options.matrix_type]

    output('matrix type:', options.matrix_type)
    output('number of variable components:',  n_c)

    output('polynomial space:', options.basis)

    output('max. order:', options.max_order)

    mesh = Mesh.from_file(data_dir + '/meshes/elements/%s_1.mesh'
                          % options.geometry)
    domain = FEDomain('domain', mesh)
    omega = domain.create_region('Omega', 'all')

    orders = nm.arange(1, options.max_order + 1, dtype=nm.int)
    conds = []

    order_fix = 0 if  options.geometry in ['2_4', '3_8'] else 1

    for order in orders:
        output('order:', order, '...')

        field = Field.from_args('fu', nm.float64, n_c, omega,
                                approx_order=order,
                                space='H1', poly_space_base=options.basis)

        to = field.approx_order
        quad_order = 2 * (max(to - order_fix, 0))
        output('quadrature order:', quad_order)

        integral = Integral('i', order=quad_order)
        qp, _ = integral.get_qp(options.geometry)
        output('number of quadrature points:', qp.shape[0])

        u = FieldVariable('u', 'unknown', field)
        v = FieldVariable('v', 'test', field, primary_var_name='u')

        m = Material('m', lam=1.0, mu=1.0)

        if options.matrix_type == 'laplace':
            term = Term.new('dw_laplace(m.mu, v, u)',
                            integral, omega, m=m, v=v, u=u)
            n_zero = 1

        else:
            assert_(options.matrix_type == 'elasticity')
            term = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u)',
                            integral, omega, m=m, v=v, u=u)
            n_zero = (dim + 1) * dim / 2

        term.setup()

        output('assembling...')
        tt = time.clock()
        mtx, iels = term.evaluate(mode='weak', diff_var='u')
        output('...done in %.2f s' % (time.clock() - tt))
        mtx = mtx[0][0, 0]

        try:
            assert_(nm.max(nm.abs(mtx - mtx.T)) < 1e-10)

        except:
            from sfepy.base.base import debug; debug()

        output('matrix shape:', mtx.shape)

        eigs = eig(mtx, method='eig.sgscipy', eigenvectors=False)
        eigs.sort()

        # Zero 'true' zeros.
        eigs[:n_zero] = 0.0

        ii = nm.where(eigs < 0.0)[0]
        if len(ii):
            output('matrix is not positive semi-definite!')

        ii = nm.where(eigs[n_zero:] < 1e-12)[0]
        if len(ii):
            output('matrix has more than %d zero eigenvalues!' % n_zero)

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

    def trace_full_callback(f):
        meigs, mvecs = eig((f**2) * mass(f), mtx_b=mtx_b,
                           eigenvectors=True, method=method)

        return meigs, mvecs

 def trace_callback(f):
     meigs = eig(mass(f), eigenvectors=False, method=method)
     return meigs,

 def find_zero_full_callback(f):
     meigs = eig((f**2) * mass(f), mtx_b=mtx_b,
                 eigenvectors=False, method=method)
     return meigs

 def find_zero_callback(f):
     meigs = eig(mass(f), eigenvectors=False, method=method)
     return meigs

    def solve_eigen_problem_1( self ):
        from sfepy.fem import Mesh

        options = self.options
        opts = self.app_options
        pb = self.problem

        dim = pb.domain.mesh.dim

        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) )

        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) )

        n_eigs = get_default( opts.n_eigs, mtx_a.shape[0] )
##         mtx_a.save( 'a.txt', format='%d %d %.12f\n' )
##         mtx_b.save( 'b.txt', format='%d %d %.12f\n' )

        output( 'computing resonance frequencies...' )
        eig = Solver.any_from_conf( pb.get_solver_conf( opts.eigen_solver ) )
        eigs, mtx_s_phi = eig( mtx_a, mtx_b, n_eigs )
        output( '...done' )

        bounding_box = pb.domain.mesh.get_bounding_box()
        # this assumes a box (3D), or a square (2D):
        a = bounding_box[1][0] - bounding_box[0][0]
        E_exact = None
        if options.hydrogen or options.boron:
            if options.hydrogen:
                Z = 1
            elif options.boron:
                Z = 5
            if dim == 2:
                E_exact = [-float(Z)**2/2/(n-0.5)**2/4
                           for n in [1]+[2]*3+[3]*5 + [4]*8 + [5]*15]
            elif dim == 3:
                E_exact = [-float(Z)**2/2/n**2 for n in [1]+[2]*2**2+[3]*3**2 ]
        if options.well:
            if dim == 2:
                E_exact = [pi**2/(2*a**2)*x
                           for x in [2, 5, 5, 8, 10, 10, 13, 13,
                                     17, 17, 18, 20, 20 ] ]
            elif dim == 3:
                E_exact = [pi**2/(2*a**2)*x
                           for x in [3, 6, 6, 6, 9, 9, 9, 11, 11,
                                     11, 12, 14, 14, 14, 14, 14,
                                     14, 17, 17, 17] ]
        if options.oscillator:
            if dim == 2:
                E_exact = [1] + [2]*2 + [3]*3 + [4]*4 + [5]*5 + [6]*6
            elif dim == 3:
                E_exact = [float(1)/2+x for x in [1]+[2]*3+[3]*6+[4]*10 ]
        if E_exact is not None:
            output("a=%f" % a)
            output("Energies:")
            output("n      exact         FEM      error")

            for i, e in enumerate(eigs):
                from numpy import NaN
                if i < len(E_exact):
                    exact = E_exact[i]
                    err = 100*abs((exact - e)/exact)
                else:
                    exact = NaN
                    err = NaN
                output("%d:  %.8f   %.8f  %5.2f%%" % (i, exact, e, err))
        else:
            output(eigs)
##         import sfepy.base.plotutils as plu
##         plu.spy( mtx_b, eps = 1e-12 )
##         plu.plt.show()
##         pause()

        mtx_phi = self.make_full( mtx_s_phi )
        self.save_results( eigs, mtx_phi )

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

    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']
#.........这里部分代码省略.........