diff --git a/tm/src/CSRMatrix.hpp b/tm/src/CSRMatrix.hpp new file mode 100644 index 0000000..abebacb --- /dev/null +++ b/tm/src/CSRMatrix.hpp @@ -0,0 +1,146 @@ +#ifndef _CSRMatrix_hpp_ +#define _CSRMatrix_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include +#include +#ifdef HAVE_MPI +#include +#endif + +namespace miniFE { + +template +struct +CSRMatrix { + CSRMatrix() + : has_local_indices(false), + rows(), row_offsets(), row_offsets_external(), + packed_cols(), packed_coefs(), + num_cols(0) +#ifdef HAVE_MPI + ,external_index(), external_local_index(), elements_to_send(), + neighbors(), recv_length(), send_length(), send_buffer(), request() +#endif + { + } + + ~CSRMatrix() + {} + + typedef Scalar ScalarType; + typedef LocalOrdinal LocalOrdinalType; + typedef GlobalOrdinal GlobalOrdinalType; + + bool has_local_indices; + std::vector rows; + std::vector row_offsets; + std::vector row_offsets_external; + std::vector packed_cols; + std::vector packed_coefs; + LocalOrdinal num_cols; + +#ifdef HAVE_MPI + std::vector external_index; + std::vector external_local_index; + std::vector elements_to_send; + std::vector neighbors; + std::vector recv_length; + std::vector send_length; + std::vector send_buffer; + std::vector request; +#endif + + size_t num_nonzeros() const + { + return row_offsets[row_offsets.size()-1]; + } + + void reserve_space(unsigned nrows, unsigned ncols_per_row) + { + rows.resize(nrows); + row_offsets.resize(nrows+1); + packed_cols.reserve(nrows * ncols_per_row); + packed_coefs.reserve(nrows * ncols_per_row); + + #pragma omp parallel for + for(MINIFE_GLOBAL_ORDINAL i = 0; i < nrows; ++i) { + rows[i] = 0; + row_offsets[i] = 0; + } + + #pragma omp parallel for + for(MINIFE_GLOBAL_ORDINAL i = 0; i < (nrows * ncols_per_row); ++i) { + packed_cols[i] = 0; + packed_coefs[i] = 0; + } + } + + void get_row_pointers(GlobalOrdinalType row, size_t& row_length, + GlobalOrdinalType*& cols, + ScalarType*& coefs) + { + ptrdiff_t local_row = -1; + //first see if we can get the local-row index using fast direct lookup: + if (rows.size() >= 1) { + ptrdiff_t idx = row - rows[0]; + if (idx < rows.size() && rows[idx] == row) { + local_row = idx; + } + } + + //if we didn't get the local-row index using direct lookup, try a + //more expensive binary-search: + if (local_row == -1) { + typename std::vector::iterator row_iter = + std::lower_bound(rows.begin(), rows.end(), row); + + //if we still haven't found row, it's not local so jump out: + if (row_iter == rows.end() || *row_iter != row) { + row_length = 0; + return; + } + + local_row = row_iter - rows.begin(); + } + + LocalOrdinalType offset = row_offsets[local_row]; + row_length = row_offsets[local_row+1] - offset; + cols = &packed_cols[offset]; + coefs = &packed_coefs[offset]; + } +}; + +}//namespace miniFE + +#endif + diff --git a/tm/src/ELLMatrix.hpp b/tm/src/ELLMatrix.hpp new file mode 100644 index 0000000..f405f6d --- /dev/null +++ b/tm/src/ELLMatrix.hpp @@ -0,0 +1,139 @@ +#ifndef _ELLMatrix_hpp_ +#define _ELLMatrix_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include +#include +#ifdef HAVE_MPI +#include +#endif + +namespace miniFE { + +template +struct +ELLMatrix { + ELLMatrix() + : has_local_indices(false), + rows(), + cols(), coefs(), + num_cols(0), + num_cols_per_row(0) +#ifdef HAVE_MPI + ,external_index(), external_local_index(), elements_to_send(), + neighbors(), recv_length(), send_length(), send_buffer(), request() +#endif + { + } + + ~ELLMatrix() + {} + + typedef Scalar ScalarType; + typedef LocalOrdinal LocalOrdinalType; + typedef GlobalOrdinal GlobalOrdinalType; + + bool has_local_indices; + std::vector rows; + std::vector cols; + std::vector coefs; + LocalOrdinal num_cols; + LocalOrdinal num_cols_per_row; + +#ifdef HAVE_MPI + std::vector external_index; + std::vector external_local_index; + std::vector elements_to_send; + std::vector neighbors; + std::vector recv_length; + std::vector send_length; + std::vector send_buffer; + std::vector request; +#endif + + size_t num_nonzeros() const + { + return rows.size()*num_cols_per_row; + } + + void reserve_space(unsigned nrows, unsigned ncols_per_row) + { + rows.resize(nrows); + cols.resize(nrows * ncols_per_row); + coefs.resize(nrows * ncols_per_row); + num_cols_per_row = ncols_per_row; + } + + void get_row_pointers(GlobalOrdinalType row, size_t& row_length, + GlobalOrdinalType*& cols_ptr, + ScalarType*& coefs_ptr) + { + ptrdiff_t local_row = -1; + //first see if we can get the local-row index using fast direct lookup: + if (rows.size() >= 1) { + ptrdiff_t idx = row - rows[0]; + if (idx < rows.size() && rows[idx] == row) { + local_row = idx; + } + } + + //if we didn't get the local-row index using direct lookup, try a + //more expensive binary-search: + if (local_row == -1) { + typename std::vector::iterator row_iter = + std::lower_bound(rows.begin(), rows.end(), row); + + //if we still haven't found row, it's not local so jump out: + if (row_iter == rows.end() || *row_iter != row) { + row_length = 0; + return; + } + + local_row = row_iter - rows.begin(); + } + + cols_ptr = &cols[local_row*num_cols_per_row]; + coefs_ptr = &coefs[local_row*num_cols_per_row]; + + int idx = num_cols_per_row-1; + while(idx>=0) { + if (cols_ptr[idx] != 0) break; + --idx; + } + row_length = idx+1; + } +}; + +}//namespace miniFE + +#endif + diff --git a/tm/src/GetNodesCoords.hpp b/tm/src/GetNodesCoords.hpp new file mode 100644 index 0000000..5278dd1 --- /dev/null +++ b/tm/src/GetNodesCoords.hpp @@ -0,0 +1,51 @@ +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#ifndef _GETNODESCOORDS_HPP_ +#define _GETNODESCOORDS_HPP_ + +#include +#include + +template +struct GetNodesCoords { + const miniFE::simple_mesh_description* mesh; + GlobalOrdinal* elemIDs; + GlobalOrdinal* node_ordinals; + Scalar* elem_node_coords; + +inline void operator()(int i) +{ + unsigned nnodes = miniFE::Hex8::numNodesPerElem; + GlobalOrdinal elemID = elemIDs[i]; + GlobalOrdinal* node_ords = node_ordinals+i*nnodes; + Scalar* node_coords = elem_node_coords+i*nnodes*miniFE::Hex8::spatialDim; + get_elem_nodes_and_coords(*mesh, elemID, node_ords, node_coords); +} +}; + +#endif diff --git a/tm/src/Hex8_box_utils.hpp b/tm/src/Hex8_box_utils.hpp new file mode 100644 index 0000000..da38684 --- /dev/null +++ b/tm/src/Hex8_box_utils.hpp @@ -0,0 +1,173 @@ +#ifndef _Hex8_box_utils_hpp_ +#define _Hex8_box_utils_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include + +#include +#include +#include +#include + +namespace miniFE { + + +template +void get_hex8_node_ids(int nx, int ny, + GlobalOrdinal node0, + GlobalOrdinal* elem_node_ids) +{ +//Given box dimensions nx and ny, and a starting node +//(local-node-0 for a hex8), compute the other nodes +//of the hex8 using the exodus ordering convention. + elem_node_ids[0] = node0; + elem_node_ids[1] = node0 + 1; + elem_node_ids[2] = node0 + nx + 1; + elem_node_ids[3] = node0 + nx; + elem_node_ids[4] = node0 + nx*ny; + elem_node_ids[5] = node0 + 1 + nx*ny; + elem_node_ids[6] = node0 + nx + nx*ny + 1; + elem_node_ids[7] = node0 + nx + nx*ny; +} + +template +void get_hex8_node_coords_3d(Scalar x, Scalar y, Scalar z, + Scalar hx, Scalar hy, Scalar hz, + Scalar* elem_node_coords) +{ + //Input: x,y,z are the coordinates of local-node 0 for a Hex8. + //'hx', 'hy', 'hz' are the lengths of the sides of the element + //in each direction. + + elem_node_coords[0] = x; + elem_node_coords[1] = y; + elem_node_coords[2] = z; + + elem_node_coords[3] = x + hx; + elem_node_coords[4] = y; + elem_node_coords[5] = z; + + elem_node_coords[6] = x + hx; + elem_node_coords[7] = y + hy; + elem_node_coords[8] = z; + + elem_node_coords[9] = x; + elem_node_coords[10] = y + hy; + elem_node_coords[11] = z; + + elem_node_coords[12] = x; + elem_node_coords[13] = y; + elem_node_coords[14] = z + hz; + + elem_node_coords[15] = x + hx; + elem_node_coords[16] = y; + elem_node_coords[17] = z + hz; + + elem_node_coords[18] = x + hx; + elem_node_coords[19] = y + hy; + elem_node_coords[20] = z + hz; + + elem_node_coords[21] = x; + elem_node_coords[22] = y + hy; + elem_node_coords[23] = z + hz; +} + +template +void +get_elem_nodes_and_coords(const simple_mesh_description& mesh, + GlobalOrdinal elemID, + GlobalOrdinal* node_ords, Scalar* node_coords) +{ + int global_nodes_x = mesh.global_box[0][1]+1; + int global_nodes_y = mesh.global_box[1][1]+1; + int global_nodes_z = mesh.global_box[2][1]+1; + + if (elemID < 0) { + //I don't think this can happen, but check for the sake of paranoia... + throw std::runtime_error("get_elem_nodes_and_coords ERROR, negative elemID"); + } + + int elem_int_x, elem_int_y, elem_int_z; + get_int_coords(elemID, global_nodes_x-1, global_nodes_y-1, global_nodes_z-1, + elem_int_x, elem_int_y, elem_int_z); + GlobalOrdinal nodeID = get_id(global_nodes_x, global_nodes_y, global_nodes_z, elem_int_x, elem_int_y, elem_int_z); + +#ifdef MINIFE_DEBUG_VERBOSE + std::cout<<"\nelemID: "<(nodeID, global_nodes_x,global_nodes_y,global_nodes_z, + ix,iy,iz); + Scalar hx = 1.0/global_elems_x; + Scalar hy = 1.0/global_elems_y; + Scalar hz = 1.0/global_elems_z; + get_hex8_node_coords_3d(ix, iy, iz, hx, hy, hz, node_coords); +#ifdef MINIFE_DEBUG_VERBOSE + int offset = 0; + for(int i=0; i +void +get_elem_nodes_and_coords(const simple_mesh_description& mesh, + GlobalOrdinal elemID, + ElemData& elem_data) +{ + get_elem_nodes_and_coords(mesh, elemID, elem_data.elem_node_ids, elem_data.elem_node_coords); +} + +}//namespace miniFE + +#endif + diff --git a/tm/src/Makefile b/tm/src/Makefile new file mode 100644 index 0000000..e15b028 --- /dev/null +++ b/tm/src/Makefile @@ -0,0 +1,44 @@ + +#SHELL = /home/projects/x86-64/time/1.9.0/bin/time sh + +#----------------------------------------------------------------------- +# This file compiles for OpenMP and MPI hybrid operations using the GNU +# compile chain. + +MINIFE_TYPES = \ + -DMINIFE_SCALAR=double \ + -DMINIFE_LOCAL_ORDINAL=int \ + -DMINIFE_GLOBAL_ORDINAL=int + +MINIFE_MATRIX_TYPE = -DMINIFE_CSR_MATRIX +#MINIFE_MATRIX_TYPE = -DMINIFE_ELL_MATRIX + +#----------------------------------------------------------------------- + +CFLAGS = -O3 -fopenmp -xcore-avx512 +CXXFLAGS = $(CFLAGS) + +CPPFLAGS = -I. -I../utils -I../fem $(MINIFE_TYPES) \ + $(MINIFE_MATRIX_TYPE) #\ +# -DHAVE_MPI -DMPICH_IGNORE_CXX_SEEK \ +# -DMINIFE_REPORT_RUSAGE + +LDFLAGS=$(CFLAGS) +LIBS= + +# The MPICH_IGNORE_CXX_SEEK macro is required for some mpich versions, +# such as the one on my cygwin machine. + +#CXX=mpiicpc +#CC=mpiicc + +#CXX=g++ +#CC=g++ + +CXX=icpc +CC=icc + +#CXX=mpicxx +#CC=mpicc + +include make_targets diff --git a/tm/src/Makefile.cray.xc30 b/tm/src/Makefile.cray.xc30 new file mode 100644 index 0000000..2710a45 --- /dev/null +++ b/tm/src/Makefile.cray.xc30 @@ -0,0 +1,35 @@ +#----------------------------------------------------------------------- +# ATTENTION: +# +# This file does not enable MPI in the binary, this is just OpenMP +# only. To enable MPI calls add -DHAVE_MPI to CPPFLAGS +# +#----------------------------------------------------------------------- + +MINIFE_TYPES = \ + -DMINIFE_SCALAR=double \ + -DMINIFE_LOCAL_ORDINAL=int \ + -DMINIFE_GLOBAL_ORDINAL=int + +MINIFE_MATRIX_TYPE = -DMINIFE_CSR_MATRIX +# MINIFE_MATRIX_TYPE = -DMINIFE_ELL_MATRIX + +#----------------------------------------------------------------------- + +MPIDIR=/opt/cray/mpt/default/gni/mpich2-cray64/81 + +CFLAGS = -O3 +CXXFLAGS = $(CFLAGS) -I$(MPIDIR)/include + +CPPFLAGS = -I. -I../utils -I../fem $(MINIFE_TYPES) $(MINIFE_MATRIX_TYPE) \ + -DHAVE_MPI + +LDFLAGS=-L$(MPIDIR)/lib +LIBS=-lmpichcxx + +CXX=CC +CC=cc + +#----------------------------------------------------------------------- + +include make_targets diff --git a/tm/src/Makefile.gnu.openmp b/tm/src/Makefile.gnu.openmp new file mode 100644 index 0000000..8fabd36 --- /dev/null +++ b/tm/src/Makefile.gnu.openmp @@ -0,0 +1,32 @@ +#----------------------------------------------------------------------- +# ATTENTION: +# +# This file does not enable MPI in the binary, this is just OpenMP +# only. To enable MPI calls add -DHAVE_MPI to CPPFLAGS +# +#----------------------------------------------------------------------- + +MINIFE_TYPES = \ + -DMINIFE_SCALAR=double \ + -DMINIFE_LOCAL_ORDINAL=int \ + -DMINIFE_GLOBAL_ORDINAL=int + +MINIFE_MATRIX_TYPE = -DMINIFE_CSR_MATRIX +# MINIFE_MATRIX_TYPE = -DMINIFE_ELL_MATRIX + +#----------------------------------------------------------------------- + +CFLAGS = -O3 -fopenmp +CXXFLAGS = $(CFLAGS) + +CPPFLAGS = -I. -I../utils -I../fem $(MINIFE_TYPES) $(MINIFE_MATRIX_TYPE) + +LDFLAGS= +LIBS= + +CXX=g++ +CC=gcc + +#----------------------------------------------------------------------- + +include make_targets diff --git a/tm/src/Makefile.intel.openmp b/tm/src/Makefile.intel.openmp new file mode 100644 index 0000000..390f163 --- /dev/null +++ b/tm/src/Makefile.intel.openmp @@ -0,0 +1,32 @@ +#----------------------------------------------------------------------- +# ATTENTION: +# +# This file does not enable MPI in the binary, this is just OpenMP +# only. To enable MPI calls add -DHAVE_MPI to CPPFLAGS +# +#----------------------------------------------------------------------- + +MINIFE_TYPES = \ + -DMINIFE_SCALAR=double \ + -DMINIFE_LOCAL_ORDINAL=int \ + -DMINIFE_GLOBAL_ORDINAL=int + +MINIFE_MATRIX_TYPE = -DMINIFE_CSR_MATRIX +# MINIFE_MATRIX_TYPE = -DMINIFE_ELL_MATRIX + +#----------------------------------------------------------------------- + +CFLAGS = -O3 -fopenmp +CXXFLAGS = $(CFLAGS) + +CPPFLAGS = -I. -I../utils -I../fem $(MINIFE_TYPES) $(MINIFE_MATRIX_TYPE) + +LDFLAGS= +LIBS= + +CXX=icpc +CC=icc + +#----------------------------------------------------------------------- + +include make_targets diff --git a/tm/src/MatrixCopyOp.hpp b/tm/src/MatrixCopyOp.hpp new file mode 100644 index 0000000..9d24be0 --- /dev/null +++ b/tm/src/MatrixCopyOp.hpp @@ -0,0 +1,60 @@ +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#ifndef _MatrixCopyOp_hpp_ +#define _MatrixCopyOp_hpp_ + +template +struct MatrixCopyOp { + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinalType; + typedef typename MatrixType::LocalOrdinalType LocalOrdinalType; + typedef typename MatrixType::ScalarType ScalarType; + + const GlobalOrdinalType* src_rows; + const LocalOrdinalType* src_rowoffsets; + const GlobalOrdinalType* src_cols; + const ScalarType* src_coefs; + + GlobalOrdinalType* dest_rows; + LocalOrdinalType* dest_rowoffsets; + GlobalOrdinalType* dest_cols; + ScalarType* dest_coefs; + int n; + + inline void operator()(int i) + { + dest_rows[i] = src_rows[i]; + dest_rowoffsets[i] = src_rowoffsets[i]; + for(int j=src_rowoffsets[i]; j +#include + +#include +#include + +#include + +template +void sort_if_needed(GlobalOrdinal* list, + GlobalOrdinal list_len) +{ + bool need_to_sort = false; + for(GlobalOrdinal i=list_len-1; i>=1; --i) { + if (list[i] < list[i-1]) { + need_to_sort = true; + break; + } + } + + if (need_to_sort) { + std::sort(list,list+list_len); + } +} + +template +struct MatrixInitOp { +}; + +template<> +struct MatrixInitOp > { + MatrixInitOp(const std::vector& rows_vec, + const std::vector& row_offsets_vec, + const std::vector& row_coords_vec, + int global_nx, int global_ny, int global_nz, + MINIFE_GLOBAL_ORDINAL global_n_rows, + const miniFE::simple_mesh_description& input_mesh, + miniFE::CSRMatrix& matrix) + : rows(&rows_vec[0]), + row_offsets(&row_offsets_vec[0]), + row_coords(&row_coords_vec[0]), + global_nodes_x(global_nx), + global_nodes_y(global_ny), + global_nodes_z(global_nz), + global_nrows(global_n_rows), + mesh(&input_mesh), + dest_rows(&matrix.rows[0]), + dest_rowoffsets(&matrix.row_offsets[0]), + dest_cols(&matrix.packed_cols[0]), + dest_coefs(&matrix.packed_coefs[0]), + n(matrix.rows.size()) + { + if (matrix.packed_cols.capacity() != matrix.packed_coefs.capacity()) { + std::cout<<"Warning, packed_cols.capacity ("<* mesh; + + inline void operator()(int i) + { + dest_rows[i] = rows[i]; + int offset = row_offsets[i]; + dest_rowoffsets[i] = offset; + int ix = row_coords[i*3]; + int iy = row_coords[i*3+1]; + int iz = row_coords[i*3+2]; + GlobalOrdinalType nnz = 0; + for(int sz=-1; sz<=1; ++sz) { + for(int sy=-1; sy<=1; ++sy) { + for(int sx=-1; sx<=1; ++sx) { + GlobalOrdinalType col_id = + miniFE::get_id(global_nodes_x, global_nodes_y, global_nodes_z, + ix+sx, iy+sy, iz+sz); + if (col_id >= 0 && col_id < global_nrows) { + GlobalOrdinalType col = mesh->map_id_to_row(col_id); + dest_cols[offset+nnz] = col; + dest_coefs[offset+nnz] = 0; + ++nnz; + } + } + } + } + + sort_if_needed(&dest_cols[offset], nnz); + } +}; + +template<> +struct MatrixInitOp > { + MatrixInitOp(const std::vector& rows_vec, + const std::vector& /*row_offsets_vec*/, + const std::vector& row_coords_vec, + int global_nx, int global_ny, int global_nz, + MINIFE_GLOBAL_ORDINAL global_n_rows, + const miniFE::simple_mesh_description& input_mesh, + miniFE::ELLMatrix& matrix) + : rows(&rows_vec[0]), + row_coords(&row_coords_vec[0]), + global_nodes_x(global_nx), + global_nodes_y(global_ny), + global_nodes_z(global_nz), + global_nrows(global_n_rows), + mesh(&input_mesh), + dest_rows(&matrix.rows[0]), + dest_cols(&matrix.cols[0]), + dest_coefs(&matrix.coefs[0]), + n(matrix.rows.size()), + ncols_per_row(matrix.num_cols_per_row) + { + } + + typedef MINIFE_GLOBAL_ORDINAL GlobalOrdinalType; + typedef MINIFE_LOCAL_ORDINAL LocalOrdinalType; + typedef MINIFE_SCALAR ScalarType; + + const GlobalOrdinalType* rows; + const int* row_coords; + + int global_nodes_x; + int global_nodes_y; + int global_nodes_z; + + GlobalOrdinalType global_nrows; + + GlobalOrdinalType* dest_rows; + GlobalOrdinalType* dest_cols; + ScalarType* dest_coefs; + int n; + int ncols_per_row; + + const miniFE::simple_mesh_description* mesh; + + inline void operator()(int i) + { + dest_rows[i] = rows[i]; + int offset = i*ncols_per_row; + int ix = row_coords[i*3]; + int iy = row_coords[i*3+1]; + int iz = row_coords[i*3+2]; + GlobalOrdinalType nnz = 0; + for(int sz=-1; sz<=1; ++sz) + for(int sy=-1; sy<=1; ++sy) + for(int sx=-1; sx<=1; ++sx) { + GlobalOrdinalType col_id = + miniFE::get_id(global_nodes_x, global_nodes_y, global_nodes_z, + ix+sx, iy+sy, iz+sz); + if (col_id >= 0 && col_id < global_nrows) { + GlobalOrdinalType col = mesh->map_id_to_row(col_id); + dest_cols[offset+nnz] = col; + dest_coefs[offset+nnz] = 0; + ++nnz; + } + } + + sort_if_needed(&dest_cols[offset], nnz); + } +}; + +#endif + diff --git a/tm/src/SparseMatrix_functions.hpp b/tm/src/SparseMatrix_functions.hpp new file mode 100644 index 0000000..87cb6a4 --- /dev/null +++ b/tm/src/SparseMatrix_functions.hpp @@ -0,0 +1,662 @@ +#ifndef _SparseMatrix_functions_hpp_ +#define _SparseMatrix_functions_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#ifdef MINIFE_HAVE_TBB +#include +#endif + +#ifdef HAVE_MPI +#include +#endif + +namespace miniFE { + +template +void init_matrix(MatrixType& M, + const std::vector& rows, + const std::vector& row_offsets, + const std::vector& row_coords, + int global_nodes_x, + int global_nodes_y, + int global_nodes_z, + typename MatrixType::GlobalOrdinalType global_nrows, + const simple_mesh_description& mesh) +{ + MatrixInitOp mat_init(rows, row_offsets, row_coords, + global_nodes_x, global_nodes_y, global_nodes_z, + global_nrows, mesh, M); + + for(int i=0; i +void sort_with_companions(ptrdiff_t len, T* array, U* companions) +{ + ptrdiff_t i, j, index; + U companion; + + for (i=1; i < len; i++) { + index = array[i]; + companion = companions[i]; + j = i; + while ((j > 0) && (array[j-1] > index)) + { + array[j] = array[j-1]; + companions[j] = companions[j-1]; + j = j - 1; + } + array[j] = index; + companions[j] = companion; + } +} + +template +void write_matrix(const std::string& filename, + MatrixType& mat) +{ + typedef typename MatrixType::LocalOrdinalType LocalOrdinalType; + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinalType; + typedef typename MatrixType::ScalarType ScalarType; + + int numprocs = 1, myproc = 0; +#ifdef HAVE_MPI + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myproc); +#endif + + std::ostringstream osstr; + osstr << filename << "." << numprocs << "." << myproc; + std::string full_name = osstr.str(); + std::ofstream ofs(full_name.c_str()); + + size_t nrows = mat.rows.size(); + size_t nnz = mat.num_nonzeros(); + + for(int p=0; p +void +sum_into_row(int row_len, + GlobalOrdinal* row_indices, + Scalar* row_coefs, + int num_inputs, + const GlobalOrdinal* input_indices, + const Scalar* input_coefs) +{ + while( true ) { + const int tm_status = _xbegin(); + + switch( tm_status ) { + case _XBEGIN_STARTED: + for(size_t i=0; i +void +sum_into_row(typename MatrixType::GlobalOrdinalType row, + size_t num_indices, + const typename MatrixType::GlobalOrdinalType* col_inds, + const typename MatrixType::ScalarType* coefs, + MatrixType& mat) +{ + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + typedef typename MatrixType::ScalarType Scalar; + + size_t row_len = 0; + GlobalOrdinal* mat_row_cols = NULL; + Scalar* mat_row_coefs = NULL; + + mat.get_row_pointers(row, row_len, mat_row_cols, mat_row_coefs); + if (row_len == 0) return; + + sum_into_row(row_len, mat_row_cols, mat_row_coefs, num_indices, col_inds, coefs); +} + +template +void +sum_in_symm_elem_matrix(size_t num, + const typename MatrixType::GlobalOrdinalType* indices, + const typename MatrixType::ScalarType* coefs, + MatrixType& mat) +{ + typedef typename MatrixType::ScalarType Scalar; + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + +//indices is length num (which should be nodes-per-elem) +//coefs is the upper triangle of the element diffusion matrix +//which should be length num*(num+1)/2 +//std::cout< +void +sum_in_elem_matrix(size_t num, + const typename MatrixType::GlobalOrdinalType* indices, + const typename MatrixType::ScalarType* coefs, + MatrixType& mat) +{ + size_t offset = 0; + + for(size_t i=0; i +void +sum_into_global_linear_system(ElemData& elem_data, + MatrixType& A, VectorType& b) +{ + sum_in_symm_elem_matrix(elem_data.nodes_per_elem, elem_data.elem_node_ids, + elem_data.elem_diffusion_matrix, A); + sum_into_vector(elem_data.nodes_per_elem, elem_data.elem_node_ids, + elem_data.elem_source_vector, b); +} + +#ifdef MINIFE_HAVE_TBB +template +void +sum_in_elem_matrix(size_t num, + const typename MatrixType::GlobalOrdinalType* indices, + const typename MatrixType::ScalarType* coefs, + LockingMatrix& mat) +{ + size_t offset = 0; + + for(size_t i=0; i +void +sum_into_global_linear_system(ElemData& elem_data, + LockingMatrix& A, LockingVector& b) +{ + sum_in_elem_matrix(elem_data.nodes_per_elem, elem_data.elem_node_ids, + elem_data.elem_diffusion_matrix, A); + sum_into_vector(elem_data.nodes_per_elem, elem_data.elem_node_ids, + elem_data.elem_source_vector, b); +} +#endif + +template +void +add_to_diagonal(typename MatrixType::ScalarType value, MatrixType& mat) +{ + for(size_t i=0; i +double +parallel_memory_overhead_MB(const MatrixType& A) +{ + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + typedef typename MatrixType::LocalOrdinalType LocalOrdinal; + double mem_MB = 0; + +#ifdef HAVE_MPI + double invMB = 1.0/(1024*1024); + mem_MB = invMB*A.external_index.size()*sizeof(GlobalOrdinal); + mem_MB += invMB*A.external_local_index.size()*sizeof(GlobalOrdinal); + mem_MB += invMB*A.elements_to_send.size()*sizeof(GlobalOrdinal); + mem_MB += invMB*A.neighbors.size()*sizeof(int); + mem_MB += invMB*A.recv_length.size()*sizeof(LocalOrdinal); + mem_MB += invMB*A.send_length.size()*sizeof(LocalOrdinal); + + double tmp = mem_MB; + MPI_Allreduce(&tmp, &mem_MB, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#endif + + return mem_MB; +} + +template +void rearrange_matrix_local_external(MatrixType& A) +{ + //This function will rearrange A so that local entries are contiguous at the front + //of A's memory, and external entries are contiguous at the back of A's memory. + // + //A.row_offsets will describe where the local entries occur, and + //A.row_offsets_external will describe where the external entries occur. + + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + typedef typename MatrixType::LocalOrdinalType LocalOrdinal; + typedef typename MatrixType::ScalarType Scalar; + + size_t nrows = A.rows.size(); + std::vector tmp_row_offsets(nrows*2); + std::vector tmp_row_offsets_external(nrows*2); + + LocalOrdinal num_local_nz = 0; + LocalOrdinal num_extern_nz = 0; + + //First sort within each row of A, so that local entries come + //before external entries within each row. + //tmp_row_offsets describe the locations of the local entries, and + //tmp_row_offsets_external describe the locations of the external entries. + // + for(size_t i=0; i ext_cols(num_extern_nz); + std::vector ext_coefs(num_extern_nz); + std::vector ext_offsets(nrows+1); + LocalOrdinal offset = 0; + for(size_t i=0; i +void +zero_row_and_put_1_on_diagonal(MatrixType& A, typename MatrixType::GlobalOrdinalType row) +{ + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + typedef typename MatrixType::LocalOrdinalType LocalOrdinal; + typedef typename MatrixType::ScalarType Scalar; + + size_t row_len = 0; + GlobalOrdinal* cols = NULL; + Scalar* coefs = NULL; + A.get_row_pointers(row, row_len, cols, coefs); + + for(size_t i=0; i +void +impose_dirichlet(typename MatrixType::ScalarType prescribed_value, + MatrixType& A, + VectorType& b, + int global_nx, + int global_ny, + int global_nz, + const std::set& bc_rows) +{ + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + typedef typename MatrixType::LocalOrdinalType LocalOrdinal; + typedef typename MatrixType::ScalarType Scalar; + + GlobalOrdinal first_local_row = A.rows.size()>0 ? A.rows[0] : 0; + GlobalOrdinal last_local_row = A.rows.size()>0 ? A.rows[A.rows.size()-1] : -1; + + typename std::set::const_iterator + bc_iter = bc_rows.begin(), bc_end = bc_rows.end(); + for(; bc_iter!=bc_end; ++bc_iter) { + GlobalOrdinal row = *bc_iter; + if (row >= first_local_row && row <= last_local_row) { + size_t local_row = row - first_local_row; + b.coefs[local_row] = prescribed_value; + zero_row_and_put_1_on_diagonal(A, row); + } + } + + const int ROW_COUNT = A.rows.size(); + + #pragma omp parallel for + for(MINIFE_GLOBAL_ORDINAL i=0; i < ROW_COUNT; ++i) { + GlobalOrdinal row = A.rows[i]; + + if (bc_rows.find(row) != bc_rows.end()) continue; + + size_t row_length = 0; + GlobalOrdinal* cols = NULL; + Scalar* coefs = NULL; + A.get_row_pointers(row, row_length, cols, coefs); + + Scalar sum = 0; + for(size_t j=0; j +struct matvec_std { +void operator()(MatrixType& A, + VectorType& x, + VectorType& y) +{ + exchange_externals(A, x); + + typedef typename MatrixType::ScalarType ScalarType; + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinalType; + typedef typename MatrixType::LocalOrdinalType LocalOrdinalType; + + const MINIFE_GLOBAL_ORDINAL rows_size = A.rows.size(); + const LocalOrdinalType* const Arowoffsets = &A.row_offsets[0]; + const GlobalOrdinalType* const Acols = &A.packed_cols[0]; + const ScalarType* const Acoefs = &A.packed_coefs[0]; + const ScalarType* const xcoefs = &x.coefs[0]; + ScalarType* ycoefs = &y.coefs[0]; + const ScalarType beta = 0; + + #pragma omp parallel for + for(MINIFE_GLOBAL_ORDINAL row = 0; row < rows_size; ++row) { + const MINIFE_GLOBAL_ORDINAL row_start = Arowoffsets[row]; + const MINIFE_GLOBAL_ORDINAL row_end = Arowoffsets[row+1]; + + MINIFE_SCALAR sum = 0; + + #pragma loop_count(15) + for(MINIFE_GLOBAL_ORDINAL i = row_start; i < row_end; ++i) { + sum += Acoefs[i] * xcoefs[Acols[i]]; + } + + ycoefs[row] = sum; + } +} +}; +#elif defined(MINIFE_ELL_MATRIX) +template +struct matvec_std { +void operator()(MatrixType& A, + VectorType& x, + VectorType& y) +{ + exchange_externals(A, x); + + typedef typename MatrixType::ScalarType ScalarType; + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinalType; + typedef typename MatrixType::LocalOrdinalType LocalOrdinalType; + + int row_len = A.num_cols_per_row; + int n = A.rows.size(); + const GlobalOrdinalType* Acols = &A.cols[0]; + const ScalarType* Acoefs = &A.coefs[0]; + const ScalarType* xcoefs = &x.coefs[0]; + ScalarType* ycoefs = &y.coefs[0]; + ScalarType beta = 0; + + #pragma omp parallel for + for(int row=0; row +void matvec(MatrixType& A, VectorType& x, VectorType& y) +{ + matvec_std mv; + mv(A, x, y); +} + +template +struct matvec_overlap { +void operator()(MatrixType& A, + VectorType& x, + VectorType& y) +{ +#ifdef HAVE_MPI + begin_exchange_externals(A, x); +#endif + + typedef typename MatrixType::ScalarType ScalarType; + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinalType; + typedef typename MatrixType::LocalOrdinalType LocalOrdinalType; + + + int n = A.rows.size(); + const LocalOrdinalType* Arowoffsets = &A.row_offsets[0]; + const GlobalOrdinalType* Acols = &A.packed_cols[0]; + const ScalarType* Acoefs = &A.packed_coefs[0]; + const ScalarType* xcoefs = &x.coefs[0]; + ScalarType* ycoefs = &y.coefs[0]; + ScalarType beta = 0; + + for(int row=0; row + +namespace miniFE { + + +template +struct Vector { + typedef Scalar ScalarType; + typedef LocalOrdinal LocalOrdinalType; + typedef GlobalOrdinal GlobalOrdinalType; + + Vector(GlobalOrdinal startIdx, LocalOrdinal local_sz) + : startIndex(startIdx), + local_size(local_sz), + coefs(local_size) + { + #pragma omp parallel for + for(MINIFE_LOCAL_ORDINAL i=0; i < local_size; ++i) { + coefs[i] = 0; + } + } + + ~Vector() + { + } + + GlobalOrdinal startIndex; + LocalOrdinal local_size; + std::vector coefs; +}; + + +}//namespace miniFE + +#endif + diff --git a/tm/src/Vector_functions.hpp b/tm/src/Vector_functions.hpp new file mode 100644 index 0000000..fd5f2f1 --- /dev/null +++ b/tm/src/Vector_functions.hpp @@ -0,0 +1,341 @@ +#ifndef _Vector_functions_hpp_ +#define _Vector_functions_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include + +#include +#include +#include + +#ifdef HAVE_MPI +#include +#endif + +#ifdef MINIFE_HAVE_TBB +#include +#endif + +#include +#include + +#define MINIFE_MIN(X, Y) ((X) < (Y) ? (X) : (Y)) + +namespace miniFE { + + +template +void write_vector(const std::string& filename, + const VectorType& vec) +{ + int numprocs = 1, myproc = 0; +#ifdef HAVE_MPI + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myproc); +#endif + + std::ostringstream osstr; + osstr << filename << "." << numprocs << "." << myproc; + std::string full_name = osstr.str(); + std::ofstream ofs(full_name.c_str()); + + typedef typename VectorType::ScalarType ScalarType; + + const std::vector& coefs = vec.coefs; + for(int p=0; p +void sum_into_vector(size_t num_indices, + const typename VectorType::GlobalOrdinalType* indices, + const typename VectorType::ScalarType* coefs, + VectorType& vec) +{ + typedef typename VectorType::GlobalOrdinalType GlobalOrdinal; + typedef typename VectorType::ScalarType Scalar; + + GlobalOrdinal first = vec.startIndex; + GlobalOrdinal last = first + vec.local_size - 1; + + std::vector& vec_coefs = vec.coefs; + + while( true ) { + const int tm_status = _xbegin(); + + switch( tm_status ) { + case _XBEGIN_STARTED: + for(size_t i=0; i last) + continue; + + size_t idx = indices[i] - first; + vec_coefs[idx] += coefs[i]; + } + + _xend(); + break; + case _XABORT_RETRY: + continue; + + default: + for(size_t i=0; i last) + continue; + + size_t idx = indices[i] - first; + #pragma omp atomic + vec_coefs[idx] += coefs[i]; + } + break; + } + + break; + } +} + +#ifdef MINIFE_HAVE_TBB +template +void sum_into_vector(size_t num_indices, + const typename VectorType::GlobalOrdinalType* indices, + const typename VectorType::ScalarType* coefs, + LockingVector& vec) +{ + vec.sum_in(num_indices, indices, coefs); +} +#endif + +//------------------------------------------------------------ +//Compute the update of a vector with the sum of two scaled vectors where: +// +// w = alpha*x + beta*y +// +// x,y - input vectors +// +// alpha,beta - scalars applied to x and y respectively +// +// w - output vector +// +template +void + waxpby(typename VectorType::ScalarType alpha, const VectorType& x, + typename VectorType::ScalarType beta, const VectorType& y, + VectorType& w) +{ + typedef typename VectorType::ScalarType ScalarType; + +#ifdef MINIFE_DEBUG_OPENMP + std::cout << "Starting WAXPBY..." << std::endl; +#endif + +#ifdef MINIFE_DEBUG + if (y.local_size < x.local_size || w.local_size < x.local_size) { + std::cerr << "miniFE::waxpby ERROR, y and w must be at least as long as x." << std::endl; + return; + } +#endif + + const int n = x.coefs.size(); + const ScalarType* xcoefs = &x.coefs[0]; + const ScalarType* ycoefs = &y.coefs[0]; + ScalarType* wcoefs = &w.coefs[0]; + + if(beta == 0.0) { + if(alpha == 1.0) { + #pragma omp parallel for + for(int i=0; i +void + daxpby(const MINIFE_SCALAR alpha, + const VectorType& x, + const MINIFE_SCALAR beta, + VectorType& y) +{ + + const MINIFE_LOCAL_ORDINAL n = MINIFE_MIN(x.coefs.size(), y.coefs.size()); + const MINIFE_SCALAR* xcoefs = &x.coefs[0]; + MINIFE_SCALAR* ycoefs = &y.coefs[0]; + + if(alpha == 1.0 && beta == 1.0) { + #pragma omp parallel for + for(int i = 0; i < n; ++i) { + ycoefs[i] += xcoefs[i]; + } + } else if (beta == 1.0) { + #pragma omp parallel for + for(int i = 0; i < n; ++i) { + ycoefs[i] += alpha * xcoefs[i]; + } + } else if (alpha == 1.0) { + #pragma omp parallel for + for(int i = 0; i < n; ++i) { + ycoefs[i] = xcoefs[i] + beta * ycoefs[i]; + } + } else if (beta == 0.0) { + #pragma omp parallel for + for(int i = 0; i < n; ++i) { + ycoefs[i] = alpha * xcoefs[i]; + } + } else { + #pragma omp parallel for + for(int i = 0; i < n; ++i) { + ycoefs[i] = alpha * xcoefs[i] + beta * ycoefs[i]; + } + } + +} + +//----------------------------------------------------------- +//Compute the dot product of two vectors where: +// +// x,y - input vectors +// +// result - return-value +// +template +typename TypeTraits::magnitude_type + dot(const Vector& x, + const Vector& y) +{ + const MINIFE_LOCAL_ORDINAL n = x.coefs.size(); + + typedef typename Vector::ScalarType Scalar; + typedef typename TypeTraits::magnitude_type magnitude; + + const Scalar* xcoefs = &x.coefs[0]; + const Scalar* ycoefs = &y.coefs[0]; + MINIFE_SCALAR result = 0; + + #pragma omp parallel for reduction(+:result) + for(int i=0; i::mpi_type(); + MPI_Allreduce(&local_dot, &global_dot, 1, mpi_dtype, MPI_SUM, MPI_COMM_WORLD); + return global_dot; +#else + return result; +#endif +} + +template +typename TypeTraits::magnitude_type + dot_r2(const Vector& x) +{ +#ifdef MINIFE_DEBUG_OPENMP + int myrank; + MPI_Comm_rank(MPI_COMM_WORLD, &myrank); + std::cout << "[" << myrank << "] Starting dot..." << std::endl; +#endif + + const MINIFE_LOCAL_ORDINAL n = x.coefs.size(); + +#ifdef MINIFE_DEBUG + if (y.local_size < n) { + std::cerr << "miniFE::dot ERROR, y must be at least as long as x."<::magnitude_type magnitude; + + const MINIFE_SCALAR* xcoefs = &x.coefs[0]; + MINIFE_SCALAR result = 0; + + #pragma omp parallel for reduction(+:result) + for(int i=0; i::mpi_type(); + MPI_Allreduce(&local_dot, &global_dot, 1, mpi_dtype, MPI_SUM, MPI_COMM_WORLD); + +#ifdef MINIFE_DEBUG_OPENMP + std::cout << "[" << myrank << "] Completed dot." << std::endl; +#endif + + return global_dot; +#else +#ifdef MINIFE_DEBUG_OPENMP + std::cout << "[" << myrank << "] Completed dot." << std::endl; +#endif + return result; +#endif +} + +}//namespace miniFE + +#endif + diff --git a/tm/src/YAML_Doc.cpp b/tm/src/YAML_Doc.cpp new file mode 100644 index 0000000..382c8cd --- /dev/null +++ b/tm/src/YAML_Doc.cpp @@ -0,0 +1,102 @@ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include +#include +#include +#include +#include +#include +#ifdef REDSTORM +#include +#include +#include +#endif +#include "YAML_Doc.hpp" +using namespace std; + +//set the microapp_name and version which will become part of the YAML doc. +YAML_Doc::YAML_Doc(const std::string& miniApp_Name, const std::string& miniApp_Version, const std::string& destination_Directory, const std::string& destination_FileName){ + miniAppName = miniApp_Name; + miniAppVersion = miniApp_Version; + destinationDirectory = destination_Directory; + destinationFileName = destination_FileName; +} + +//inherits the destructor from YAML_Element +YAML_Doc::~YAML_Doc(void){ +} + +/* +* generates YAML from the elements of the document and saves it +* to a file +*/ +string YAML_Doc::generateYAML(){ + string yaml; + yaml = yaml + "Mini-Application Name: " + miniAppName + "\n"; + yaml = yaml + "Mini-Application Version: " + miniAppVersion + "\n"; + for(size_t i=0; iprintYAML(""); + } + + time_t rawtime; + tm * ptm; + time ( &rawtime ); + ptm = localtime(&rawtime); + char sdate[25]; + //use tm_mon+1 because tm_mon is 0 .. 11 instead of 1 .. 12 + sprintf (sdate,"%04d:%02d:%02d-%02d:%02d:%02d",ptm->tm_year + 1900, ptm->tm_mon+1, + ptm->tm_mday, ptm->tm_hour, ptm->tm_min,ptm->tm_sec); + + string filename; + if (destinationFileName=="") + filename = miniAppName + "-" + miniAppVersion + "_"; + else + filename = destinationFileName; + filename = filename + string(sdate) + ".yaml"; + if (destinationDirectory!="" && destinationDirectory!=".") { + string mkdir_cmd = "mkdir " + destinationDirectory; +#ifdef REDSTORM + mkdir(destinationDirectory.c_str(),0755); +#else + system(mkdir_cmd.c_str()); +#endif + filename = destinationDirectory + "/" + destinationFileName; + } + else + filename = "./" + filename; + + ofstream myfile; + myfile.open(filename.c_str()); + myfile << yaml; + myfile.close(); + return yaml; +} + + diff --git a/tm/src/YAML_Doc.hpp b/tm/src/YAML_Doc.hpp new file mode 100644 index 0000000..8ab6961 --- /dev/null +++ b/tm/src/YAML_Doc.hpp @@ -0,0 +1,115 @@ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#ifndef YAML_DOC_H +#define YAML_DOC_H +#include +#include +#include "YAML_Element.hpp" + +//! The Mantevo YAML_Doc class for the uniform collecting and reporting of performance data for mini-applications + +/*! + +The YAML_Doc class works in conjuction with the YAML_Element class to facilitate easy collecting and reporting of YAML-formatted +data that can be then registered with the Mantevo results collection website. + +\code + +//EXAMPLE CODE FOR GENERATING YAML + + YAML_Doc doc("hpccg","1.0"); + doc.add("final_residual",1.4523e-13); + doc.add("time","4.893"); + +//note: the following line will remove the data (4.890) associated with "time" + doc.get("time")->add("total",4.243); + +//note: the following line will likewise remove the data (1.243) associated with "time" + doc.get("time")->get("total")->add("time",2.457); + doc.get("time")->get("total")->add("flops",4.88e5); + doc.get("time")->add("ddot",1.243); + doc.get("time")->add("sparsemv",""); + doc.get("time")->get("sparsemv")->add("time",0.3445); + doc.get("time")->get("sparsemv")->add("overhead",""); + doc.get("time")->get("sparsemv")->get("overhead")->add("time",0.0123); + doc.get("time")->get("sparsemv")->get("overhead")->add("percentage",0.034); + cout << doc.generateYAML() << endl; + return 0; + +\endcode + +Below is the output generated by the above code: + +\verbatim + +final_residual: 1.4523e-13 +time: + total: + time: 2.457 + flops: 4.88e5 + ddot: 1.243 + sparsemv: + time: 0.3445 + overhead: + time: 0.0123 + percentage: 0.034 + +\endverbatim + +\note {No value is allowed to be attached to a key that has children. If children are added to a key, the value is simply set to "".} + +*/ +class YAML_Doc: public YAML_Element { + public: + //! Constructor: accepts mini-application name and version as strings, optionally accepts directory and file name for printing results. + /*! + The sole constructor for this class accepts and name and version number for the mini-application as well as optional directory + and file name information for results that are generated by the generateYAML() method. + \param miniApp_Name (in) string containing name of the mini-application + \param miniApp_Version (in) string containing the version of the mini-application + \param destination_Directory (in, optional) path of diretory where results file will be stored, relative to current working directory. + If this value is not supplied, the results file will be stored in the current working directory. If the directory does not exist + it will be created. + \param destination_FileName (in, optional) root name of the results file. A suffix of ".yaml" will be automatically appended. If no + file name is specified the filename will be constructed by concatenating the miniAppName + miniAppVersion + ".yaml" strings. + */ + YAML_Doc(const std::string& miniApp_Name, const std::string& miniApp_Version, const std::string& destination_Directory = "", const std::string& destination_FileName = ""); + //! Destructor + ~YAML_Doc(); + //! Generate YAML results to standard out and to a file using specified directory and filename, using current directory and miniAppName + miniAppVersion + ".yaml" by default + std::string generateYAML(); + +protected: + std::string miniAppName; + std::string miniAppVersion; + std::string destinationDirectory; + std::string destinationFileName; +}; +#endif /* YAML_DOC_H */ + diff --git a/tm/src/YAML_Element.cpp b/tm/src/YAML_Element.cpp new file mode 100644 index 0000000..26ba04b --- /dev/null +++ b/tm/src/YAML_Element.cpp @@ -0,0 +1,148 @@ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include +#include +#include +#include "YAML_Element.hpp" +using namespace std; +YAML_Element::YAML_Element(const std::string& key_arg, const std::string& value_arg){ + key = key_arg; + value = value_arg; +} + +YAML_Element::~YAML_Element(){ + for (size_t i=0; ivalue = ""; + string converted_value = convert_double_to_string(value_arg); + YAML_Element* element = new YAML_Element(key_arg,converted_value); + children.push_back(element); + return element; +} + +YAML_Element* YAML_Element::add(const std::string& key_arg, int value_arg) { + this->value = ""; + string converted_value = convert_int_to_string(value_arg); + YAML_Element* element = new YAML_Element(key_arg,converted_value); + children.push_back(element); + return element; +} + +#ifndef MINIFE_NO_LONG_LONG + +YAML_Element* YAML_Element::add(const std::string& key_arg, long long value_arg) { + this->value = ""; + string converted_value = convert_long_long_to_string(value_arg); + YAML_Element* element = new YAML_Element(key_arg,converted_value); + children.push_back(element); + return element; +} + +#endif + +YAML_Element* YAML_Element::add(const std::string& key_arg, size_t value_arg) { + this->value = ""; + string converted_value = convert_size_t_to_string(value_arg); + YAML_Element* element = new YAML_Element(key_arg,converted_value); + children.push_back(element); + return element; +} + +YAML_Element* YAML_Element::add(const std::string& key_arg, const std::string& value_arg) { + this->value = ""; + YAML_Element* element = new YAML_Element(key_arg, value_arg); + children.push_back(element); + return element; +} + +/* +* returns pointer to the YAML_Element for the given key. +* I, cam, believe an exception should be thrown if there is no +* element in the vector for the specified key +*/ +YAML_Element* YAML_Element::get(const std::string& key_arg) { + for (size_t i=0; igetKey() == key_arg){ + return children[i]; + } + } + return 0; +} + +/* +* prints a line of a YAML document. Correct YAML depends on +* correct spacing; the parameter space should be the proper +* amount of space for the parent element +*/ +string YAML_Element::printYAML(std::string space){ + string yaml_line = space + key + ": " + value + "\n"; + for(int i=0; i<2; i++) space = space + " "; + for(size_t i=0; iprintYAML(space); + } + return yaml_line; +} + +string YAML_Element::convert_double_to_string(double value_arg){ + stringstream strm; + strm << value_arg; + return strm.str(); +} +string YAML_Element::convert_int_to_string(int value_arg){ + stringstream strm; + strm << value_arg; + return strm.str(); +} + +#ifndef MINIFE_NO_LONG_LONG + +string YAML_Element::convert_long_long_to_string(long long value_arg){ + stringstream strm; + strm << value_arg; + return strm.str(); +} + +#endif + +string YAML_Element::convert_size_t_to_string(size_t value_arg){ + stringstream strm; + strm << value_arg; + return strm.str(); +} diff --git a/tm/src/YAML_Element.hpp b/tm/src/YAML_Element.hpp new file mode 100644 index 0000000..8ab1e5f --- /dev/null +++ b/tm/src/YAML_Element.hpp @@ -0,0 +1,79 @@ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#ifndef YAML_ELEMENT_H +#define YAML_ELEMENT_H +#include +#include +//! The Mantevo YAML_Element class for registering key-value pairs of performance data + +/*! + Mantevo mini-applications generate a collection of performance data for each run of the executable. YAML_Element, and + the related YAML_Doc class, provide a uniform facility for gathering and reporting this data using the YAML text format. +*/ +class YAML_Element { + public: + + //! Default constructor. + YAML_Element (){key="";value="";} + //! Construct with known key-value pair + YAML_Element (const std::string& key_arg, const std::string& value_arg); + //! Destructor + ~YAML_Element (); + //! Key accessor method + std::string getKey(){return key;} + //! Add a child element to an element list associated with this element, value of type double + YAML_Element* add(const std::string& key_arg, double value_arg); + //! Add a child element to an element list associated with this element, value of type int + YAML_Element* add(const std::string& key_arg, int value_arg); +#ifndef MINIFE_NO_LONG_LONG + //! Add a child element to an element list associated with this element, value of type long long + YAML_Element* add(const std::string& key_arg, long long value_arg); +#endif + //! Add a child element to an element list associated with this element, value of type size_t + YAML_Element* add(const std::string& key_arg, size_t value_arg); + //! Add a child element to an element list associated with this element, value of type string + YAML_Element* add(const std::string& key_arg, const std::string& value_arg); + //! get the element in the list with the given key + YAML_Element* get(const std::string& key_arg); + std::string printYAML(std::string space); + +protected: + std::string key; + std::string value; + std::vector children; + +private: + std::string convert_double_to_string(double value_arg); + std::string convert_int_to_string(int value_arg); +#ifndef MINIFE_NO_LONG_LONG + std::string convert_long_long_to_string(long long value_arg); +#endif + std::string convert_size_t_to_string(size_t value_arg); +}; +#endif /* YAML_ELEMENT_H */ diff --git a/tm/src/assemble_FE_data.hpp b/tm/src/assemble_FE_data.hpp new file mode 100644 index 0000000..a60b29c --- /dev/null +++ b/tm/src/assemble_FE_data.hpp @@ -0,0 +1,78 @@ +#ifndef _assemble_FE_data_hpp_ +#define _assemble_FE_data_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include + +#include + +namespace miniFE { + +template +void +assemble_FE_data(const simple_mesh_description& mesh, + MatrixType& A, + VectorType& b, + Parameters& params) +{ + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + + int global_elems_x = mesh.global_box[0][1]; + int global_elems_y = mesh.global_box[1][1]; + int global_elems_z = mesh.global_box[2][1]; + + Box local_elem_box; + copy_box(mesh.local_box, local_elem_box); + + if (get_num_ids(local_elem_box) < 1) { + return; + } + + // + //We want the element-loop to loop over our (processor-local) domain plus a + //ghost layer, so we can assemble the complete linear-system without doing + //any communication. + // + int ghost = 1; + if (local_elem_box[0][0] > 0) local_elem_box[0][0] -= ghost; + if (local_elem_box[1][0] > 0) local_elem_box[1][0] -= ghost; + if (local_elem_box[2][0] > 0) local_elem_box[2][0] -= ghost; + if (local_elem_box[0][1] < global_elems_x) local_elem_box[0][1] += ghost; + if (local_elem_box[1][1] < global_elems_y) local_elem_box[1][1] += ghost; + if (local_elem_box[2][1] < global_elems_z) local_elem_box[2][1] += ghost; + + perform_element_loop(mesh, local_elem_box, A, b, params); +} + +}//namespace miniFE + +#endif + diff --git a/tm/src/cg_solve.hpp b/tm/src/cg_solve.hpp new file mode 100644 index 0000000..69812c4 --- /dev/null +++ b/tm/src/cg_solve.hpp @@ -0,0 +1,215 @@ +#ifndef _cg_solve_hpp_ +#define _cg_solve_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include + +#include +#include + +#include + +namespace miniFE { + +template +void print_vec(const std::vector& vec, const std::string& name) +{ + for(size_t i=0; i +bool breakdown(typename VectorType::ScalarType inner, + const VectorType& v, + const VectorType& w) +{ + typedef typename VectorType::ScalarType Scalar; + typedef typename TypeTraits::magnitude_type magnitude; + +//This is code that was copied from Aztec, and originally written +//by my hero, Ray Tuminaro. +// +//Assuming that inner = (inner product of v and w), +//v and w are considered orthogonal if +// |inner| < 100 * ||v||_2 * ||w||_2 * epsilon + + magnitude vnorm = std::sqrt(dot(v,v)); + magnitude wnorm = std::sqrt(dot(w,w)); + return std::abs(inner) <= 100*vnorm*wnorm*std::numeric_limits::epsilon(); +} + +template +void +cg_solve(OperatorType& A, + const VectorType& b, + VectorType& x, + Matvec matvec, + typename OperatorType::LocalOrdinalType max_iter, + typename TypeTraits::magnitude_type& tolerance, + typename OperatorType::LocalOrdinalType& num_iters, + typename TypeTraits::magnitude_type& normr, + timer_type* my_cg_times) +{ + typedef typename OperatorType::ScalarType ScalarType; + typedef typename OperatorType::GlobalOrdinalType GlobalOrdinalType; + typedef typename OperatorType::LocalOrdinalType LocalOrdinalType; + typedef typename TypeTraits::magnitude_type magnitude_type; + + timer_type t0 = 0, tWAXPY = 0, tDOT = 0, tMATVEC = 0, tMATVECDOT = 0; + timer_type total_time = mytimer(); + + int myproc = 0; +#ifdef HAVE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &myproc); +#endif + + if (!A.has_local_indices) { + std::cerr << "miniFE::cg_solve ERROR, A.has_local_indices is false, needs to be true. This probably means " + << "miniFE::make_local_matrix(A) was not called prior to calling miniFE::cg_solve." + << std::endl; + return; + } + + size_t nrows = A.rows.size(); + LocalOrdinalType ncols = A.num_cols; + + VectorType r(b.startIndex, nrows); + VectorType p(0, ncols); + VectorType Ap(b.startIndex, nrows); + + normr = 0; + magnitude_type rtrans = 0; + magnitude_type oldrtrans = 0; + + LocalOrdinalType print_freq = max_iter/10; + if (print_freq>50) print_freq = 50; + if (print_freq<1) print_freq = 1; + + ScalarType one = 1.0; + ScalarType zero = 0.0; + + TICK(); waxpby(one, x, zero, x, p); TOCK(tWAXPY); + +// print_vec(p.coefs, "p"); + + TICK(); + matvec(A, p, Ap); + TOCK(tMATVEC); + + TICK(); waxpby(one, b, -one, Ap, r); TOCK(tWAXPY); + + TICK(); rtrans = dot_r2(r); TOCK(tDOT); + +//std::cout << "rtrans="< tolerance; ++k) { + if (k == 1) { + TICK(); waxpby(one, r, zero, r, p); TOCK(tWAXPY); + } + else { + oldrtrans = rtrans; + TICK(); rtrans = dot_r2(r); TOCK(tDOT); + magnitude_type beta = rtrans/oldrtrans; + TICK(); daxpby(one, r, beta, p); TOCK(tWAXPY); + } + + normr = sqrt(rtrans); + + if (myproc == 0 && (k%print_freq==0 || k==max_iter)) { + std::cout << "Iteration = "< +int +driver(const Box& global_box, Box& my_box, + Parameters& params, YAML_Doc& ydoc) +{ + int global_nx = global_box[0][1]; + int global_ny = global_box[1][1]; + int global_nz = global_box[2][1]; + + int numprocs = 1, myproc = 0; +#ifdef HAVE_MPI + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myproc); +#endif + + if (params.load_imbalance > 0) { + add_imbalance(global_box, my_box, params.load_imbalance, ydoc); + } + + float largest_imbalance = 0, std_dev = 0; + compute_imbalance(global_box, my_box, largest_imbalance, + std_dev, ydoc, true); + + + //Create a representation of the mesh: + //Note that 'simple_mesh_description' is a virtual or conceptual + //mesh that doesn't actually store mesh data. + + if (myproc==0) { + std::cout.width(30); + std::cout << "creating/filling mesh..."; + std::cout.flush(); + } + + timer_type t_start = mytimer(); + timer_type t0 = mytimer(); + + simple_mesh_description mesh(global_box, my_box); + + timer_type mesh_fill = mytimer() - t0; + timer_type t_total = mytimer() - t_start; + + if (myproc==0) { + std::cout << mesh_fill << "s, total time: " << t_total << std::endl; + } + + //next we will generate the matrix structure. + + //Declare matrix object: + +#if defined(MINIFE_ELL_MATRIX) + typedef ELLMatrix MatrixType; +#else + typedef CSRMatrix MatrixType; +#endif + + MatrixType A; + + timer_type gen_structure; + RUN_TIMED_FUNCTION("generating matrix structure...", + generate_matrix_structure(mesh, A), + gen_structure, t_total); + + GlobalOrdinal local_nrows = A.rows.size(); + GlobalOrdinal my_first_row = local_nrows > 0 ? A.rows[0] : -1; + + Vector b(my_first_row, local_nrows); + Vector x(my_first_row, local_nrows); + + //Assemble finite-element sub-matrices and sub-vectors into the global + //linear system: + + timer_type fe_assembly; + RUN_TIMED_FUNCTION("assembling FE data...", + assemble_FE_data(mesh, A, b, params), + fe_assembly, t_total); + + if (myproc == 0) { + ydoc.add("Matrix structure generation",""); + ydoc.get("Matrix structure generation")->add("Mat-struc-gen Time",gen_structure); + ydoc.add("FE assembly",""); + ydoc.get("FE assembly")->add("FE assembly Time",fe_assembly); + } + +#ifdef MINIFE_DEBUG + write_matrix("A_prebc.mtx", A); + write_vector("b_prebc.vec", b); +#endif + + //Now apply dirichlet boundary-conditions + //(Apply the 0-valued surfaces first, then the 1-valued surface last.) + + timer_type dirbc_time; + RUN_TIMED_FUNCTION("imposing Dirichlet BC...", + impose_dirichlet(0.0, A, b, global_nx+1, global_ny+1, global_nz+1, mesh.bc_rows_0), dirbc_time, t_total); + RUN_TIMED_FUNCTION("imposing Dirichlet BC...", + impose_dirichlet(1.0, A, b, global_nx+1, global_ny+1, global_nz+1, mesh.bc_rows_1), dirbc_time, t_total); + +#ifdef MINIFE_DEBUG + write_matrix("A.mtx", A); + write_vector("b.vec", b); +#endif + + //Transform global indices to local, set up communication information: + + timer_type make_local_time; + RUN_TIMED_FUNCTION("making matrix indices local...", + make_local_matrix(A), + make_local_time, t_total); + +#ifdef MINIFE_DEBUG + write_matrix("A_local.mtx", A); + write_vector("b_local.vec", b); +#endif + + size_t global_nnz = compute_matrix_stats(A, myproc, numprocs, ydoc); + + //Prepare to perform conjugate gradient solve: + + LocalOrdinal max_iters = 200; + LocalOrdinal num_iters = 0; + typedef typename TypeTraits::magnitude_type magnitude; + magnitude rnorm = 0; + magnitude tol = std::numeric_limits::epsilon(); + + timer_type cg_times[NUM_TIMERS]; + + typedef Vector VectorType; + + t_total = mytimer() - t_start; + + bool matvec_with_comm_overlap = params.mv_overlap_comm_comp==1; + + int verify_result = 0; + +#if MINIFE_KERNELS != 0 + if (myproc==0) { + std::cout.width(30); + std::cout << "Starting kernel timing loops ..." << std::endl; + } + + max_iters = 500; + x.coefs[0] = 0.9; + if (matvec_with_comm_overlap) { + time_kernels(A, b, x, matvec_overlap(), max_iters, rnorm, cg_times); + } + else { + time_kernels(A, b, x, matvec_std(), max_iters, rnorm, cg_times); + } + num_iters = max_iters; + std::string title("Kernel timings"); +#else + if (myproc==0) { + std::cout << "Starting CG solver ... " << std::endl; + } + + if (matvec_with_comm_overlap) { +#ifdef MINIFE_CSR_MATRIX + rearrange_matrix_local_external(A); + cg_solve(A, b, x, matvec_overlap(), max_iters, tol, + num_iters, rnorm, cg_times); +#else + std::cout << "ERROR, matvec with overlapping comm/comp only works with CSR matrix."<(), max_iters, tol, + num_iters, rnorm, cg_times); + if (myproc == 0) { + std::cout << "Final Resid Norm: " << rnorm << std::endl; + } + + if (params.verify_solution > 0) { + double tolerance = 0.06; + bool verify_whole_domain = false; + #ifdef MINIFE_DEBUG + verify_whole_domain = true; + #endif + if (myproc == 0) { + if (verify_whole_domain) std::cout << "verifying solution..." << std::endl; + else std::cout << "verifying solution at ~ (0.5, 0.5, 0.5) ..." << std::endl; + } + verify_result = verify_solution(mesh, x, tolerance, verify_whole_domain); + } + } + +#ifdef MINIFE_DEBUG + write_vector("x.vec", x); +#endif + std::string title("CG solve"); +#endif + + if (myproc == 0) { + ydoc.get("Global Run Parameters")->add("ScalarType",TypeTraits::name()); + ydoc.get("Global Run Parameters")->add("GlobalOrdinalType",TypeTraits::name()); + ydoc.get("Global Run Parameters")->add("LocalOrdinalType",TypeTraits::name()); + ydoc.add(title,""); + ydoc.get(title)->add("Iterations",num_iters); + ydoc.get(title)->add("Final Resid Norm",rnorm); + + GlobalOrdinal global_nrows = global_nx; + global_nrows *= global_ny*global_nz; + + //flops-per-mv, flops-per-dot, flops-per-waxpy: + double mv_flops = global_nnz*2.0; + double dot_flops = global_nrows*2.0; + double waxpy_flops = global_nrows*3.0; + +#if MINIFE_KERNELS == 0 +//if MINIFE_KERNELS == 0 then we did a CG solve, and in that case +//there were num_iters+1 matvecs, num_iters*2 dots, and num_iters*3+2 waxpys. + mv_flops *= (num_iters+1); + dot_flops *= (2*num_iters); + waxpy_flops *= (3*num_iters+2); +#else +//if MINIFE_KERNELS then we did one of each operation per iteration. + mv_flops *= num_iters; + dot_flops *= num_iters; + waxpy_flops *= num_iters; +#endif + + double total_flops = mv_flops + dot_flops + waxpy_flops; + + double mv_mflops = -1; + if (cg_times[MATVEC] > 1.e-4) + mv_mflops = 1.e-6 * (mv_flops/cg_times[MATVEC]); + + double dot_mflops = -1; + if (cg_times[DOT] > 1.e-4) + dot_mflops = 1.e-6 * (dot_flops/cg_times[DOT]); + + double waxpy_mflops = -1; + if (cg_times[WAXPY] > 1.e-4) + waxpy_mflops = 1.e-6 * (waxpy_flops/cg_times[WAXPY]); + + double total_mflops = -1; + if (cg_times[TOTAL] > 1.e-4) + total_mflops = 1.e-6 * (total_flops/cg_times[TOTAL]); + + ydoc.get(title)->add("WAXPY Time",cg_times[WAXPY]); + ydoc.get(title)->add("WAXPY Flops",waxpy_flops); + if (waxpy_mflops >= 0) + ydoc.get(title)->add("WAXPY Mflops",waxpy_mflops); + else + ydoc.get(title)->add("WAXPY Mflops","inf"); + + ydoc.get(title)->add("DOT Time",cg_times[DOT]); + ydoc.get(title)->add("DOT Flops",dot_flops); + if (dot_mflops >= 0) + ydoc.get(title)->add("DOT Mflops",dot_mflops); + else + ydoc.get(title)->add("DOT Mflops","inf"); + + ydoc.get(title)->add("MATVEC Time",cg_times[MATVEC]); + ydoc.get(title)->add("MATVEC Flops",mv_flops); + if (mv_mflops >= 0) + ydoc.get(title)->add("MATVEC Mflops",mv_mflops); + else + ydoc.get(title)->add("MATVEC Mflops","inf"); + +#ifdef MINIFE_FUSED + ydoc.get(title)->add("MATVECDOT Time",cg_times[MATVECDOT]); + ydoc.get(title)->add("MATVECDOT Flops",mv_flops); + if (mv_mflops >= 0) + ydoc.get(title)->add("MATVECDOT Mflops",mv_mflops); + else + ydoc.get(title)->add("MATVECDOT Mflops","inf"); +#endif + +#if MINIFE_KERNELS == 0 + ydoc.get(title)->add("Total",""); + ydoc.get(title)->get("Total")->add("Total CG Time",cg_times[TOTAL]); + ydoc.get(title)->get("Total")->add("Total CG Flops",total_flops); + if (total_mflops >= 0) + ydoc.get(title)->get("Total")->add("Total CG Mflops",total_mflops); + else + ydoc.get(title)->get("Total")->add("Total CG Mflops","inf"); + ydoc.get(title)->add("Time per iteration",cg_times[TOTAL]/num_iters); +#endif + } + + return verify_result; +} + +}//namespace miniFE + +#endif + diff --git a/tm/src/exchange_externals.hpp b/tm/src/exchange_externals.hpp new file mode 100644 index 0000000..7faf242 --- /dev/null +++ b/tm/src/exchange_externals.hpp @@ -0,0 +1,270 @@ +#ifndef _exchange_externals_hpp_ +#define _exchange_externals_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include + +#ifdef HAVE_MPI +#include +#endif + +#include + +#include + +namespace miniFE { + +template +void +exchange_externals(MatrixType& A, + VectorType& x) +{ +#ifdef HAVE_MPI +#ifdef MINIFE_DEBUG + std::ostream& os = outstream(); + os << "entering exchange_externals\n"; +#endif + + int numprocs = 1; + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + + if (numprocs < 2) return; + + typedef typename MatrixType::ScalarType Scalar; + typedef typename MatrixType::LocalOrdinalType LocalOrdinal; + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + + // Extract Matrix pieces + + int local_nrow = A.rows.size(); + int num_neighbors = A.neighbors.size(); + const std::vector& recv_length = A.recv_length; + const std::vector& send_length = A.send_length; + const std::vector& neighbors = A.neighbors; + const std::vector& elements_to_send = A.elements_to_send; + + std::vector& send_buffer = A.send_buffer; + + // + // first post receives, these are immediate receives + // Do not wait for result to come, will do that at the + // wait call below. + // + + int MPI_MY_TAG = 99; + + std::vector& request = A.request; + + // + // Externals are at end of locals + // + + std::vector& x_coefs = x.coefs; + Scalar* x_external = &(x_coefs[local_nrow]); + + MPI_Datatype mpi_dtype = TypeTraits::mpi_type(); + + // Post receives first + for(int i=0; i x.coefs.size()) { + os << "error, out-of-range. x.coefs.size()=="< exch_ext_requests; +#endif + +template +void +begin_exchange_externals(MatrixType& A, + VectorType& x) +{ +#ifdef HAVE_MPI + + int numprocs = 1, myproc = 0; + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myproc); + + if (numprocs < 2) return; + + typedef typename MatrixType::ScalarType Scalar; + typedef typename MatrixType::LocalOrdinalType LocalOrdinal; + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + + // Extract Matrix pieces + + int local_nrow = A.rows.size(); + int num_neighbors = A.neighbors.size(); + const std::vector& recv_length = A.recv_length; + const std::vector& send_length = A.send_length; + const std::vector& neighbors = A.neighbors; + const std::vector& elements_to_send = A.elements_to_send; + + std::vector send_buffer(elements_to_send.size(), 0); + + // + // first post receives, these are immediate receives + // Do not wait for result to come, will do that at the + // wait call below. + // + + int MPI_MY_TAG = 99; + + exch_ext_requests.resize(num_neighbors); + + // + // Externals are at end of locals + // + + std::vector& x_coefs = x.coefs; + Scalar* x_external = &(x_coefs[local_nrow]); + + MPI_Datatype mpi_dtype = TypeTraits::mpi_type(); + + // Post receives first + for(int i=0; i +# example: +# % generate_info_header g++ -O3 miniFE MINIFE +# this will cause the appropriate info to be put in a +# header named miniFE_info.hpp and the info will be in macros +# that start with MINIFE. +# +# an example of usage can be seen in miniFE/make_targets +# +if [ $# != 4 ] ; then +echo "error, need 4 arguments."; +exit 1; +fi + +cxx=`which ${1}` +errcode="$?" +if [ ${errcode} != "0" ] ; then +cxx="unknown"; +fi +echo "CXX: ${cxx}" + +cxx_ver=`${1} --version 2>&1` +errcode="$?" +if [ ${errcode} != "0" ] ; then +cxx_ver=`${1} -V 2>&1`; +errcode="$?" +if [ ${errcode} != "0" ] ; then +cxx_ver="unknown"; +fi +fi + +cxx_ver=${cxx_ver// /@} +cxx_version="" +for i in $(echo ${cxx_ver}); +do + if [ "$cxx_version" == "" ]; then + cxx_version=$i; + fi +done +cxx_version=${cxx_version//@/ } +echo "Compiler version: ${cxx_version}" + +cxxflags=${2} +hostname=`uname -n` +errcode="$?" +if [ ${errcode} != "0" ] ; then +hostname="unknown"; +fi + +kern_name=`uname -s` +errcode="$?" +if [ ${errcode} != "0" ] ; then +kern_name="unknown"; +fi + +kern_rel=`uname -r` +errcode="$?" +if [ ${errcode} != "0" ] ; then +kern_rel="unknown"; +fi + +proc=`uname -p` +errcode="$?" +if [ ${errcode} != "0" ] ; then +proc="unknown"; +fi + +header_prefix=${3} +macro_prefix=${4} + +cat << END_CAT > ${header_prefix}_info.hpp +#ifndef ${header_prefix}_info_hpp +#define ${header_prefix}_info_hpp + +#define ${macro_prefix}_HOSTNAME "${hostname}" +#define ${macro_prefix}_KERNEL_NAME "'${kern_name}'" +#define ${macro_prefix}_KERNEL_RELEASE "'${kern_rel}'" +#define ${macro_prefix}_PROCESSOR "'${proc}'" + +#define ${macro_prefix}_CXX "'${cxx}'" +#define ${macro_prefix}_CXX_VERSION "'${cxx_version}'" +#define ${macro_prefix}_CXXFLAGS "'${cxxflags}'" + +#endif +END_CAT diff --git a/tm/src/generate_matrix_structure.hpp b/tm/src/generate_matrix_structure.hpp new file mode 100644 index 0000000..76de40c --- /dev/null +++ b/tm/src/generate_matrix_structure.hpp @@ -0,0 +1,165 @@ +#ifndef _generate_matrix_structure_hpp_ +#define _generate_matrix_structure_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include +#include +#include + +#include +#include +#include +#include + +#ifdef HAVE_MPI +#include +#endif + +namespace miniFE { + +template +int +generate_matrix_structure(const simple_mesh_description& mesh, + MatrixType& A) +{ + int myproc = 0; +#ifdef HAVE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &myproc); +#endif + + int threw_exc = 0; + try { + + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + typedef typename MatrixType::LocalOrdinalType LocalOrdinal; + + const int global_nodes_x = mesh.global_box[0][1]+1; + const int global_nodes_y = mesh.global_box[1][1]+1; + const int global_nodes_z = mesh.global_box[2][1]+1; + Box box; + copy_box(mesh.local_box, box); + + //num-owned-nodes in each dimension is num-elems+1 + //only if num-elems > 0 in that dimension *and* + //we are at the high end of the global range in that dimension: + if (box[0][1] > box[0][0] && box[0][1] == mesh.global_box[0][1]) ++box[0][1]; + if (box[1][1] > box[1][0] && box[1][1] == mesh.global_box[1][1]) ++box[1][1]; + if (box[2][1] > box[2][0] && box[2][1] == mesh.global_box[2][1]) ++box[2][1]; + + GlobalOrdinal global_nrows = global_nodes_x; + global_nrows *= global_nodes_y*global_nodes_z; + + GlobalOrdinal nrows = get_num_ids(box); + try { + A.reserve_space(nrows, 27); + } + catch(std::exception& exc) { + std::ostringstream osstr; + osstr << "One of A.rows.resize, A.row_offsets.resize, A.packed_cols.reserve or A.packed_coefs.reserve: nrows=" < rows(nrows); + std::vector row_offsets(nrows+1); + std::vector row_coords(nrows*3); + + const MINIFE_GLOBAL_ORDINAL z_width = box[2][1] - box[2][0]; + const MINIFE_GLOBAL_ORDINAL y_width = box[1][1] - box[1][0]; + const MINIFE_GLOBAL_ORDINAL x_width = box[0][1] - box[0][0]; + const MINIFE_GLOBAL_ORDINAL r_n = (box[2][1] - box[2][0]) * + (box[1][1] - box[1][0]) * + (box[0][1] - box[0][0]); + const MINIFE_GLOBAL_ORDINAL xy_width = x_width * y_width; + MINIFE_GLOBAL_ORDINAL* const row_ptr = &rows[0]; + MINIFE_LOCAL_ORDINAL* const row_offset_ptr = &row_offsets[0]; + MINIFE_LOCAL_ORDINAL* const row_coords_ptr = &row_coords[0]; + + #pragma omp parallel for + for(int r = 0; r < r_n; ++r) { + int iz = r / (xy_width) + box[2][0]; + int iy = (r / x_width) % y_width + box[1][0]; + int ix = r % x_width + box[0][0]; + + GlobalOrdinal row_id = + get_id(global_nodes_x, global_nodes_y, global_nodes_z, + ix, iy, iz); + row_ptr[r] = mesh.map_id_to_row(row_id); + row_coords_ptr[r*3] = ix; + row_coords_ptr[r*3+1] = iy; + row_coords_ptr[r*3+2] = iz; + + MINIFE_LOCAL_ORDINAL nnz = 0; + for(int sz=-1; sz<=1; ++sz) { + for(int sy=-1; sy<=1; ++sy) { + for(int sx=-1; sx<=1; ++sx) { + GlobalOrdinal col_id = +get_id(global_nodes_x, global_nodes_y, global_nodes_z, + ix+sx, iy+sy, iz+sz); + + if (col_id >= 0 && col_id < global_nrows) { + ++nnz; + } + } + } + } + row_offset_ptr[r+1] = nnz; + + } + + const MINIFE_GLOBAL_ORDINAL n = row_offsets.size() - 1; + for(int i = 0; i < n; ++i) { + row_offset_ptr[i+1] += row_offset_ptr[i]; + } + + init_matrix(A, rows, row_offsets, row_coords, + global_nodes_x, global_nodes_y, global_nodes_z, global_nrows, mesh); + } + catch(...) { + std::cout << "proc " << myproc << " threw an exception in generate_matrix_structure, probably due to running out of memory." << std::endl; + threw_exc = 1; + } +#ifdef HAVE_MPI + int global_throw = 0; + MPI_Allreduce(&threw_exc, &global_throw, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + threw_exc = global_throw; +#endif + if (threw_exc) { + return 1; + } + + return 0; +} + +}//namespace miniFE + +#endif + diff --git a/tm/src/get_common_files b/tm/src/get_common_files new file mode 100755 index 0000000..e2448e3 --- /dev/null +++ b/tm/src/get_common_files @@ -0,0 +1,15 @@ +#!/bin/bash + +if [ -f ./YAML_Doc.cpp ] ; then + exit 0; +fi + +dir=../../common + +cp ${dir}/YAML_Doc.cpp . +cp ${dir}/YAML_Doc.hpp . +cp ${dir}/YAML_Element.cpp . +cp ${dir}/YAML_Element.hpp . + +cp ${dir}/generate_info_header . + diff --git a/tm/src/main.cpp b/tm/src/main.cpp new file mode 100644 index 0000000..5a38b7f --- /dev/null +++ b/tm/src/main.cpp @@ -0,0 +1,276 @@ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include +#include +#include + +#include + +#include + +#ifdef _OPENMP +#include +#endif + +#ifdef MINIFE_REPORT_RUSAGE +#include +#include +#endif + +#ifdef HAVE_MPI +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +#if MINIFE_INFO != 0 +#include +#else +#include +#endif + +//The following macros should be specified as compile-macros in the +//makefile. They are defaulted here just in case... +#ifndef MINIFE_SCALAR +#define MINIFE_SCALAR double +#endif +#ifndef MINIFE_LOCAL_ORDINAL +#define MINIFE_LOCAL_ORDINAL int +#endif +#ifndef MINIFE_GLOBAL_ORDINAL +#define MINIFE_GLOBAL_ORDINAL int +#endif + +// ************************************************************************ + +void add_params_to_yaml(YAML_Doc& doc, miniFE::Parameters& params); +void add_configuration_to_yaml(YAML_Doc& doc, int numprocs, int numthreads); +void add_timestring_to_yaml(YAML_Doc& doc); + +// +//We will create a 'box' of size nx X ny X nz, partition it among processors, +//then call miniFE::driver which will use the partitioned box as the domain +//from which to assemble finite-element matrices into a global matrix and +//vector, then solve the linear-system using Conjugate Gradients. +// + +int main(int argc, char** argv) { + miniFE::Parameters params; + miniFE::get_parameters(argc, argv, params); + + int numprocs = 1, myproc = 0; + miniFE::initialize_mpi(argc, argv, numprocs, myproc); + + miniFE::timer_type start_time = miniFE::mytimer(); + +#ifdef MINIFE_DEBUG + outstream(numprocs, myproc); +#endif + + if(myproc==0) { + std::cout << "MiniFE Mini-App, OpenMP Peer Implementation" << std::endl; + std::cout << "Creating OpenMP Thread Pool..." << std::endl; + } + int value = 0; + const int thread_count = omp_get_max_threads(); +#pragma omp parallel for reduction(+:value) + for(int i = 0; i < thread_count; ++i) { + value += 1; + } + double global_threadcount; + double local_threadcount = value; + +#ifdef HAVE_MPI + MPI_Allreduce(&local_threadcount,&global_threadcount,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); +#else + global_threadcount = local_threadcount; +#endif + if(myproc==0) { + std::cout << "Counted: " << global_threadcount << " threads." << std::endl; + std::cout << "Running MiniFE Mini-App..." << std::endl; + } + + //make sure each processor has the same parameters: + miniFE::broadcast_parameters(params); + + + Box global_box = { 0, params.nx, 0, params.ny, 0, params.nz }; + std::vector local_boxes(numprocs); + + box_partition(0, numprocs, 2, global_box, &local_boxes[0]); + + Box& my_box = local_boxes[myproc]; + + MINIFE_GLOBAL_ORDINAL num_my_ids = miniFE::get_num_ids(my_box); + MINIFE_GLOBAL_ORDINAL min_ids = num_my_ids; + +#ifdef HAVE_MPI + MPI_Datatype mpi_dtype = miniFE::TypeTraits::mpi_type(); + MPI_Allreduce(&num_my_ids, &min_ids, 1, mpi_dtype, MPI_MIN, MPI_COMM_WORLD); +#endif + + if (min_ids == 0) { + std::cout<<"One or more processors have 0 equations. Not currently supported. Exiting."<. + //To run miniFE with float instead of double, or 'long long' instead of int, + //etc., change these template-parameters by changing the macro definitions in + //the makefile or on the make command-line. + + int return_code = + miniFE::driver< MINIFE_SCALAR, MINIFE_LOCAL_ORDINAL, MINIFE_GLOBAL_ORDINAL>(global_box, my_box, params, doc); + + miniFE::timer_type total_time = miniFE::mytimer() - start_time; + +#ifdef MINIFE_REPORT_RUSAGE + struct rusage get_mem; + getrusage(RUSAGE_SELF, &get_mem); + + long long int rank_rss = get_mem.ru_maxrss; + long long int global_rss = 0; + long long int max_rss = 0; + +#ifdef HAVE_MPI + MPI_Reduce(&rank_rss, &global_rss, 1, + MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD); + MPI_Reduce(&rank_rss, &max_rss, 1, + MPI_LONG_LONG, MPI_MAX, 0, MPI_COMM_WORLD); + if (myproc == 0) { + doc.add("Global All-RSS (kB)", global_rss); + doc.add("Global Max-RSS (kB)", max_rss); + } +#else + doc.add("RSS (kB)", rank_rss); +#endif +#endif + + if (myproc == 0) { + doc.add("Total Program Time",total_time); + doc.generateYAML(); + } + + + miniFE::finalize_mpi(); + + return return_code; +} + +void add_params_to_yaml(YAML_Doc& doc, miniFE::Parameters& params) +{ + doc.add("Global Run Parameters",""); + doc.get("Global Run Parameters")->add("dimensions",""); + doc.get("Global Run Parameters")->get("dimensions")->add("nx",params.nx); + doc.get("Global Run Parameters")->get("dimensions")->add("ny",params.ny); + doc.get("Global Run Parameters")->get("dimensions")->add("nz",params.nz); + doc.get("Global Run Parameters")->add("load_imbalance", params.load_imbalance); + if (params.mv_overlap_comm_comp == 1) { + std::string val("1 (yes)"); + doc.get("Global Run Parameters")->add("mv_overlap_comm_comp", val); + } + else { + std::string val("0 (no)"); + doc.get("Global Run Parameters")->add("mv_overlap_comm_comp", val); + } +#ifdef _OPENMP + doc.get("Global Run Parameters")->add("OpenMP Max Threads:", omp_get_max_threads()); +#endif +} + +void add_configuration_to_yaml(YAML_Doc& doc, int numprocs, int numthreads) +{ + doc.get("Global Run Parameters")->add("number of processors", numprocs); + + doc.add("Platform",""); + doc.get("Platform")->add("hostname",MINIFE_HOSTNAME); + doc.get("Platform")->add("kernel name",MINIFE_KERNEL_NAME); + doc.get("Platform")->add("kernel release",MINIFE_KERNEL_RELEASE); + doc.get("Platform")->add("processor",MINIFE_PROCESSOR); + + doc.add("Build",""); + doc.get("Build")->add("CXX",MINIFE_CXX); +#if MINIFE_INFO != 0 + doc.get("Build")->add("compiler version",MINIFE_CXX_VERSION); +#endif + doc.get("Build")->add("CXXFLAGS",MINIFE_CXXFLAGS); + std::string using_mpi("no"); +#ifdef HAVE_MPI + using_mpi = "yes"; +#endif + doc.get("Build")->add("using MPI",using_mpi); +} + +void add_timestring_to_yaml(YAML_Doc& doc) +{ + std::time_t rawtime; + struct tm * timeinfo; + std::time(&rawtime); + timeinfo = std::localtime(&rawtime); + std::ostringstream osstr; + osstr.fill('0'); + osstr << timeinfo->tm_year+1900 << "-"; + osstr.width(2); osstr << timeinfo->tm_mon+1 << "-"; + osstr.width(2); osstr << timeinfo->tm_mday << ", "; + osstr.width(2); osstr << timeinfo->tm_hour << "-"; + osstr.width(2); osstr << timeinfo->tm_min << "-"; + osstr.width(2); osstr << timeinfo->tm_sec; + std::string timestring = osstr.str(); + doc.add("Run Date/Time",timestring); +} + diff --git a/tm/src/make_local_matrix.hpp b/tm/src/make_local_matrix.hpp new file mode 100644 index 0000000..76970d0 --- /dev/null +++ b/tm/src/make_local_matrix.hpp @@ -0,0 +1,447 @@ +#ifndef _make_local_matrix_hpp_ +#define _make_local_matrix_hpp_ +#include + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include + +#include + +#ifdef HAVE_MPI +#include +#endif + +namespace miniFE { + +template +void +make_local_matrix(MatrixType& A) +{ +#ifdef HAVE_MPI + int numprocs = 1, myproc = 0; + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myproc); + + if (numprocs < 2) { + A.num_cols = A.rows.size(); + A.has_local_indices = true; + return; + } + + typedef typename MatrixType::GlobalOrdinalType GlobalOrdinal; + typedef typename MatrixType::LocalOrdinalType LocalOrdinal; + typedef typename MatrixType::ScalarType Scalar; + + std::map externals; + LocalOrdinal num_external = 0; + + //Extract Matrix pieces + + size_t local_nrow = A.rows.size(); + GlobalOrdinal start_row = local_nrow>0 ? A.rows[0] : -1; + GlobalOrdinal stop_row = local_nrow>0 ? A.rows[local_nrow-1] : -1; + + // We need to convert the index values for the rows on this processor + // to a local index space. We need to: + // - Determine if each index reaches to a local value or external value + // - If local, subtract start_row from index value to get local index + // - If external, find out if it is already accounted for. + // - If so, then do nothing, + // - otherwise + // - add it to the list of external indices, + // - find out which processor owns the value. + // - Set up communication for sparse MV operation + + /////////////////////////////////////////// + // Scan the indices and transform to local + /////////////////////////////////////////// + + std::vector& external_index = A.external_index; + + for(size_t i=0; i tmp_buffer(numprocs, 0); // Temp buffer space needed below + + // Build list of global index offset + + std::vector global_index_offsets(numprocs, 0); + + tmp_buffer[myproc] = start_row; // This is my start row + + // This call sends the start_row of each ith processor to the ith + // entry of global_index_offsets on all processors. + // Thus, each processor knows the range of indices owned by all + // other processors. + // Note: There might be a better algorithm for doing this, but this + // will work... + + MPI_Datatype mpi_dtype = TypeTraits::mpi_type(); + MPI_Allreduce(&tmp_buffer[0], &global_index_offsets[0], numprocs, mpi_dtype, + MPI_SUM, MPI_COMM_WORLD); + + // Go through list of externals and find the processor that owns each + std::vector external_processor(num_external); + + for(LocalOrdinal i=0; i=0; --j) { + if (global_index_offsets[j] <= cur_ind && global_index_offsets[j] >= 0) { + external_processor[i] = j; + break; + } + } + } + + ///////////////////////////////////////////////////////////////////////// + // Sift through the external elements. For each newly encountered external + // point assign it the next index in the sequence. Then look for other + // external elements who are updated by the same node and assign them the next + // set of index numbers in the sequence (ie. elements updated by the same node + // have consecutive indices). + ///////////////////////////////////////////////////////////////////////// + + size_t count = local_nrow; + std::vector& external_local_index = A.external_local_index; + external_local_index.assign(num_external, -1); + + for(LocalOrdinal i=0; i new_external_processor(num_external, 0); + + for(int i=0; i No external elements are updated by + // processor i. + // tmp_neighbors[i] = x ==> (x-1)/numprocs elements are updated from + // processor i. + /// + //////////////////////////////////////////////////////////////////////// + + std::vector tmp_neighbors(numprocs, 0); + + int num_recv_neighbors = 0; + int length = 1; + + for(LocalOrdinal i=0; i recv_list; + recv_list.push_back(new_external_processor[0]); + for(LocalOrdinal i=1; i send_list(num_send_neighbors, 0); + + // + // first post receives, these are immediate receives + // Do not wait for result to come, will do that at the + // wait call below. + // + int MPI_MY_TAG = 99; + + std::vector request(num_send_neighbors); + for(int i=0; i new_external(num_external); + for(LocalOrdinal i=0; i lengths(num_recv_neighbors); + + ++MPI_MY_TAG; + + // First post receives + + for(int i=0; i& neighbors = A.neighbors; + std::vector& recv_length = A.recv_length; + std::vector& send_length = A.send_length; + + neighbors.resize(num_recv_neighbors, 0); + A.request.resize(num_recv_neighbors); + recv_length.resize(num_recv_neighbors, 0); + send_length.resize(num_recv_neighbors, 0); + + LocalOrdinal j = 0; + for(int i=0; i= A.rows.size()) { +//std::cout<<"start_row: "< +#include +#include +#include +#include +#include + +#ifdef _OPENMP +#include +#endif + +//#include "advisor-annotate.h" + +namespace miniFE { + +template +void +perform_element_loop(const simple_mesh_description& mesh, + const Box& local_elem_box, + MatrixType& A, VectorType& b, + Parameters& /*params*/) +{ + typedef typename MatrixType::ScalarType Scalar; + + int global_elems_x = mesh.global_box[0][1]; + int global_elems_y = mesh.global_box[1][1]; + int global_elems_z = mesh.global_box[2][1]; + + //We will iterate the local-element-box (local portion of the mesh), and + //get element-IDs in preparation for later assembling the FE operators + //into the global sparse linear-system. + + GlobalOrdinal num_elems = get_num_ids(local_elem_box); + std::vector elemIDs(num_elems); + + BoxIterator iter = BoxIterator::begin(local_elem_box); + BoxIterator end = BoxIterator::end(local_elem_box); + + for(size_t i=0; iter != end; ++iter, ++i) { + elemIDs[i] = get_id(global_elems_x, global_elems_y, global_elems_z, + iter.x, iter.y, iter.z); + } + + timer_type t_gn = 0, t_ce = 0, t_si = 0; + timer_type t0 = 0; + + const MINIFE_GLOBAL_ORDINAL elemID_size = elemIDs.size(); + + #pragma omp parallel for shared (elemIDs) + for(MINIFE_GLOBAL_ORDINAL i=0; i < elemID_size; ++i) { + ElemData elem_data; + compute_gradient_values(elem_data.grad_vals); + + get_elem_nodes_and_coords(mesh, elemIDs[i], elem_data); + compute_element_matrix_and_vector(elem_data); + sum_into_global_linear_system(elem_data, A, b); + } + +} + +}//namespace miniFE + +#endif + diff --git a/tm/src/simple_mesh_description.hpp b/tm/src/simple_mesh_description.hpp new file mode 100644 index 0000000..b34d44a --- /dev/null +++ b/tm/src/simple_mesh_description.hpp @@ -0,0 +1,248 @@ + +#ifndef _simple_mesh_description_hpp_ +#define _simple_mesh_description_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include +#include +#include + +namespace miniFE { + +template +class simple_mesh_description { +public: + simple_mesh_description(const Box& global_box_in, const Box& local_box_in) + { + Box local_node_box; + for(int i=0; i<3; ++i) { + global_box[i][0] = global_box_in[i][0]; + global_box[i][1] = global_box_in[i][1]; + local_box[i][0] = local_box_in[i][0]; + local_box[i][1] = local_box_in[i][1]; + local_node_box[i][0] = local_box_in[i][0]; + local_node_box[i][1] = local_box_in[i][1]; + //num-owned-nodes == num-elems+1 in this dimension if the elem box is not empty + //and we are at the high end of the global range in that dimension: + if (local_box_in[i][1] > local_box_in[i][0] && local_box_in[i][1] == global_box[i][1]) local_node_box[i][1] += 1; + } + + int max_node_x = global_box[0][1]+1; + int max_node_y = global_box[1][1]+1; + int max_node_z = global_box[2][1]+1; + create_map_id_to_row(max_node_x, max_node_y, max_node_z, local_node_box, + map_ids_to_rows); + + //As described in analytic_soln.hpp, + //we will impose a 0 boundary-condition on faces x=0, y=0, z=0, y=1, z=1 + //we will impose a 1 boundary-condition on face x=1 + +#ifdef MINIFE_DEBUG +std::cout< 0) --miny; + if (local_node_box[Z][0] > 0) --minz; + if (local_node_box[Y][1] < max_node_y) ++maxy; + if (local_node_box[Z][1] < max_node_z) ++maxz; + + for(int iz=minz; iz(max_node_x, max_node_y, max_node_z, + 0, iy, iz); +#ifdef MINIFE_DEBUG +std::cout<<"x=0 BC, node "< 0) --minx; + if (local_node_box[Z][0] > 0) --minz; + if (local_node_box[X][1] < max_node_x) ++maxx; + if (local_node_box[Z][1] < max_node_z) ++maxz; + + for(int iz=minz; iz(max_node_x, max_node_y, max_node_z, + ix, 0, iz); +#ifdef MINIFE_DEBUG +std::cout<<"y=0 BC, node "< 0) --minz; + if (local_node_box[Y][0] > 0) --miny; + if (local_node_box[Z][1] < max_node_z) ++maxz; + if (local_node_box[Y][1] < max_node_y) ++maxy; + + for(int iy=miny; iy(max_node_x, max_node_y, max_node_z, + x1, iy, iz); + GlobalOrdinal row = map_id_to_row(nodeID); +#ifdef MINIFE_DEBUG +std::cout<<"x=1 BC, node "< 0) --minz; + if (local_node_box[X][0] > 0) --minx; + if (local_node_box[Z][1] < max_node_z) ++maxz; + if (local_node_box[X][1] < max_node_x) ++maxx; + + for(int ix=minx; ix(max_node_x, max_node_y, max_node_z, + ix, y1, iz); +#ifdef MINIFE_DEBUG +std::cout<<"y=1 BC, node "< 0) --miny; + if (local_node_box[X][0] > 0) --minx; + if (local_node_box[Y][1] < max_node_y) ++maxy; + if (local_node_box[X][1] < max_node_x) ++maxx; + + for(int ix=minx; ix(max_node_x, max_node_y, max_node_z, + ix, iy, z1); +#ifdef MINIFE_DEBUG +std::cout<<"z=1 BC, node "<::const_iterator mend = map_ids_to_rows.end(); + --mend; + return mend->second; + } + std::set bc_rows_0; + std::set bc_rows_1; + std::map map_ids_to_rows; + Box global_box; + Box local_box; +};//class simple_mesh_description + +}//namespace miniFE + +#endif diff --git a/tm/src/time_kernels.hpp b/tm/src/time_kernels.hpp new file mode 100644 index 0000000..e14ff09 --- /dev/null +++ b/tm/src/time_kernels.hpp @@ -0,0 +1,139 @@ +#ifndef _time_kernels_hpp_ +#define _time_kernels_hpp_ + +//@HEADER +// ************************************************************************ +// +// MiniFE: Simple Finite Element Assembly and Solve +// Copyright (2006-2013) Sandia Corporation +// +// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive +// license for use of this work by or on behalf of the U.S. Government. +// +// This library 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 2.1 of the +// License, or (at your option) any later version. +// +// This library 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 this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +// USA +// +// ************************************************************************ +//@HEADER + +#include + +#include +#include + +#ifdef MINIFE_HAVE_CUDA +#include +#endif + +namespace miniFE { + +template +void +time_kernels(OperatorType& A, + const VectorType& b, + VectorType& x, + Matvec matvec, + typename OperatorType::LocalOrdinalType max_iter, + typename OperatorType::ScalarType& xdotp, + timer_type* my_kern_times) +{ + typedef typename OperatorType::ScalarType ScalarType; + typedef typename OperatorType::LocalOrdinalType OrdinalType; + typedef typename TypeTraits::magnitude_type magnitude_type; + + timer_type t0 = 0, tWAXPY = 0, tDOT = 0, tMATVEC = 0; + + int myproc = 0; +#ifdef HAVE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &myproc); +#endif + + if (!A.has_local_indices) { + std::cerr << "miniFE::time_kernels ERROR, A.has_local_indices is false, needs to be true. This probably means " + << "miniFE::make_local_matrix(A) was not called prior to calling miniFE::time_kernels." + << std::endl; + return; + } + + OrdinalType nrows = A.rows.size(); + OrdinalType ncols = A.num_cols; + + VectorType p(0, ncols, b.compute_node); + + ScalarType one = 1.0; + ScalarType zero = 0.0; + + typedef typename VectorType::ComputeNodeType ComputeNodeType; + ComputeNodeType& compute_node = x.compute_node; + + //The following lines that create and initialize buffers are no-ops in many + //cases, but perform actual allocations and copies if a off-cpu device such as + //a GPU is being used by compute_node. + + //Do any required allocations for buffers that will be needed during CG: + ScalarType* d_x = compute_node.get_buffer(&x.coefs[0], x.coefs.size()); + ScalarType* d_p = compute_node.get_buffer(&p.coefs[0], p.coefs.size()); + ScalarType* d_b = compute_node.get_buffer(&b.coefs[0], b.coefs.size()); + OrdinalType* d_Arowoff = compute_node.get_buffer(&A.row_offsets[0], A.row_offsets.size()); + OrdinalType* d_Acols = compute_node.get_buffer(&A.packed_cols[0], A.packed_cols.size()); + ScalarType* d_Acoefs = compute_node.get_buffer(&A.packed_coefs[0], A.packed_coefs.size()); + + //Copy data to buffers that need to be initialized from input data: + compute_node.copy_to_buffer(&x.coefs[0], x.coefs.size(), d_x); + compute_node.copy_to_buffer(&b.coefs[0], b.coefs.size(), d_b); + compute_node.copy_to_buffer(&A.row_offsets[0], A.row_offsets.size(), d_Arowoff); + compute_node.copy_to_buffer(&A.packed_cols[0], A.packed_cols.size(), d_Acols); + compute_node.copy_to_buffer(&A.packed_coefs[0], A.packed_coefs.size(), d_Acoefs); + + TICK(); + for(OrdinalType i=0; i