operator_spectrum.py 4.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society

import numpy as np
import scipy.sparse.linalg as ssl

from .domain_tuple import DomainTuple
from .domains.unstructured_domain import UnstructuredDomain
from .multi_domain import MultiDomain
from .operators.linear_operator import LinearOperator
from .operators.sandwich_operator import SandwichOperator
from .sugar import from_global_data, makeDomain


class _DomRemover(LinearOperator):
    """Operator which transforms between a structured MultiDomain
    and an unstructured domain.

    Parameters
    ----------
    domain: MultiDomain
        the full input domain of the operator.

    Notes
    -----
    The operator converts the full domain of its input domain to an
    UnstructuredDomain
    """

    def __init__(self, domain):
        self._domain = makeDomain(domain)
        if isinstance(self._domain, MultiDomain):
            self._size_array = np.array([0] +
                                        [d.size for d in domain.values()])
        else:
            self._size_array = np.array([0, domain.size])
        np.cumsum(self._size_array, out=self._size_array)
        target = UnstructuredDomain(self._size_array[-1])
        self._target = makeDomain(target)
        self._capability = self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
        x = x.to_global_data()
        if isinstance(self._domain, DomainTuple):
            res = x.ravel() if mode == self.TIMES else x.reshape(
                self._domain.shape)
        else:
Philipp Arras's avatar
Philipp Arras committed
60
            res = np.empty(self.target.shape, dtype=x.dtype) if mode == self.TIMES else {}
61 62 63 64 65 66 67 68 69
            for ii, (kk, dd) in enumerate(self.domain.items()):
                i0, i1 = self._size_array[ii:ii + 2]
                if mode == self.TIMES:
                    res[i0:i1] = x[kk].ravel()
                else:
                    res[kk] = x[i0:i1].reshape(dd.shape)
        return from_global_data(self._tgt(mode), res)


70 71 72 73 74 75
def operator_spectrum(A, k, hermitian, which='LM', tol=0):
    '''
    Find k eigenvalues and eigenvectors of the endomorphism A.

    Parameters
    ----------
Philipp Arras's avatar
Philipp Arras committed
76
    A : LinearOperator
77 78 79 80 81 82 83 84 85
        Operator of which eigenvalues shall be computed.
    k : int
        The number of eigenvalues and eigenvectors desired. `k` must be
        smaller than N-1. It is not possible to compute all eigenvectors of a
        matrix.
    hermitian: bool
        Specifies whether A is hermitian or not.
    which : str, ['LM' | 'SM' | 'LR' | 'SR' | 'LI' | 'SI'], optional
        Which `k` eigenvectors and eigenvalues to find:
Philipp Arras's avatar
Philipp Arras committed
86

87
            'LM' : largest magnitude
Philipp Arras's avatar
Philipp Arras committed
88

89
            'SM' : smallest magnitude
Philipp Arras's avatar
Philipp Arras committed
90

91
            'LR' : largest real part
Philipp Arras's avatar
Philipp Arras committed
92

93
            'SR' : smallest real part
Philipp Arras's avatar
Philipp Arras committed
94

95
            'LI' : largest imaginary part
Philipp Arras's avatar
Philipp Arras committed
96

97
            'SI' : smallest imaginary part
Philipp Arras's avatar
Philipp Arras committed
98

99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
    tol : float, optional
        Relative accuracy for eigenvalues (stopping criterion)
        The default value of 0 implies machine precision.

    Returns
    -------
    w : ndarray
        Array of k eigenvalues.

    Raises
    ------
    ArpackNoConvergence
        When the requested convergence is not obtained.
        The currently converged eigenvalues and eigenvectors can be found
        as ``eigenvalues`` and ``eigenvectors`` attributes of the exception
        object.
    '''
    if not isinstance(A, LinearOperator):
117
        raise TypeError('Operator needs to be linear.')
118
    if A.domain is not A.target:
119
        raise TypeError('Operator needs to be endomorphism.')
120
    size = A.domain.size
Philipp Arras's avatar
Philipp Arras committed
121 122 123 124
    Ar = SandwichOperator.make(_DomRemover(A.domain).adjoint, A)
    M = ssl.LinearOperator(
        shape=2*(size,),
        matvec=lambda x: Ar(from_global_data(Ar.domain, x)).to_global_data())
125
    f = ssl.eigsh if hermitian else ssl.eigs
126
    eigs = f(M, k=k, tol=tol, return_eigenvectors=False, which=which)
127
    return np.flip(np.sort(eigs), axis=0)