#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