operator_spectrum.py 4.93 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 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
Philipp Arras's avatar
Philipp Arras committed
21
from .field import Field
22
from .multi_domain import MultiDomain
Philipp Arras's avatar
Philipp Arras committed
23
from .multi_field import MultiField
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
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):
Philipp Arras's avatar
Philipp Arras committed
57
58
        self._check_input(x, mode)
        self._check_float_dtype(x)
59
60
61
62
63
        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
64
            res = np.empty(self.target.shape) if mode == self.TIMES else {}
65
66
67
68
69
70
71
72
            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)

Philipp Arras's avatar
Philipp Arras committed
73
74
75
76
77
78
79
80
81
82
83
84
    @staticmethod
    def _check_float_dtype(fld):
        if isinstance(fld, MultiField):
            dts = [ff.local_data.dtype for ff in fld.values()]
        elif isinstance(fld, Field):
            dts = [fld.local_data.dtype]
        else:
            raise TypeError
        for dt in dts:
            if not np.issubdtype(dt, np.float64):
                raise TypeError('Operator supports only floating point dtypes')

85

86
87
88
89
90
91
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
92
    A : LinearOperator
93
94
95
96
97
98
99
100
101
        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
102

103
            'LM' : largest magnitude
Philipp Arras's avatar
Philipp Arras committed
104

105
            'SM' : smallest magnitude
Philipp Arras's avatar
Philipp Arras committed
106

107
            'LR' : largest real part
Philipp Arras's avatar
Philipp Arras committed
108

109
            'SR' : smallest real part
Philipp Arras's avatar
Philipp Arras committed
110

111
            'LI' : largest imaginary part
Philipp Arras's avatar
Philipp Arras committed
112

113
            'SI' : smallest imaginary part
Philipp Arras's avatar
Philipp Arras committed
114

115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
    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):
133
        raise TypeError('Operator needs to be linear.')
134
    if A.domain is not A.target:
135
        raise TypeError('Operator needs to be endomorphism.')
136
    size = A.domain.size
Philipp Arras's avatar
Philipp Arras committed
137
138
139
140
    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())
141
    f = ssl.eigsh if hermitian else ssl.eigs
142
    eigs = f(M, k=k, tol=tol, return_eigenvectors=False, which=which)
143
    return np.flip(np.sort(eigs), axis=0)