| #ifndef _utils_hpp_ |
| #define _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 <cstdlib> |
| #include <cmath> |
| #include <vector> |
| #include <map> |
| |
| #ifdef HAVE_MPI |
| #include <mpi.h> |
| #endif |
| |
| #include <TypeTraits.hpp> |
| #include <Parameters.hpp> |
| |
| namespace miniFE { |
| |
| void get_parameters(int argc, char** argv, Parameters& params); |
| |
| void broadcast_parameters(Parameters& params); |
| |
| void initialize_mpi(int argc, char** argv, int& numprocs, int& myproc); |
| |
| void finalize_mpi(); |
| |
| template<typename Scalar> |
| Scalar percentage_difference(Scalar value, Scalar average) |
| { |
| //result will be the difference between value and average, represented as |
| //a percentage of average. |
| //Examples: |
| // if value=100 and average=50, result is 100% |
| // if value=500 and average=400, result is 25% |
| |
| //Note: if average is 0, result is undefined. We'll return -1.0; |
| |
| Scalar result = std::abs(value-average); |
| if (std::abs(average) > 1.e-5) { |
| result /= average; |
| result *= 100; |
| } |
| else result = -1; |
| |
| return result; |
| } |
| |
| template<typename GlobalOrdinal> |
| void get_global_min_max(GlobalOrdinal local_n, |
| GlobalOrdinal& global_n, |
| GlobalOrdinal& min_n, |
| int& min_proc, |
| GlobalOrdinal& max_n, |
| int& max_proc) |
| { |
| //Given a local_n, compute global_n, min/max, etc. All computed results |
| //will be returned on all processors. |
| // |
| int numprocs = 1, myproc = 0; |
| #ifdef HAVE_MPI |
| MPI_Comm_size(MPI_COMM_WORLD, &numprocs); |
| MPI_Comm_rank(MPI_COMM_WORLD, &myproc); |
| #endif |
| |
| std::vector<GlobalOrdinal> all_n(numprocs, 0); |
| all_n[myproc] = local_n; |
| #ifdef HAVE_MPI |
| std::vector<GlobalOrdinal> tmp(all_n); |
| MPI_Datatype mpi_dtype = TypeTraits<GlobalOrdinal>::mpi_type(); |
| MPI_Allreduce(&tmp[0], &all_n[0], numprocs, mpi_dtype, MPI_MAX, MPI_COMM_WORLD); |
| #endif |
| |
| global_n = 0; |
| min_n= 5*local_n; |
| min_proc = 0; |
| max_n= 0; |
| max_proc = 0; |
| |
| for(int i=0; i<numprocs; ++i) { |
| global_n += all_n[i]; |
| //min_proc will be the lowest-numbered proc with n = min_n |
| if (all_n[i] < min_n) { |
| min_n = all_n[i]; |
| min_proc = i; |
| } |
| //max_proc will be the highest-numbered proc with n = max_n |
| if (all_n[i] >= max_n) { |
| max_n = all_n[i]; |
| max_proc = i; |
| } |
| } |
| } |
| |
| template<typename Scalar> |
| Scalar compute_std_dev_as_percentage(Scalar local_nrows, |
| Scalar avg_nrows) |
| { |
| //compute and return a standard deviation for the deviation of local_nrows from the average. |
| //the std. dev. will be expressed as a percentage of avg_nrows. |
| // |
| //Input argument local_nrows is really a integer, but taking it as a floating-point scalar is |
| //harmless. |
| // |
| #ifdef HAVE_MPI |
| int numprocs = 1, myproc = 0; |
| MPI_Comm_size(MPI_COMM_WORLD, &numprocs); |
| MPI_Comm_rank(MPI_COMM_WORLD, &myproc); |
| MPI_Datatype mpi_dtype = TypeTraits<Scalar>::mpi_type(); |
| |
| //If it's significantly more efficient, we may consider using MPI_Gather below instead of |
| //MPI_Allgather. We really only need to compute std.dev. on proc 0... |
| // |
| //(But for now, use MPI_Allgather and compute on all procs.) |
| |
| std::vector<Scalar> all_nrows(numprocs, 0); |
| MPI_Allgather(&local_nrows, 1, mpi_dtype, &all_nrows[0], 1, mpi_dtype, MPI_COMM_WORLD); |
| |
| //turn all_nrows contents into deviations, add to sum-of-squares-of-deviations: |
| Scalar sum_sqr_dev = 0; |
| for(size_t i=0; i<all_nrows.size(); ++i) { |
| all_nrows[i] -= avg_nrows; |
| all_nrows[i] *= all_nrows[i]; |
| sum_sqr_dev += all_nrows[i]; |
| } |
| Scalar tmp1 = sum_sqr_dev; |
| Scalar std_dev = numprocs>1 ? std::sqrt(tmp1/(numprocs-1)) : 0; |
| |
| //std_dev is now the standard deviation of rows-per-processor with respect |
| //to avg_nrows. |
| //Next turn std_dev into a percentage of avg_nrows: |
| std_dev /= avg_nrows; |
| std_dev *= 100; |
| return std_dev; |
| #else |
| return 0; |
| #endif |
| } |
| |
| template<typename GlobalOrdinal> |
| GlobalOrdinal find_row_for_id(GlobalOrdinal id, |
| const std::map<GlobalOrdinal,GlobalOrdinal>& ids_to_rows) |
| { |
| typename std::map<GlobalOrdinal,GlobalOrdinal>::const_iterator |
| iter = ids_to_rows.lower_bound(id); |
| |
| if (iter == ids_to_rows.end() || iter->first != id) { |
| if (ids_to_rows.size() > 0) { |
| --iter; |
| } |
| else { |
| std::cout << "ERROR, failed to map id to row."<<std::endl; |
| return -99; |
| } |
| } |
| |
| if (iter->first == id) { |
| return iter->second; |
| } |
| |
| if (iter == ids_to_rows.begin() && iter->first > id) { |
| std::cout << "ERROR, id:" << id << ", ids_to_rows.begin(): " << iter->first<<std::endl; |
| return -99; |
| } |
| |
| GlobalOrdinal offset = id - iter->first; |
| |
| if (offset < 0) { |
| std::cout << "ERROR, negative offset in find_row_for_id for id="<<id<<std::endl; |
| return -99; |
| } |
| |
| return iter->second + offset; |
| } |
| |
| }//namespace miniFE |
| |
| #endif |
| |