Mercurial > fem-fenics-eugenio
changeset 247:8ca45824938e
Add factories hierarchy for matrices' and vectors' assembly
| author | Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> |
|---|---|
| date | Sat, 12 Jul 2014 12:30:13 +0200 |
| parents | a2991e2a085b |
| children | 4522d8fd9cba |
| files | src/Makefile.in src/assemble.cc src/assemble_system.cc src/femfenics_base_factory.h src/femfenics_factory.cc src/femfenics_factory.h src/uBLAS_factory.cc src/uBLAS_factory.h |
| diffstat | 8 files changed, 288 insertions(+), 194 deletions(-) [+] |
line wrap: on
line diff
--- a/src/Makefile.in Sat Jun 28 16:22:13 2014 +0200 +++ b/src/Makefile.in Sat Jul 12 12:30:13 2014 +0200 @@ -17,6 +17,7 @@ MKOCTFILE ?= mkoctfile FFC ?= ffc +ARFLAGS=-r -s CPPFLAGS='@DENSE_CPPFLAGS@ @DOLFIN_CPPFLAGS@ @EIGEN_CPPFLAGS@ -I.' LIBS_RAW=@EIGEN_LIBS@ @DOLFIN_LIBS@ LIBS=$(patsubst %, "%", $(LIBS_RAW)) @@ -73,16 +74,18 @@ Function.o: Function.cc function.h CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c Function.cc -o $@ -assemble.oct: assemble.o - CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) assemble.o -o $@ $(LIBS) +assemble.oct: assemble.o libfemfenics_factories.a + CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) assemble.o -o $@ \ + libfemfenics_factories.a $(LIBS) -assemble.o: assemble.cc form.h boundarycondition.h +assemble.o: assemble.cc form.h boundarycondition.h femfenics_factory.h CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c assemble.cc -o $@ -assemble_system.oct: assemble_system.o - CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) assemble_system.o -o $@ $(LIBS) +assemble_system.oct: assemble_system.o libfemfenics_factories.a + CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) assemble_system.o -o $@ \ + libfemfenics_factories.a $(LIBS) -assemble_system.o: assemble_system.cc form.h boundarycondition.h +assemble_system.o: assemble_system.cc form.h boundarycondition.h femfenics_factory.h CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c assemble_system.cc -o $@ save.oct: save.o mkfunction @@ -130,8 +133,18 @@ interpolate.o: interpolate.cc function.h CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c interpolate.cc -o $@ +femfenics_factory.o: femfenics_factory.cc femfenics_factory.h femfenics_base_factory.h\ + uBLAS_factory.h + CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c $< -o $@ + +uBLAS_factory.o: uBLAS_factory.cc uBLAS_factory.h femfenics_base_factory.h + CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c $< -o $@ + +libfemfenics_factories.a: femfenics_factory.o uBLAS_factory.o + $(AR) $(ARFLAGS) $@ $^ + clean: - $(RM) *.o core octave-core *.oct *~ *.xml + $(RM) *.o core octave-core *.oct *~ *.xml *.a cleanall: $(RM) *.o core octave-core *.oct *~ *.xml *.status *.log \
--- a/src/assemble.cc Sat Jun 28 16:22:13 2014 +0200 +++ b/src/assemble.cc Sat Jul 12 12:30:13 2014 +0200 @@ -17,6 +17,7 @@ #include "form.h" #include "boundarycondition.h" +#include "femfenics_factory.h" DEFUN_DLD (assemble, args, nargout, "-*- texinfo -*-\n\ @@ -79,7 +80,8 @@ if (a.rank () == 2) { - dolfin::parameters["linear_algebra_backend"] = "uBLAS"; + femfenics_factory factory; + dolfin::Matrix A; dolfin::assemble (A, a); @@ -104,55 +106,13 @@ error ("assemble: unknown argument type"); } - // Get capacity of the dolfin sparse matrix - boost::tuples::tuple<const std::size_t*, - const std::size_t*, - const double*, int> - aa = A.data (); - - int nnz = aa.get<3> (); - std::size_t nr = A.size (0), nc = A.size (1); - std::vector<double> data_tmp; - std::vector<std::size_t> cidx_tmp; - - dim_vector dims (nnz, 1); - octave_idx_type nz = 0, ii = 0; - Array<octave_idx_type> - ridx (dims, 0), - cidx (dims, 0); - Array<double> data (dims, 0); - - octave_idx_type* orow = ridx.fortran_vec (); - octave_idx_type* oc = cidx.fortran_vec (); - double* ov = data.fortran_vec (); - - for (std::size_t i = 0; i < nr; ++i) - { - A.getrow (i, cidx_tmp, data_tmp); - nz += cidx_tmp.size (); - - for (octave_idx_type j = 0; - j < cidx_tmp.size (); ++j) - { - orow [ii + j] = i; - oc [ii + j] = cidx_tmp [j]; - ov [ii + j] = data_tmp [j]; - } - - ii = nz; - } - - dims(0) = ii; - ridx.resize (dims); - cidx.resize (dims); - data.resize (dims); - - SparseMatrix sm (data, ridx, cidx, nr, nc); - retval(0) = sm; + retval(0) = factory.matrix (A); } else if (a.rank () == 1) { + femfenics_factory factory; + dolfin::Vector A; dolfin::assemble (A, a); @@ -176,16 +136,7 @@ error ("assemble: unknown argument type"); } - dim_vector dims; - dims.resize (2); - dims(0) = A.size (); - dims(1) = 1; - Array<double> myb (dims); - - for (std::size_t i = 0; i < A.size (); ++i) - myb.xelem (i) = A[i]; - - retval(0) = myb; + retval(0) = factory.vector (A); } else if (a.rank () == 0) @@ -219,6 +170,8 @@ if (a.rank () == 1) { + femfenics_factory factory; + dolfin::Vector A; dolfin::assemble (A, a); @@ -247,20 +200,8 @@ error ("assemble: unknown argument type"); } - dim_vector dims; - dims.resize (2); - dims(0) = A.size (); - dims(1) = 1; - Array<double> myb (dims), myc (dims); - - for (std::size_t i = 0; i < A.size (); ++i) - { - myb.xelem (i) = A[i]; - myc.xelem (i) = x[i]; - } - - retval(0) = myb; - retval(1) = myc; + retval(0) = factory.vector (A); + retval(1) = factory.vector (x); } else
--- a/src/assemble_system.cc Sat Jun 28 16:22:13 2014 +0200 +++ b/src/assemble_system.cc Sat Jul 12 12:30:13 2014 +0200 @@ -17,6 +17,7 @@ #include "form.h" #include "boundarycondition.h" +#include "femfenics_factory.h" DEFUN_DLD (assemble_system, args, nargout, "-*- texinfo -*-\n\ @@ -79,7 +80,8 @@ if (a.rank () == 2 && b.rank () == 1) { - dolfin::parameters["linear_algebra_backend"] = "uBLAS"; + femfenics_factory factory; + dolfin::Matrix A; dolfin::assemble (A, a); dolfin::Vector B; @@ -105,62 +107,8 @@ error ("assemble_system: unknown argument type"); } - // Get capacity of the dolfin sparse matrix - boost::tuples::tuple<const std::size_t*, - const std::size_t*, - const double*, int> - aa = A.data (); - - int nnz = aa.get<3> (); - std::size_t nr = A.size (0), nc = A.size (1); - std::vector<double> data_tmp; - std::vector<std::size_t> cidx_tmp; - - dim_vector dims (nnz, 1); - octave_idx_type nz = 0, ii = 0; - Array<octave_idx_type> - ridx (dims, 0), - cidx (dims, 0); - Array<double> data (dims, 0); - - octave_idx_type* orow = ridx.fortran_vec (); - octave_idx_type* oc = cidx.fortran_vec (); - double* ov = data.fortran_vec (); - - for (std::size_t i = 0; i < nr; ++i) - { - A.getrow (i, cidx_tmp, data_tmp); - nz += cidx_tmp.size (); - - for (octave_idx_type j = 0; - j < cidx_tmp.size (); ++j) - { - orow [ii + j] = i; - oc [ii + j] = cidx_tmp [j]; - ov [ii + j] = data_tmp [j]; - } - - ii = nz; - } - - dims(0) = ii; - ridx.resize (dims); - cidx.resize (dims); - data.resize (dims); - - SparseMatrix sm (data, ridx, cidx, nr, nc); - retval(0) = sm; - - dim_vector dim; - dim.resize (2); - dim(0) = B.size (); - dim(1) = 1; - Array<double> myb (dim); - - for (std::size_t i = 0; i < B.size (); ++i) - myb.xelem (i) = B[i]; - - retval(1) = myb; + retval(0) = factory.matrix (A); + retval(1) = factory.vector (B); } } else @@ -192,7 +140,8 @@ if (a.rank () == 2 && b.rank () == 1) { - dolfin::parameters["linear_algebra_backend"] = "uBLAS"; + femfenics_factory factory; + dolfin::Matrix A; dolfin::assemble (A, a); dolfin::Vector B; @@ -224,65 +173,9 @@ error ("assemble_system: unknown argument type"); } - // Get capacity of the dolfin sparse matrix - boost::tuples::tuple<const std::size_t*, - const std::size_t*, - const double*, int> - aa = A.data (); - int nnz = aa.get<3> (); - std::size_t nr = A.size (0), nc = A.size (1); - std::vector<double> data_tmp; - std::vector<std::size_t> cidx_tmp; - - dim_vector dims (nnz, 1); - octave_idx_type nz = 0, ii = 0; - Array<octave_idx_type> - ridx (dims, 0), - cidx (dims, 0); - Array<double> data (dims, 0); - - octave_idx_type* orow = ridx.fortran_vec (); - octave_idx_type* oc = cidx.fortran_vec (); - double* ov = data.fortran_vec (); - - for (std::size_t i = 0; i < nr; ++i) - { - A.getrow (i, cidx_tmp, data_tmp); - nz += cidx_tmp.size (); - - for (octave_idx_type j = 0; - j < cidx_tmp.size (); ++j) - { - orow [ii + j] = i; - oc [ii + j] = cidx_tmp [j]; - ov [ii + j] = data_tmp [j]; - } - - ii = nz; - } - - dims(0) = ii; - ridx.resize (dims); - cidx.resize (dims); - data.resize (dims); - - SparseMatrix sm (data, ridx, cidx, nr, nc); - retval(0) = sm; - - dim_vector dim; - dim.resize (2); - dim(0) = B.size (); - dim(1) = 1; - Array<double> myb (dim), myc (dim); - - for (std::size_t i = 0; i < B.size (); ++i) - { - myb.xelem (i) = B[i]; - myc.xelem (i) = x[i]; - } - - retval(1) = myb; - retval(2) = myc; + retval(0) = factory.matrix (A); + retval(1) = factory.vector (B); + retval(2) = factory.vector (x); } } else
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/femfenics_base_factory.h Sat Jul 12 12:30:13 2014 +0200 @@ -0,0 +1,35 @@ +/* + Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> + + 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/>. +*/ + +#ifndef __FEMFENICS_BASE_FACTORY__ +#define __FEMFENICS_BASE_FACTORY__ + +#include <dolfin.h> +#include <octave/oct.h> + +class femfenics_base_factory +{ + public: + + femfenics_base_factory () {} + virtual ~femfenics_base_factory () {} + + virtual octave_value matrix (dolfin::Matrix const&) const =0; + virtual octave_value vector (dolfin::Vector const&) const =0; +}; + +#endif \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/femfenics_factory.cc Sat Jul 12 12:30:13 2014 +0200 @@ -0,0 +1,27 @@ +/* + Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> + + 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/>. +*/ + +#include "femfenics_factory.h" +#include "uBLAS_factory.h" + +femfenics_base_factory const& +femfenics_factory::factory (void) const + { + //FIXME: Since just one backend has its interface implemented, this method + // doesn't check which one to provide, yet + return uBLAS_factory::instance (); + } \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/femfenics_factory.h Sat Jul 12 12:30:13 2014 +0200 @@ -0,0 +1,46 @@ +/* + Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> + + 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/>. +*/ + +#ifndef __FEMFENICS_FACTORY__ +#define __FEMFENICS_FACTORY__ + +#include "femfenics_base_factory.h" + +class femfenics_factory : public femfenics_base_factory +{ + public: + + femfenics_factory () + { dolfin::parameters["linear_algebra_backend"] = linear_algebra_backend (); } + virtual ~femfenics_factory () {} + + //FIXME: just one backend implemented + inline char const * linear_algebra_backend (void) const + { return "uBLAS"; } + + virtual inline octave_value matrix (dolfin::Matrix const& A) const + { return factory ().matrix (A); } + + virtual inline octave_value vector (dolfin::Vector const& b) const + { return factory ().vector (b); } + + private: + + femfenics_base_factory const& factory (void) const; +}; + +#endif \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/uBLAS_factory.cc Sat Jul 12 12:30:13 2014 +0200 @@ -0,0 +1,98 @@ +/* + Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> + + 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/>. +*/ + +#include "uBLAS_factory.h" + +femfenics_base_factory const& +uBLAS_factory::instance (void) + { + static uBLAS_factory const theinstance; + return theinstance; + } + +octave_value +uBLAS_factory::matrix (dolfin::Matrix const& A) const + { + octave_value retval; + + // Get capacity of the dolfin sparse matrix + boost::tuples::tuple<const std::size_t*, + const std::size_t*, + const double*, int> + aa = A.data (); + + int nnz = aa.get<3> (); + std::size_t nr = A.size (0), nc = A.size (1); + std::vector<double> data_tmp; + std::vector<std::size_t> cidx_tmp; + + dim_vector dims (nnz, 1); + octave_idx_type nz = 0, ii = 0; + Array<octave_idx_type> + ridx (dims, 0), + cidx (dims, 0); + Array<double> data (dims, 0); + + octave_idx_type* orow = ridx.fortran_vec (); + octave_idx_type* oc = cidx.fortran_vec (); + double* ov = data.fortran_vec (); + + for (std::size_t i = 0; i < nr; ++i) + { + A.getrow (i, cidx_tmp, data_tmp); + nz += cidx_tmp.size (); + + for (octave_idx_type j = 0; + j < cidx_tmp.size (); ++j) + { + orow [ii + j] = i; + oc [ii + j] = cidx_tmp [j]; + ov [ii + j] = data_tmp [j]; + } + + ii = nz; + } + + dims(0) = ii; + ridx.resize (dims); + cidx.resize (dims); + data.resize (dims); + + SparseMatrix sm (data, ridx, cidx, nr, nc); + retval = sm; + + return retval; + } + +octave_value +uBLAS_factory::vector (dolfin::Vector const& b) const + { + octave_value retval; + + dim_vector dims; + dims.resize (2); + dims(0) = b.size (); + dims(1) = 1; + Array<double> myb (dims); + + for (std::size_t i = 0; i < b.size (); ++i) + myb.xelem (i) = b[i]; + + retval = myb; + + return retval; + } \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/uBLAS_factory.h Sat Jul 12 12:30:13 2014 +0200 @@ -0,0 +1,41 @@ +/* + Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> + + 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/>. +*/ + +#ifndef __FEMFENICS_UBLAS_FACTORY__ +#define __FEMFENICS_UBLAS_FACTORY__ + +#include "femfenics_base_factory.h" + +class uBLAS_factory : public femfenics_base_factory +{ + public: + + virtual ~uBLAS_factory () {} + + octave_value matrix (dolfin::Matrix const&) const; + octave_value vector (dolfin::Vector const&) const; + + static femfenics_base_factory const& instance (void); + + private: + + uBLAS_factory () {} + uBLAS_factory (uBLAS_factory const&); + uBLAS_factory & operator = (uBLAS_factory const&); +}; + +#endif \ No newline at end of file