Skip to content

/ ngsolve

Permalink
Browse files

Merge branch '64-sparsematrix-access-for-fespace-with-dim-1' into 'ma… 

…ster'


Resolve "SparseMatrix access for FESpace with dim != 1"

Closes #64

See merge request !239
  • Loading branch information
@JSchoeberl
JSchoeberl committed Sep 21, 2017
2 parents 12721c7 + 9f133b8 commit fad28e7fa220ab9b2e39f06521269f2a62320027
Unified Split
Showing with 109 additions and 94 deletions.
  1. +17 −0 basiclinalg/python_bla.cpp
  2. +57 −94 linalg/python_linalg.cpp
  3. +35 −0 tests/pytest/test_matrix.py
17 basiclinalg/python_bla.cpp
View file
@@ -477,6 +477,23 @@ void NGS_DLL_HEADER ExportNgbla(py::module & m) {
;
PyDefMatBuffer<Matrix<Complex>>(class_MC);

auto class_Mat2D = py::class_<Mat<2,2,double>>(m,"Mat2D", py::buffer_protocol());
PyDefMatBuffer<Mat<2,2,double>>(class_Mat2D);
class_Mat2D.def("__getitem__", [](Mat<2,2,double> self, py::tuple i)
{ return self(i[0].cast<size_t>(),i[1].cast<size_t>()); });
auto class_Mat2C = py::class_<Mat<2,2,Complex>>(m,"Mat2C", py::buffer_protocol());
PyDefMatBuffer<Mat<2,2,Complex>>(class_Mat2C);
class_Mat2C.def("__getitem__", [](Mat<2,2,Complex> self, py::tuple i)
{ return self(i[0].cast<size_t>(),i[1].cast<size_t>()); });
auto class_Mat3D = py::class_<Mat<3,3,double>>(m,"Mat3D", py::buffer_protocol());
PyDefMatBuffer<Mat<3,3,double>>(class_Mat3D);
class_Mat3D.def("__getitem__", [](Mat<3,3,double> self, py::tuple i)
{ return self(i[0].cast<size_t>(),i[1].cast<size_t>()); });
auto class_Mat3C = py::class_<Mat<3,3,Complex>>(m,"Mat3C", py::buffer_protocol());
PyDefMatBuffer<Mat<3,3,Complex>>(class_Mat3C);
class_Mat3C.def("__getitem__", [](Mat<3,3,Complex> self, py::tuple i)
{ return self(i[0].cast<size_t>(),i[1].cast<size_t>()); });

m.def("Matrix",
[] (int h, int w, bool is_complex) {
if(is_complex) return py::cast(Matrix<Complex>(h,w));
151 linalg/python_linalg.cpp
View file
@@ -4,6 +4,54 @@
using namespace ngla;


template<typename T>
void ExportSparseMatrix(py::module m)
{
py::class_<SparseMatrix<T>, shared_ptr<SparseMatrix<T>>, BaseSparseMatrix, S_BaseMatrix<typename mat_traits<T>::TSCAL>>
(m, (string("SparseMatrix") + typeid(T).name()).c_str(),
"a sparse matrix in CSR storage")
.def("__getitem__",
[](const SparseMatrix<T> & self, py::tuple t)
{
size_t row = t[0].cast<size_t>();
size_t col = t[1].cast<size_t>();
return self(row,col);
})
.def("__setitem__",
[](SparseMatrix<T> & self, py::tuple t, T value)
{
size_t row = t[0].cast<size_t>();
size_t col = t[1].cast<size_t>();
self(row,col) = value;
})

.def("COO", [] (SparseMatrix<T> * sp) -> py::object
{
size_t nze = sp->NZE();
Array<int> ri(nze), ci(nze);
Vector<T> vals(nze);
for (size_t i = 0, ii = 0; i < sp->Height(); i++)
{
FlatArray<int> ind = sp->GetRowIndices(i);
FlatVector<T> rv = sp->GetRowValues(i);
for (int j = 0; j < ind.Size(); j++, ii++)
{
ri[ii] = i;
ci[ii] = ind[j];
vals[ii] = rv[j];
}
}

py::object pyri = py::cast(std::move(ri));
py::object pyci = py::cast(std::move(ci));
py::object pyvals = py::cast(std::move(vals));
return py::make_tuple (pyri, pyci, pyvals);
})
;

py::class_<SparseMatrixSymmetric<T>, shared_ptr<SparseMatrixSymmetric<T>>, SparseMatrix<T>>
(m, (string("SparseMatrixSymmetric") + typeid(T).name()).c_str());
}

void NGS_DLL_HEADER ExportNgla(py::module &m) {

@@ -445,101 +493,16 @@ void NGS_DLL_HEADER ExportNgla(py::module &m) {
;

py::class_<S_BaseMatrix<double>, shared_ptr<S_BaseMatrix<double>>, BaseMatrix>
(m, "S_BaseMatrixD", "base sparse matrix double")
;
(m, "S_BaseMatrixD", "base sparse matrix");
py::class_<S_BaseMatrix<Complex>, shared_ptr<S_BaseMatrix<Complex>>, BaseMatrix>
(m, "S_BaseMatrixC", "base sparse matrix complex")
;

py::class_<SparseMatrix<double>, shared_ptr<SparseMatrix<double>>, BaseSparseMatrix, S_BaseMatrix<double> >
(m, "SparseMatrixD",
"a sparse matrix in CSR storage, entries are real")
.def("__getitem__",
[](const SparseMatrix<double> & self, py::tuple t)
{
size_t row = t[0].cast<size_t>();
size_t col = t[1].cast<size_t>();
return self(row,col);
})
.def("__setitem__",
[](SparseMatrix<double> & self, py::tuple t, double value)
{
size_t row = t[0].cast<size_t>();
size_t col = t[1].cast<size_t>();
self(row,col) = value;
})

.def("COO", [] (SparseMatrix<double> * sp) -> py::object
{
size_t nze = sp->NZE();
Array<int> ri(nze), ci(nze);
Vector<double> vals(nze);
for (size_t i = 0, ii = 0; i < sp->Height(); i++)
{
FlatArray<int> ind = sp->GetRowIndices(i);
FlatVector<double> rv = sp->GetRowValues(i);
for (int j = 0; j < ind.Size(); j++, ii++)
{
ri[ii] = i;
ci[ii] = ind[j];
vals[ii] = rv[j];
}
}

py::object pyri = py::cast(std::move(ri));
py::object pyci = py::cast(std::move(ci));
py::object pyvals = py::cast(std::move(vals));
return py::make_tuple (pyri, pyci, pyvals);
})
;

py::class_<SparseMatrix<Complex>, shared_ptr<SparseMatrix<Complex>>, BaseSparseMatrix, S_BaseMatrix<Complex>>
(m, "SparseMatrixC",
"a sparse matrix in CSR storage, entries are complex")
.def("__getitem__",
[](const SparseMatrix<Complex> & self, py::tuple t)
{
size_t row = t[0].cast<size_t>();
size_t col = t[1].cast<size_t>();
return self(row,col);
})
.def("__setitem__",
[](SparseMatrix<Complex> & self, py::tuple t, Complex value)
{
size_t row = t[0].cast<size_t>();
size_t col = t[1].cast<size_t>();
self(row,col) = value;
})
.def("COO", [] (SparseMatrix<Complex> * spc) -> py::object
{
size_t nze = spc->NZE();
Array<int> ri(nze), ci(nze);
Vector<Complex> vals(nze);
for (size_t i = 0, ii = 0; i < spc->Height(); i++)
{
FlatArray<int> ind = spc->GetRowIndices(i);
FlatVector<Complex> rv = spc->GetRowValues(i);
for (int j = 0; j < ind.Size(); j++, ii++)
{
ri[ii] = i;
ci[ii] = ind[j];
vals[ii] = rv[j];
}
}
py::object pyri = py::cast(std::move(ri));
py::object pyci = py::cast(std::move(ci));
py::object pyvals = py::cast(std::move(vals));
return py::make_tuple (pyri, pyci, pyvals);
})
;


py::class_<SparseMatrixSymmetric<double>, shared_ptr<SparseMatrixSymmetric<double>>, SparseMatrix<double>>
(m, "SparseMatrixSymmetricD");

py::class_<SparseMatrixSymmetric<Complex>, shared_ptr<SparseMatrixSymmetric<Complex>>, SparseMatrix<Complex>>
(m, "SparseMatrixSymmetricC");

(m, "S_BaseMatrixC", "base sparse matrix");

ExportSparseMatrix<double>(m);
ExportSparseMatrix<Complex>(m);
ExportSparseMatrix<Mat<2,2,double>>(m);
ExportSparseMatrix<Mat<2,2,Complex>>(m);
ExportSparseMatrix<Mat<3,3,double>>(m);
ExportSparseMatrix<Mat<3,3,Complex>>(m);

py::class_<BaseBlockJacobiPrecond, shared_ptr<BaseBlockJacobiPrecond>, BaseMatrix>
(m, "BlockSmoother",
35 tests/pytest/test_matrix.py
View file
@@ -1,5 +1,8 @@
import pytest
from ngsolve import *
from netgen.geom2d import unit_square
from netgen.csg import unit_cube
import numpy as np

def test_matrix():
n = 10
@@ -50,3 +53,35 @@ def test_matrix_numpy():
assert d[0,1] == c[0,1]
d[0,1] = 1+3j
assert d[0,1] == c[0,1]

def test_sparsematrix_access():
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
fes = H1(mesh,dim=2)
u,v = fes.TrialFunction(), fes.TestFunction()
a = BilinearForm(fes)
a += SymbolicBFI(InnerProduct(u,v))
a.Assemble()
assert np.linalg.norm(np.array(a.mat[1,1]) - np.array([[0.00550458,0],[0,0.00550458]])) < 1e-8

fes = HCurl(mesh,dim=2,complex=True)
u,v = fes.TrialFunction(), fes.TestFunction()
a = BilinearForm(fes)
a += SymbolicBFI(InnerProduct(u,v))
a.Assemble()
assert abs(a.mat[1,1][0,0] - (0.33333333333325+0j)) < 1e-8

mesh = Mesh(unit_cube.GenerateMesh(maxh=0.2))
fes = H1(mesh,dim=3)
u,v = fes.TrialFunction(), fes.TestFunction()
a = BilinearForm(fes)
a += SymbolicBFI(InnerProduct(u,v))
a.Assemble()
assert np.linalg.norm(np.array(a.mat[1,1]) - np.array([[0.0002156,0.,0.],[0.,0.0002156,0.],[0.,0.,0.0002156]])) < 1e-8

fes = HCurl(mesh,dim=3,complex=True)
u,v = fes.TrialFunction(), fes.TestFunction()
a = BilinearForm(fes)
a += SymbolicBFI(InnerProduct(u,v))
a.Assemble()
assert abs(a.mat[1,1][0,0] - (0.016666666666666604+0j)) < 1e-8

0 comments on commit fad28e7

Please sign in to comment.