← Back to branch summary

~ubuntu-branches/ubuntu/wily/dolfin/wily-proposed

« back to all changes in this revision

Viewing changes to test/unit/python/la/test_matrix.py

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2015-03-17 07:57:11 UTC
  • mfrom: (1.1.18) (19.1.24 experimental)
  • Revision ID: package-import@ubuntu.com-20150317075711-1v207zbty9qmygow
Tags: 1.5.0-1
http://bugs.debian.org/780359
* New upstream release (closes: #780359).
* debian/control:
  - Bump Standards-Version to 3.9.6 (no changes needed).
  - Bump X-Python-Version to >= 2.7.
  - Update package names for new SONAME 1.5 (libdolfin1.4 ->
    libdolfin1.5, libdolfin1.4-dbg -> libdolfin1.5-dbg and
    libdolfin1.4-dev -> libdolfin1.5-dev).
  - Bump minimum required version for python-instant, python-ufl and
    python-ffc to 1.5.0.
  - Add python-sympy and python-six to Depends for binary package
    python-dolfin.
  - Add dh-python to Build-Depends.
  - Remove libcgal-dev from {Build-}Depends.
* Remove CSGCGALMeshGenerator3D-oom.patch since CGAL is no longer used
  by DOLFIN.
* Move debian/libdolfin1.4.install -> debian/libdolfin1.5.install.
* debian/rules: No longer any non DFSG-free stuff to remove, so update
  get-orig-source target (update debian/watch accordingly).
* Update debian/copyright and debian/copyright_hints.

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/unit/python/la/test_matrix.py
 
1
#!/usr/bin/env py.test
 
2
 
 
3
"""Unit tests for the Matrix interface"""
 
4
 
 
5
# Copyright (C) 2011-2014 Garth N. Wells
 
6
#
 
7
# This file is part of DOLFIN.
 
8
#
 
9
# DOLFIN is free software: you can redistribute it and/or modify
 
10
# it under the terms of the GNU Lesser General Public License as published by
 
11
# the Free Software Foundation, either version 3 of the License, or
 
12
# (at your option) any later version.
 
13
#
 
14
# DOLFIN is distributed in the hope that it will be useful,
 
15
# but WITHOUT ANY WARRANTY; without even the implied warranty of
 
16
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 
17
# GNU Lesser General Public License for more details.
 
18
#
 
19
# You should have received a copy of the GNU Lesser General Public License
 
20
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 
21
#
 
22
# Modified by Anders Logg 2011
 
23
# Modified by Mikael Mortensen 2011
 
24
# Modified by Jan Blechta 2013
 
25
 
 
26
from __future__ import print_function
 
27
import pytest
 
28
from dolfin import *
 
29
from six.moves import xrange as range
 
30
 
 
31
from dolfin_utils.test import *
 
32
 
 
33
 
 
34
# TODO: Reuse this fixture setup code between matrix and vector tests:
 
35
 
 
36
# Lists of backends supporting or not supporting data access
 
37
data_backends = []
 
38
no_data_backends = [("PETSc", "")]
 
39
 
 
40
# Add serial only backends
 
41
if MPI.size(mpi_comm_world()) == 1:
 
42
    # TODO: What about "Dense" and "Sparse"? The sub_backend wasn't
 
43
    # used in the old test.
 
44
    data_backends += [("uBLAS", "Dense"), ("uBLAS", "Sparse")]
 
45
    no_data_backends += [("PETScCusp", "")]
 
46
 
 
47
# TODO: STL tests were disabled in old test framework, and do not work now:
 
48
# If we have PETSc, STL Vector gets typedefed to one of these and data
 
49
# test will not work. If none of these backends are available
 
50
# STLVector defaults to uBLASVEctor, which data will work
 
51
#if has_linear_algebra_backend("PETSc"):
 
52
#    no_data_backends += [("STL", "")]
 
53
#else:
 
54
#    data_backends += [("STL", "")]
 
55
 
 
56
# Remove backends we haven't built with
 
57
data_backends = [b for b in data_backends if has_linear_algebra_backend(b[0])]
 
58
no_data_backends = [b for b in no_data_backends if has_linear_algebra_backend(b[0])]
 
59
any_backends = data_backends + no_data_backends
 
60
 
 
61
 
 
62
# Fixtures setting up and resetting the global linear algebra backend for a list of backends
 
63
any_backend     = set_parameters_fixture("linear_algebra_backend", any_backends, lambda x: x[0])
 
64
data_backend    = set_parameters_fixture("linear_algebra_backend", data_backends, lambda x: x[0])
 
65
no_data_backend = set_parameters_fixture("linear_algebra_backend", no_data_backends, lambda x: x[0])
 
66
 
 
67
# With and without explicit backend choice
 
68
use_backend = true_false_fixture
 
69
 
 
70
 
 
71
class TestMatrixForAnyBackend:
 
72
 
 
73
    def assemble_matrices(self, use_backend=False, keep_diagonal=False):
 
74
        " Assemble a pair of matrices, one (square) MxM and one MxN"
 
75
        mesh = UnitSquareMesh(21, 23)
 
76
 
 
77
        V = FunctionSpace(mesh, "Lagrange", 2)
 
78
        W = FunctionSpace(mesh, "Lagrange", 1)
 
79
 
 
80
        v = TestFunction(V)
 
81
        u = TrialFunction(V)
 
82
        s = TrialFunction(W)
 
83
 
 
84
        # Forms
 
85
        a = dot(grad(u), grad(v))*ds
 
86
        b = v*s*dx
 
87
 
 
88
        if use_backend:
 
89
            backend = globals()[self.backend + self.sub_backend + 'Factory'].instance()
 
90
        else:
 
91
            backend = None
 
92
 
 
93
        A = assemble(a, backend=backend, keep_diagonal=keep_diagonal)
 
94
        B = assemble(b, backend=backend, keep_diagonal=keep_diagonal)
 
95
        return A, B
 
96
 
 
97
    def test_basic_la_operations(self, use_backend, any_backend):
 
98
        # Hack to make old tests work in new framework. The original
 
99
        # setup was a bit exoteric...
 
100
        # TODO: Removing use of self in this class will make it
 
101
        # clearer what happens in this test.
 
102
        self.backend, self.sub_backend = any_backend
 
103
 
 
104
        from numpy import ndarray, array, ones, sum
 
105
 
 
106
        # Tests bailout for this choice
 
107
        if self.backend == "uBLAS" and not use_backend:
 
108
            return
 
109
 
 
110
        A, B = self.assemble_matrices(use_backend)
 
111
        unit_norm = A.norm('frobenius')
 
112
 
 
113
        def wrong_getitem(type):
 
114
            if type == 0:
 
115
                A["0,1"]
 
116
            elif type == 1:
 
117
                A[0]
 
118
            elif type == 2:
 
119
                A[0, 0, 0]
 
120
 
 
121
        # Test wrong getitem
 
122
        with pytest.raises(TypeError):
 
123
            wrong_getitem(0)
 
124
        with pytest.raises(TypeError):
 
125
            wrong_getitem(1)
 
126
        with pytest.raises(TypeError):
 
127
            wrong_getitem(2)
 
128
 
 
129
        # Test __imul__ operator
 
130
        A *= 2
 
131
        assert round(A.norm('frobenius') - 2*unit_norm, 7) == 0
 
132
 
 
133
        # Test __idiv__ operator
 
134
        A /= 2
 
135
        assert round(A.norm('frobenius') - unit_norm, 7) == 0
 
136
 
 
137
        # Test __mul__ operator
 
138
        C = 4*A
 
139
        assert round(C.norm('frobenius') - 4*unit_norm, 7) == 0
 
140
 
 
141
        # Test __iadd__ operator
 
142
        A += C
 
143
        assert round(A.norm('frobenius') - 5*unit_norm, 7) == 0
 
144
 
 
145
        # Test __isub__ operator
 
146
        A -= C
 
147
        assert round(A.norm('frobenius') - unit_norm, 7) == 0
 
148
 
 
149
        # Test __mul__ and __add__ operator
 
150
        D = (C+A)*0.2
 
151
        assert round(D.norm('frobenius') - unit_norm, 7) == 0
 
152
 
 
153
        # Test __div__ and __sub__ operator
 
154
        F = (C-A)/3
 
155
        assert round(F.norm('frobenius') - unit_norm, 7) == 0
 
156
 
 
157
        # Test axpy
 
158
        A.axpy(10,C,True)
 
159
        assert round(A.norm('frobenius') - 41*unit_norm, 7) == 0
 
160
 
 
161
        # Test expected size of rectangular array
 
162
        assert A.size(0) == B.size(0)
 
163
        assert B.size(1) == 528
 
164
 
 
165
        # Test setitem/getitem
 
166
        #A[5,5] = 15
 
167
        #assert A[5,5] == 15
 
168
 
 
169
    @skip_in_parallel
 
170
    def test_numpy_array(self, use_backend, any_backend):
 
171
        self.backend, self.sub_backend = any_backend
 
172
 
 
173
        from numpy import ndarray, array, ones, sum
 
174
 
 
175
        # Tests bailout for this choice
 
176
        if self.backend == "uBLAS" and not use_backend:
 
177
            return
 
178
 
 
179
        # Assemble matrices
 
180
        A, B = self.assemble_matrices(use_backend)
 
181
 
 
182
        # Test to NumPy array
 
183
        A2 = A.array()
 
184
        assert isinstance(A2,ndarray)
 
185
        assert A2.shape == (2021, 2021)
 
186
        assert round(sqrt(sum(A2**2)) - A.norm('frobenius'), 7) == 0
 
187
 
 
188
        if self.backend == 'uBLAS' and self.sub_backend == 'Sparse':
 
189
            try:
 
190
                import scipy.sparse
 
191
                import numpy.linalg
 
192
                A3 = A.sparray()
 
193
                assert isinstance(A3, scipy.sparse.csr_matrix)
 
194
                assert round(numpy.linalg.norm(A3.todense() - A2) - 0.0, 7) == 0
 
195
            except ImportError:
 
196
                pass
 
197
 
 
198
    #def create_sparsity_pattern(self):
 
199
    #    "Create a sparsity pattern"
 
200
    #    mesh = UnitSquareMesh(34, 33)
 
201
    #
 
202
    #    V = FunctionSpace(mesh, "Lagrange", 2)
 
203
    #    W = FunctionSpace(mesh, "Lagrange", 1)
 
204
    #
 
205
    #    v = TestFunction(V)
 
206
    #    u = TrialFunction(V)
 
207
    #    s = TrialFunction(W)
 
208
    #
 
209
    #    # Forms
 
210
    #    a = dot(grad(u), grad(v))*dx
 
211
    #    b = v*s*dx
 
212
 
 
213
    def test_create_empty_matrix(self, any_backend):
 
214
        A = Matrix()
 
215
        assert A.size(0) == 0
 
216
        assert A.size(1) == 0
 
217
        info(A)
 
218
        info(A, True)
 
219
 
 
220
    def test_copy_empty_matrix(self, any_backend):
 
221
        A = Matrix()
 
222
        B = Matrix(A)
 
223
        assert B.size(0) == 0
 
224
        assert B.size(1) == 0
 
225
 
 
226
    def test_copy_matrix(self, any_backend):
 
227
        A0, B0 = self.assemble_matrices()
 
228
 
 
229
        A1 = Matrix(A0)
 
230
        assert A0.size(0) == A1.size(0)
 
231
        assert A0.size(1) == A1.size(1)
 
232
        assert A0.norm("frobenius") == A1.norm("frobenius")
 
233
 
 
234
        B1 = Matrix(B0)
 
235
        assert B0.size(0) == B1.size(0)
 
236
        assert B0.size(1) == B1.size(1)
 
237
        assert round(B0.norm("frobenius") - B1.norm("frobenius"), 7) == 0
 
238
 
 
239
    def test_compress_matrix(self, any_backend):
 
240
        self.backend, self.sub_backend = any_backend
 
241
 
 
242
        A, B = self.assemble_matrices()
 
243
        A_norm = A.norm('frobenius')
 
244
        C = Matrix()
 
245
        A.compressed(C)
 
246
        C_norm = C.norm('frobenius')
 
247
        assert round(A_norm - C_norm, 7) == 0
 
248
 
 
249
    @pytest.mark.skipif(MPI.size(mpi_comm_world()) > 1, reason="Disabled because it tends to crash the tests in parallel.")
 
250
    @pytest.mark.slow
 
251
    def test_ident_zeros(self, use_backend, any_backend):
 
252
        self.backend, self.sub_backend = any_backend
 
253
 
 
254
        # Check that PETScMatrix::ident_zeros() rethrows PETSc error
 
255
        if self.backend[0:5] == "PETSc":
 
256
            A, B = self.assemble_matrices(use_backend=use_backend)
 
257
            with pytest.raises(Exception):
 
258
                A.ident_zeros()
 
259
 
 
260
        # Assemble matrix A with diagonal entries
 
261
        A, B = self.assemble_matrices(use_backend=use_backend, keep_diagonal=True)
 
262
 
 
263
        # Find zero rows
 
264
        zero_rows = []
 
265
        for i in range(A.local_range(0)[0], A.local_range(0)[1]):
 
266
            row = A.getrow(i)[1]
 
267
            if sum(abs(row)) < DOLFIN_EPS:
 
268
                zero_rows.append(i)
 
269
 
 
270
        # Set zero rows to (0,...,0, 1, 0,...,0)
 
271
        A.ident_zeros()
 
272
 
 
273
        # Check it
 
274
        for i in zero_rows:
 
275
            cols = A.getrow(i)[0]
 
276
            row  = A.getrow(i)[1]
 
277
            for j in range(cols.size + 1):
 
278
                if i == cols[j]:
 
279
                    assert round(row[j] - 1.0, 7) == 0
 
280
                    break
 
281
            assert j < cols.size
 
282
            assert round(sum(abs(row)) - 1.0, 7) == 0
 
283
 
 
284
    def test_setting_diagonal(self, use_backend, any_backend):
 
285
        self.backend, self.sub_backend = any_backend
 
286
 
 
287
        mesh = UnitSquareMesh(21, 23)
 
288
 
 
289
        V = FunctionSpace(mesh, "Lagrange", 2)
 
290
        W = FunctionSpace(mesh, "Lagrange", 1)
 
291
 
 
292
        v = TestFunction(V)
 
293
        u = TrialFunction(V)
 
294
 
 
295
        if use_backend:
 
296
            backend = globals()[self.backend + self.sub_backend + 'Factory'].instance()
 
297
        else:
 
298
            backend = None
 
299
 
 
300
        B = assemble(u*v*dx(), backend=backend, keep_diagonal=True)
 
301
 
 
302
        b = assemble(action(u*v*dx(), Constant(1)))
 
303
        A = B.copy()
 
304
        A.zero()
 
305
        A.set_diagonal(b)
 
306
 
 
307
        resultsA = Vector()
 
308
        resultsB = Vector()
 
309
        A.init_vector(resultsA, 1)
 
310
        B.init_vector(resultsB, 1)
 
311
 
 
312
        ones = b.copy()
 
313
        ones[:] = 1.0
 
314
 
 
315
        A.mult(ones, resultsA)
 
316
        B.mult(ones, resultsB)
 
317
        assert round(resultsA.norm("l2") - resultsB.norm("l2"), 7) == 0
 
318
 
 
319
    #def test_create_from_sparsity_pattern(self):
 
320
 
 
321
    #def test_size(self):
 
322
 
 
323
    #def test_local_range(self):
 
324
 
 
325
    #def test_zero(self):
 
326
 
 
327
    #def test_apply(self):
 
328
 
 
329
    #def test_str(self):
 
330
 
 
331
    #def test_resize(self):
 
332
 
 
333
 
 
334
    # Test the access of the raw data through pointers
 
335
    # This is only available for uBLAS and MTL4 backends
 
336
    def test_matrix_data(self, use_backend, data_backend):
 
337
        """ Test for ordinary Matrix"""
 
338
        self.backend, self.sub_backend = data_backend
 
339
 
 
340
        # Tests bailout for this choice
 
341
        if self.backend == "uBLAS" and \
 
342
               (not use_backend or self.sub_backend =="Dense"):
 
343
            return
 
344
 
 
345
        A, B = self.assemble_matrices(use_backend)
 
346
        array = A.array()
 
347
        rows, cols, values = A.data()
 
348
        i = 0
 
349
        for row in range(A.size(0)):
 
350
            for col in range(rows[row], rows[row+1]):
 
351
                assert array[row, cols[col]] == values[i]
 
352
                i += 1
 
353
 
 
354
        # Test none writeable of a shallow copy of the data
 
355
        rows, cols, values = A.data(False)
 
356
        def write_data(data):
 
357
            data[0] = 1
 
358
        with pytest.raises(Exception):
 
359
            write_data(rows)
 
360
        with pytest.raises(Exception):
 
361
            write_data(cols)
 
362
        with pytest.raises(Exception):
 
363
            write_data(values)
 
364
 
 
365
        # Test for as_backend_typeed Matrix
 
366
        A = as_backend_type(A)
 
367
        rows, cols, values = A.data()
 
368
        for row in range(A.size(0)):
 
369
            for k in range(rows[row], rows[row+1]):
 
370
                assert array[row,cols[k]] == values[k]
 
371
 
 
372
 
 
373
    def test_matrix_no_data(self, no_data_backend):
 
374
        self.backend, self.sub_backend = no_data_backend
 
375
 
 
376
        A, B = self.assemble_matrices()
 
377
        with pytest.raises(RuntimeError):
 
378
            A.data()
 
379
 
 
380
        A = as_backend_type(A)
 
381
        with pytest.raises(RuntimeError):
 
382
            A.data()
 
383
 
 
384
 
 
385
    def test_matrix_nnz(self, any_backend):
 
386
        A, B = self.assemble_matrices()
 
387
        assert A.nnz() == 2992
 
388
        assert B.nnz() == 9398