eigenanalysis.py 文件源码

def eigs(As):
    """Solve eigenvalues from state matrix"""

    # if (not SLEPC) or (not SLEPC):
    return numpy.linalg.eig(As)
    # try:

    As = sparse(As)
    As = scipy.sparse.csc_matrix((numpy.array(As.CCS[2]).flatten(),
                                      numpy.array(As.CCS[1]).flatten(),
                                      numpy.array(As.CCS[0]).flatten()
                                      ), shape=As.size)
    # return linalg.eigs(As, k=1326)
    As_csr = As.tocsr()

    opts = PETSc.Options()
    A = PETSc.Mat().createAIJ(size=As_csr.shape,
                              csr=(As_csr.indptr, As_csr.indices, As_csr.data))
    A.setFromOptions()
    A.setUp()

    xr, tmp = A.getVecs()
    xi, tmp = A.getVecs()

    # Setup the eigensolver
    E = SLEPc.EPS().create()
    E.setOperators(A, None)
    E.setDimensions(500, PETSc.DECIDE)
    E.setProblemType( SLEPc.EPS.ProblemType.GNHEP )
    E.setFromOptions()

    # Solve the eigensystem
    E.solve()

    Print = PETSc.Sys.Print
    Print("")
    its = E.getIterationNumber()
    Print("Number of iterations of the method: %i" % its)
    sol_type = E.getType()
    Print("Solution method: %s" % sol_type)
    nev, ncv, mpd = E.getDimensions()
    Print("Number of requested eigenvalues: %i" % nev)
    tol, maxit = E.getTolerances()
    Print("Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit))
    nconv = E.getConverged()
    Print("Number of converged eigenpairs: %d" % nconv)
    if nconv > 0:
        Print("")
        Print("        k          ||Ax-kx||/||kx|| ")
        Print("----------------- ------------------")
        for i in range(nconv):
            k = E.getEigenpair(i, xr, xi)
            error = E.computeError(i)
            if k.imag != 0.0:
              Print(" %9f%+9f j  %12g" % (k.real, k.imag, error))
            else:
              Print(" %12f       %12g" % (k.real, error))
        Print("")
    # except:
    #     return numpy.linalg.eigvals(As)