MathOperators_array.tpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit 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 bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
// ================================================================================== //
//                                 MATH OPERATORS                                     //
//                                                                                    //
// Basic math operators for std::array.                                               //
// ================================================================================== //
// INFO                                                                               //
// ================================================================================== //
// Author     : Alessandro Alaia                                                      //
// Version    : v3.0                                                                  //
//                                                                                    //
// All rights reserved.                                                               //
// ================================================================================== //
 
// ================================================================================== //
// TEMPLATE IMPLEMENTATIONS                                                           //
// ================================================================================== //
 
// Operator "min" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
std::array<T, d> min(
    const std::array<T, d>                      &x,
    const std::array<T, d>                      &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<T, d>          z;
 
// ================================================================================== //
// COMPARE VECTORS                                                                    //
// ================================================================================== //
for (size_t i = 0; i < d; i++) {
    z[i] = min(x[i], y[i]);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
std::array<T, d> min(
    const std::array<T, d>                      &x,
    const T                                     &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<T, d>    z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
for (size_t i = 0; i < d; i++) {
    z[i] = min(x[i], y);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
std::array<T, d> min(
    const T                                     &x,
    const std::array<T, d>                      &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<T, d>    z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
z = min(y, x);
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d, size_t n>
std::array<std::array<T, n>, d> min(
    const std::array<std::array<T, n>, d>       &x,
    const T                                     &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<std::array<T, n>, d>    z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
for (size_t i = 0; i < d; i++) {
    z[i] = min(x[i], y);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d, size_t n>
std::array<std::array<T, n>, d> min(
    const T                                     &x,
    const std::array<std::array<T, n>, d>       &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<std::array<T, n>, d>    z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
z = min(y, x);
 
return(z); };
 
// Operator "minval" ================================================================ //
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d, class T1>
void minval(
    const std::array<T, d>                      &x,
    T1                                          &min_value
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
T1               value;
 
// Counters
size_t           i;
 
// ================================================================================== //
// FIND THE MIN-VALUE                                                                 //
// ================================================================================== //
if (d > 0) {
    minval(x[0], min_value);
    for (i = 1; i < d; i++) {
        minval(x[i], value);
        if (value < min_value) {
            min_value = value;
        }
    } //next i
}
 
return; };
 
// Operator "max" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
std::array<T, d> max(
    const std::array<T, d>                      &x,
    const std::array<T, d>                      &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<T, d>    z;
 
// ================================================================================== //
// COMPARE VECTORS                                                                    //
// ================================================================================== //
for (size_t i = 0; i < d; i++) {
    z[i] = max(x[i], y[i]);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
std::array<T, d> max(
    const std::array<T, d>                      &x,
    const T                                     &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<T, d>      z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
for (size_t i = 0; i < d; i++) {
    z[i] = max(x[i], y);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
std::array<T, d> max(
    const T                                     &x,
    const std::array<T, d>                      &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<T, d>       z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
z = max(y, x);
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d, size_t n>
std::array<std::array<T, n>, d> max(
    const std::array<std::array<T, n>, d>       &x,
    const T                                     &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<std::array<T, n>, d>      z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
for (size_t i = 0; i < d; i++) {
    z[i] = max(x[i], y);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d, size_t n>
std::array<std::array<T, n>, d> max(
    const T                                     &x,
    const std::array<std::array<T, n>, d>       &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::array<std::array<T, n>, d>       z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
z = max(y, x);
 
return(z); };
 
 
// Operator "maxval" ================================================================ //
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d, class T1>
void maxval(
    const std::array<T, d>                      &x,
    T1                                          &max_value
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
T1               value;
 
// Counters
size_t           i;
 
// ================================================================================== //
// FIND THE MIN-VALUE                                                                 //
// ================================================================================== //
if (d > 0) {
    maxval(x[0], max_value);
    for (i = 1; i < d; i++) {
        maxval(x[i], value);
        if (value > max_value) {
            max_value = value;
        }
    } //next i
}
 
return; };
 
// Operator "sum" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d, class T1>
void sum(
    const std::array<T, d>                      &x,
    T1                                          &s
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
T1            value;
 
// Counters
size_t        i;
 
// ================================================================================== //
// PERFORM SUMMATION                                                                  //
// ================================================================================== //
if (d > 0) {
    sum(x[0], s);
    for (i = 1; i < d; i++) {
        sum(x[i], value);
        s += value;
    } //next i
}
 
return; };
 
// Operator "abs" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
std::array<T, d> abs(
    const std::array<T, d>                      &x
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
std::array<T, d>             z;
 
// Counters
size_t                       i;
 
// ================================================================================== //
// COMPUTE THE ABSOLUTE VALUE OF A VECTOR                                             //
// ================================================================================== //
if (d > 0) {
    for (i = 0; i < d; ++i) {
        z[i] = abs(x[i]);
    } //next i
}
 
return(z); };
 
// Operator "pow" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
std::array<T, d> pow(
    std::array<T, d>                            &x,
    double                                       p
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
std::array<T, d>      y;
 
// Counters
size_t                i;
 
// ================================================================================== //
// COMPUTE ELEMENT-WISE POWER                                                         //
// ================================================================================== //
for (i = 0; i < d; i++) {
    y[i] = static_cast<T>(pow(x[i], p));
} //next i
 
return(y); };
 
// Operator "norm" ================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
double norm1(
    const std::array<T, d>                      &x
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
double          z = 0.0;
 
// Counters
size_t          i;
 
// ================================================================================== //
// COMPUTE THE P-NORM                                                                 //
// ================================================================================== //
if (d > 0) {
    for (i = 0; i < d; i++) {
        z += abs(x[i]);
    } //next i
}
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
double norm2(
    const std::array<T, d>                      &x
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
double          z = 0.0;
 
// Counters
size_t          i;
 
// ================================================================================== //
// COMPUTE THE P-NORM                                                                 //
// ================================================================================== //
if (d > 0) {
    for (i = 0; i < d; i++) {
        z += x[i]*x[i];
    } //next i
}
 
return(sqrt(z)); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
double norm(
    const std::array<T, d>                      &x,
    int                                          p
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
double          z = 0.0;
double          t, y;
 
// Counters
size_t          i;
int             j;
 
// ================================================================================== //
// COMPUTE THE P-NORM                                                                 //
// ================================================================================== //
if (p == 1) { return(norm1(x)); }
if (p == 2) { return(norm2(x)); }
 
if (d > 0) {
    for (i = 0; i < d; i++) {
        y = 1.0;
        t = abs(x[i]);
        for (j = 1; j <= p; j++) {
            y = y*t;
        } //next j
        //z += pow(abs(x[i]), p);
        z += y;
    } //next i
}
 
return(std::exp(std::log(std::max(z, 1.0e-307))/((double) p))); };
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
double normInf(
    const std::array<T, d>                      &x
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
double          z = 0.0, y;
 
// Counters
size_t          i;
 
// ================================================================================== //
// COMPUTE THE inf-NORM                                                               //
// ================================================================================== //
if (d > 0) {
    z = abs(x[0]);
    for (i = 1; i < d; i++) {
        y = abs(x[i]);
        if (y > z) {
            z = y;
        }
    } //next i
}
 
return(z); };
 
// Operator "dotProduct" =========================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T, size_t d>
T dotProduct(
    const std::array<T, d>                      &x,
    const std::array<T, d>                      &y
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
T                   dp = ((T) 0.0);
 
// Counters
size_t              i;
 
// ================================================================================== //
// COMPUTE THE DOT PRODUCT                                                            //
// ================================================================================== //
if (d > 0) {
    for (i = 0; i < d; i++) {
        dp += x[i]*y[i];
    } //next i
}
 
return(dp); };
 
// Operator "crossProduct" ========================================================= //
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::array<T, 3> crossProduct(
    const std::array<T, 3>                      &x,
    const std::array<T, 3>                      &y
) {
 
// =================================================================================== //
// VARIABLES DECLARATION                                                               //
// =================================================================================== //
std::array<T, 3>      z;
 
// =================================================================================== //
// COMPUTE THE EXTERNAL PRODUCT                                                        //
// =================================================================================== //
z[0] = x[1] * y[2] - x[2] * y[1];
z[1] = x[2] * y[0] - x[0] * y[2];
z[2] = x[0] * y[1] - x[1] * y[0];
 
return (z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
T crossProduct(
    const std::array<T, 2>                      &x,
    const std::array<T, 2>                      &y
) {
 
// =================================================================================== //
// VARIABLES DECLARATION                                                               //
// =================================================================================== //
T      z;
 
// =================================================================================== //
// COMPUTE THE EXTERNAL PRODUCT                                                        //
// =================================================================================== //
z = x[0] * y[1] - x[1] * y[0];
 
return (z); };