Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/tests/utilities/scipy_tensorize.py
ee51b65c7f42master

Last Change
Raw File

scipy_tensorize.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy 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 Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
import scipy.sparse as scsp
from rrompy.solver.scipy import tensorizeLS, detensorizeLS
def test_dense():
N = 5
A1 = np.random.rand(N, N)
x1 = np.random.rand(N)
b1 = A1.dot(x1)
A2 = np.diag(np.arange(1, N + 1) * 2.)
x2 = np.ones(N)
b2 = A2.dot(x2)
A3 = np.eye(N)
x3 = np.ones(N)
b3 = A3.dot(x3)
A, b = tensorizeLS([A1, A2, A3], [b1, b2, b3])
assert np.allclose(A.shape, (3 * N, 3 * N))
assert np.allclose(b.shape, (3 * N,))
x = scsp.linalg.spsolve(A, b)
x1O, x2O, x3O = detensorizeLS(x, 3)
assert np.allclose(x1O, x1, rtol = 1e-8)
assert np.allclose(x2O, x2, rtol = 1e-8)
assert np.allclose(x3O, x3, rtol = 1e-8)
def test_sparse():
N = 5
A1 = scsp.eye(N, format = "csr")
x1 = np.random.rand(N)
b1 = A1.dot(x1)
A2 = scsp.diags(np.arange(1, N + 1) * 2., 0, format = "csr")
x2 = np.ones(N)
b2 = A2.dot(x2)
A, b = tensorizeLS([A1, A2], [b1, b2])
assert np.allclose(A.shape, (2 * N, 2 * N))
assert np.allclose(b.shape, (2 * N,))
x = scsp.linalg.spsolve(A, b)
x1O, x2O = detensorizeLS(x, sizes = [5, 5])
assert np.allclose(x1O, x1, rtol = 1e-8)
assert np.allclose(x2O, x2, rtol = 1e-8)
c4science · Help