Skip to content

fixed_dof_solve

fixed_dof_solve_precompute

Source code in src/gpytoolbox/fixed_dof_solve.py 
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
class fixed_dof_solve_precompute:
    def __init__(self, A, k=None):
        """Prepare a precomputation object to efficiently solve a linear system
        while fixing certain degrees of freedom for which the linear system is to be
        ignored during solution.

        For the linear system `A*u = b`, the linear system will be enforced at all
        the rows that don't correspond to fixed degrees of freedom, and the fixed
        degrees of freedom will be enforced on their respective rows.

        This can be used to implement finite differences with Dirichlet boundary
        conditions by fixing the boundary degrees of freedom to the appropriate
        boundary values.

        Parameters
        ----------
        A : (n,n) scipy csc matrix
            square matrix for the linear system
        k : (o,) numpy int array
            index vector of fixed degrees of freedom


        Returns
        -------
        precomputed : precomputation object that can be used to solve the problem



        See Also
        --------
        min_quad_with_fixed


        Examples
        --------
        ```python
        >>> import gpytoolbox as gpy
        >>> import numpy as np
        >>> import scipy as sp
        >>> 
        >>> # This matrix is not symmetric!
        >>> A = sp.sparse.csc_matrix(np.array([[1.,0.,1.],[0.,1.,0.],[0.,0.,1.]]))
        >>> k = np.array([2])
        >>> precomp = gpy.fixed_dof_solve_precompute(A, k=k)
        >>> 
        >>> b = np.array([0.,0.,0.])
        >>> y = np.array([5.])
        >>> u = precomp.solve(b=b, y=y)
        >>> u
        array([-5.,  0.,  5.])
        >>> 
        >>> # A*u matches b, except for indices where constraints apply.
        >>> A*u - b
        array([0., 0., 5.])
        ```
        """

        self.n = A.shape[0]
        assert A.shape[1] == self.n
        assert self.n>0
        # self.A = A.copy()

        if k is None:
            self.o = 0
            self.k = None
        else:
            self.o = k.shape[0]
            assert k.shape == (self.o,)
            assert np.min(k)>=0 and np.max(k)<self.n
            assert np.unique(k).shape == k.shape, "No duplicate indices"
            self.k = k.copy()

        # If k is provided, remove these degrees of freedom.
        if self.o==0:
            self.ki = None
            self.n_reduced = self.n
            self.Ared = A
            self.A_for_extra_b = None
        else:
            self.ki = np.setdiff1d(np.arange(0, self.n), self.k)
            self.n_reduced = self.ki.shape[0]
            assert self.n_reduced == self.n - self.k.shape[0]
            self.Ared = A[self.ki, :][:, self.ki]
            self.A_for_extra_b = A[self.ki, :][:, self.k]

        splu = sp.sparse.linalg.splu(self.Ared)
        self.solver = lambda x : splu.solve(x)
    def solve(self, b=None, y=None):
        """Solve the prefactored linear system
        while fixing certain degrees of freedom for which the linear system is to be
        ignored during solution.

        For the linear system `A*u = b`, the linear system will be enforced at all
        the rows that don't correspond to fixed degrees of freedom, and the fixed
        degrees of freedom will be enforced on their respective rows.

        This can be used to implement finite differences with Dirichlet boundary
        conditions by fixing the boundary degrees of freedom to the appropriate
        boundary values.

        Parameters
        ----------
        b : (n,) or (n,p) numpy float array
            right-hand side of the linear system
        y : (o,) or (o,p) numpy float array
            what the degrees of freedom are fixed to, `u[k] == y` or `u[k,:] == y`


        Returns
        -------
        u : (n,) or (n,p) numpy float array such that
            `A[not k,:] * u == b[not k]` or `A[not k,:] * u == b[not k,:]`
            and `u[k] == y` or `u[k,:] == y`.


        See Also
        --------
        min_quad_with_fixed


        Examples
        --------
        ```python
        >>> import gpytoolbox as gpy
        >>> import numpy as np
        >>> import scipy as sp
        >>> 
        >>> # This matrix is not symmetric!
        >>> A = sp.sparse.csc_matrix(np.array([[1.,0.,1.],[0.,1.,0.],[0.,0.,1.]]))
        >>> k = np.array([2])
        >>> precomp = gpy.fixed_dof_solve_precompute(A, k=k)
        >>> 
        >>> b = np.array([0.,0.,0.])
        >>> y = np.array([5.])
        >>> u = precomp.solve(b=b, y=y)
        >>> u
        array([-5.,  0.,  5.])
        >>> 
        >>> # A*u matches b, except for indices where constraints apply.
        >>> A*u - b
        array([0., 0., 5.])
        ```

        """

        def cp(x):
            if x is None or np.isscalar(x):
                return 0
            if len(x.shape)==1:
                return 0
            return x.shape[1]
        p = max([cp(b), cp(y)])

        assert b is None or np.isscalar(b) or (p==0 and b.shape==(self.n,)) or (p>0 and b.shape==(self.n,p))
        assert y is None or (np.isscalar(y) and self.o>0) or (p==0 and y.shape==(self.o,)) or (p>0 and y.shape==(self.o,p))

        # Get everything to full dimensions
        if b is None:
            b = 0.
        if np.isscalar(b):
            if p==0:
                b = np.full(self.n, b)
            else:
                b = np.full((self.n,p), b)
        if y is None and self.o>0:
            y = 0.
        if np.isscalar(y) and self.o>0:
            if p==0:
                y = np.full(self.o, y)
            else:
                y = np.full((self.o,p), y)

        # Modified rhs based on known values
        if self.o==0:
            rhs = b
        else:
            b_reduced = b[self.ki] if p==0 else b[self.ki,:]
            rhs = b_reduced - self.A_for_extra_b*y
        ured = self.solver(rhs)

        if self.o==0:
            u = ured
        else:
            if p==0:
                u = np.empty(self.n, dtype=np.float64)
                u[self.ki] = ured
                u[self.k] = y
            else:
                u = np.empty((self.n,p), dtype=np.float64)
                u[self.ki,:] = ured
                u[self.k,:] = y

        return u

__init__(A, k=None)

Prepare a precomputation object to efficiently solve a linear system while fixing certain degrees of freedom for which the linear system is to be ignored during solution.

For the linear system A*u = b, the linear system will be enforced at all the rows that don't correspond to fixed degrees of freedom, and the fixed degrees of freedom will be enforced on their respective rows.

This can be used to implement finite differences with Dirichlet boundary conditions by fixing the boundary degrees of freedom to the appropriate boundary values.

Parameters:

Name Type Description Default
A (n,n) scipy csc matrix

square matrix for the linear system

 required
k (o,) numpy int array

index vector of fixed degrees of freedom

 None

Returns:

Name Type Description
precomputed precomputation object that can be used to solve the problem
See Also

min_quad_with_fixed

Examples:

>>> import gpytoolbox as gpy
>>> import numpy as np
>>> import scipy as sp
>>> 
>>> # This matrix is not symmetric!
>>> A = sp.sparse.csc_matrix(np.array([[1.,0.,1.],[0.,1.,0.],[0.,0.,1.]]))
>>> k = np.array([2])
>>> precomp = gpy.fixed_dof_solve_precompute(A, k=k)
>>> 
>>> b = np.array([0.,0.,0.])
>>> y = np.array([5.])
>>> u = precomp.solve(b=b, y=y)
>>> u
array([-5.,  0.,  5.])
>>> 
>>> # A*u matches b, except for indices where constraints apply.
>>> A*u - b
array([0., 0., 5.])
Source code in src/gpytoolbox/fixed_dof_solve.py 
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def __init__(self, A, k=None):
    """Prepare a precomputation object to efficiently solve a linear system
    while fixing certain degrees of freedom for which the linear system is to be
    ignored during solution.

    For the linear system `A*u = b`, the linear system will be enforced at all
    the rows that don't correspond to fixed degrees of freedom, and the fixed
    degrees of freedom will be enforced on their respective rows.

    This can be used to implement finite differences with Dirichlet boundary
    conditions by fixing the boundary degrees of freedom to the appropriate
    boundary values.

    Parameters
    ----------
    A : (n,n) scipy csc matrix
        square matrix for the linear system
    k : (o,) numpy int array
        index vector of fixed degrees of freedom


    Returns
    -------
    precomputed : precomputation object that can be used to solve the problem



    See Also
    --------
    min_quad_with_fixed


    Examples
    --------
    ```python
    >>> import gpytoolbox as gpy
    >>> import numpy as np
    >>> import scipy as sp
    >>> 
    >>> # This matrix is not symmetric!
    >>> A = sp.sparse.csc_matrix(np.array([[1.,0.,1.],[0.,1.,0.],[0.,0.,1.]]))
    >>> k = np.array([2])
    >>> precomp = gpy.fixed_dof_solve_precompute(A, k=k)
    >>> 
    >>> b = np.array([0.,0.,0.])
    >>> y = np.array([5.])
    >>> u = precomp.solve(b=b, y=y)
    >>> u
    array([-5.,  0.,  5.])
    >>> 
    >>> # A*u matches b, except for indices where constraints apply.
    >>> A*u - b
    array([0., 0., 5.])
    ```
    """

    self.n = A.shape[0]
    assert A.shape[1] == self.n
    assert self.n>0
    # self.A = A.copy()

    if k is None:
        self.o = 0
        self.k = None
    else:
        self.o = k.shape[0]
        assert k.shape == (self.o,)
        assert np.min(k)>=0 and np.max(k)<self.n
        assert np.unique(k).shape == k.shape, "No duplicate indices"
        self.k = k.copy()

    # If k is provided, remove these degrees of freedom.
    if self.o==0:
        self.ki = None
        self.n_reduced = self.n
        self.Ared = A
        self.A_for_extra_b = None
    else:
        self.ki = np.setdiff1d(np.arange(0, self.n), self.k)
        self.n_reduced = self.ki.shape[0]
        assert self.n_reduced == self.n - self.k.shape[0]
        self.Ared = A[self.ki, :][:, self.ki]
        self.A_for_extra_b = A[self.ki, :][:, self.k]

    splu = sp.sparse.linalg.splu(self.Ared)
    self.solver = lambda x : splu.solve(x)

solve(b=None, y=None)

Solve the prefactored linear system while fixing certain degrees of freedom for which the linear system is to be ignored during solution.

For the linear system A*u = b, the linear system will be enforced at all the rows that don't correspond to fixed degrees of freedom, and the fixed degrees of freedom will be enforced on their respective rows.

This can be used to implement finite differences with Dirichlet boundary conditions by fixing the boundary degrees of freedom to the appropriate boundary values.

Parameters:

Name Type Description Default
b (n,) or (n,p) numpy float array

right-hand side of the linear system

 None
y (o,) or (o,p) numpy float array

what the degrees of freedom are fixed to, u[k] == y or u[k,:] == y

 None

Returns:

Name Type Description
u (n,) or (n,p) numpy float array such that

A[not k,:] * u == b[not k] or A[not k,:] * u == b[not k,:] and u[k] == y or u[k,:] == y.
See Also

min_quad_with_fixed

Examples:

>>> import gpytoolbox as gpy
>>> import numpy as np
>>> import scipy as sp
>>> 
>>> # This matrix is not symmetric!
>>> A = sp.sparse.csc_matrix(np.array([[1.,0.,1.],[0.,1.,0.],[0.,0.,1.]]))
>>> k = np.array([2])
>>> precomp = gpy.fixed_dof_solve_precompute(A, k=k)
>>> 
>>> b = np.array([0.,0.,0.])
>>> y = np.array([5.])
>>> u = precomp.solve(b=b, y=y)
>>> u
array([-5.,  0.,  5.])
>>> 
>>> # A*u matches b, except for indices where constraints apply.
>>> A*u - b
array([0., 0., 5.])
Source code in src/gpytoolbox/fixed_dof_solve.py 
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
def solve(self, b=None, y=None):
    """Solve the prefactored linear system
    while fixing certain degrees of freedom for which the linear system is to be
    ignored during solution.

    For the linear system `A*u = b`, the linear system will be enforced at all
    the rows that don't correspond to fixed degrees of freedom, and the fixed
    degrees of freedom will be enforced on their respective rows.

    This can be used to implement finite differences with Dirichlet boundary
    conditions by fixing the boundary degrees of freedom to the appropriate
    boundary values.

    Parameters
    ----------
    b : (n,) or (n,p) numpy float array
        right-hand side of the linear system
    y : (o,) or (o,p) numpy float array
        what the degrees of freedom are fixed to, `u[k] == y` or `u[k,:] == y`


    Returns
    -------
    u : (n,) or (n,p) numpy float array such that
        `A[not k,:] * u == b[not k]` or `A[not k,:] * u == b[not k,:]`
        and `u[k] == y` or `u[k,:] == y`.


    See Also
    --------
    min_quad_with_fixed


    Examples
    --------
    ```python
    >>> import gpytoolbox as gpy
    >>> import numpy as np
    >>> import scipy as sp
    >>> 
    >>> # This matrix is not symmetric!
    >>> A = sp.sparse.csc_matrix(np.array([[1.,0.,1.],[0.,1.,0.],[0.,0.,1.]]))
    >>> k = np.array([2])
    >>> precomp = gpy.fixed_dof_solve_precompute(A, k=k)
    >>> 
    >>> b = np.array([0.,0.,0.])
    >>> y = np.array([5.])
    >>> u = precomp.solve(b=b, y=y)
    >>> u
    array([-5.,  0.,  5.])
    >>> 
    >>> # A*u matches b, except for indices where constraints apply.
    >>> A*u - b
    array([0., 0., 5.])
    ```

    """

    def cp(x):
        if x is None or np.isscalar(x):
            return 0
        if len(x.shape)==1:
            return 0
        return x.shape[1]
    p = max([cp(b), cp(y)])

    assert b is None or np.isscalar(b) or (p==0 and b.shape==(self.n,)) or (p>0 and b.shape==(self.n,p))
    assert y is None or (np.isscalar(y) and self.o>0) or (p==0 and y.shape==(self.o,)) or (p>0 and y.shape==(self.o,p))

    # Get everything to full dimensions
    if b is None:
        b = 0.
    if np.isscalar(b):
        if p==0:
            b = np.full(self.n, b)
        else:
            b = np.full((self.n,p), b)
    if y is None and self.o>0:
        y = 0.
    if np.isscalar(y) and self.o>0:
        if p==0:
            y = np.full(self.o, y)
        else:
            y = np.full((self.o,p), y)

    # Modified rhs based on known values
    if self.o==0:
        rhs = b
    else:
        b_reduced = b[self.ki] if p==0 else b[self.ki,:]
        rhs = b_reduced - self.A_for_extra_b*y
    ured = self.solver(rhs)

    if self.o==0:
        u = ured
    else:
        if p==0:
            u = np.empty(self.n, dtype=np.float64)
            u[self.ki] = ured
            u[self.k] = y
        else:
            u = np.empty((self.n,p), dtype=np.float64)
            u[self.ki,:] = ured
            u[self.k,:] = y

    return u

fixed_dof_solve(A, b=None, k=None, y=None)

Solves a linear system while fixing certain degrees of freedom for which the linear system is to be ignored during solution.

For the linear system A*u = b, the linear system will be enforced at all the rows that don't correspond to fixed degrees of freedom, and the fixed degrees of freedom will be enforced on their respective rows.

This can be used to implement finite differences with Dirichlet boundary conditions by fixing the boundary degrees of freedom to the appropriate boundary values.

Parameters:

Name Type Description Default
A (n,n) scipy csc matrix

square matrix for the linear system

 required
b (n,) or (n,p) numpy float array

right-hand side of the linear system

 None
k (o,) numpy int array

index vector of fixed degrees of freedom

 None
y (o,) or (o,p) numpy float array

what the degrees of freedom are fixed to, u[k] == y or u[k,:] == y

 None

Returns:

Name Type Description
u (n,) or (n,p) numpy float array such that

A[not k,:] * u == b[not k] or A[not k,:] * u == b[not k,:] and u[k] == y or u[k,:] == y.
See Also

min_quad_with_fixed for solving a quadratic programming problem with linear constraints and a symmetric matrix.

Examples:

>>> import gpytoolbox as gpy
>>> import numpy as np
>>> import scipy as sp
>>> 
>>> # This matrix is not symmetric!
>>> A = sp.sparse.csc_matrix(np.array([[1.,0.,1.],[0.,1.,0.],[0.,0.,1.]]))
>>> b = np.array([0.,0.,0.])
>>> k = np.array([2])
>>> y = np.array([5.])
>>> u = gpy.fixed_dof_solve(A, b, k, y)
>>> u
array([-5.,  0.,  5.])
>>> 
>>> # A*u matches b, except for indices where constraints apply.
>>> A*u - b
array([0., 0., 5.])
Source code in src/gpytoolbox/fixed_dof_solve.py 
 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
60
61
62
63
64
def fixed_dof_solve(A, b=None, k=None, y=None):
    """Solves a linear system while fixing certain degrees of freedom for which
    the linear system is to be ignored during solution.

    For the linear system `A*u = b`, the linear system will be enforced at all
    the rows that don't correspond to fixed degrees of freedom, and the fixed
    degrees of freedom will be enforced on their respective rows.

    This can be used to implement finite differences with Dirichlet boundary
    conditions by fixing the boundary degrees of freedom to the appropriate
    boundary values.

    Parameters
    ----------
    A : (n,n) scipy csc matrix
        square matrix for the linear system
    b : (n,) or (n,p) numpy float array
        right-hand side of the linear system
    k : (o,) numpy int array
        index vector of fixed degrees of freedom
    y : (o,) or (o,p) numpy float array
        what the degrees of freedom are fixed to, `u[k] == y` or `u[k,:] == y`


    Returns
    -------
    u : (n,) or (n,p) numpy float array such that
        `A[not k,:] * u == b[not k]` or `A[not k,:] * u == b[not k,:]`
        and `u[k] == y` or `u[k,:] == y`.


    See Also
    --------
    `min_quad_with_fixed` for solving a quadratic programming problem with
    linear constraints and a symmetric matrix.


    Examples
    --------
    ```python
    >>> import gpytoolbox as gpy
    >>> import numpy as np
    >>> import scipy as sp
    >>> 
    >>> # This matrix is not symmetric!
    >>> A = sp.sparse.csc_matrix(np.array([[1.,0.,1.],[0.,1.,0.],[0.,0.,1.]]))
    >>> b = np.array([0.,0.,0.])
    >>> k = np.array([2])
    >>> y = np.array([5.])
    >>> u = gpy.fixed_dof_solve(A, b, k, y)
    >>> u
    array([-5.,  0.,  5.])
    >>> 
    >>> # A*u matches b, except for indices where constraints apply.
    >>> A*u - b
    array([0., 0., 5.])
    ```
    """

    return fixed_dof_solve_precompute(A, k).solve(b, y)