MathOperators_vector.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 vectors.                                                  //
// ================================================================================== //
// INFO                                                                               //
// ================================================================================== //
// Author     : Alessandro Alaia                                                      //
// Version    : v3.0                                                                  //
//                                                                                    //
// All rights reserved.                                                               //
// ================================================================================== //
 
// ================================================================================== //
// TEMPLATE IMPLEMENTATIONS                                                           //
// ================================================================================== //
// Operator "min" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< T > min(
    const std::vector< T >                      &x,
    const std::vector< T >                      &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
int            n = std::min(x.size(), y.size());
int            m = std::max(x.size(), y.size());
std::vector< T >    z(m);
 
// ================================================================================== //
// COMPARE VECTORS                                                                    //
// ================================================================================== //
for (int i = 0; i < n; i++) {
    z[i] = min(x[i], y[i]);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< T > min(
    const std::vector< T >                      &x,
    const T                                     &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
int            n = x.size();
std::vector< T >    z(n);
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
for (int i = 0; i < n; i++) {
    z[i] = min(x[i], y);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< T > min(
    const T                                     &x,
    const std::vector< T >                      &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::vector< T >    z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
z = min(y, x);
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< std::vector< T > > min(
    const std::vector< std::vector < T > >      &x,
    const T                                     &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
int                      n = x.size();
std::vector< std::vector< T > >    z(n);
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
for (int i = 0; i < n; i++) {
    z[i] = min(x[i], y);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< std::vector < T > > min(
    const T                                     &x,
    const std::vector< std::vector< T > >       &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::vector< std::vector< T > >    z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
z = min(y, x);
 
return(z); };
 
// Operator "minval" ================================================================ //
 
// ---------------------------------------------------------------------------------- //
template <typename T, typename std::enable_if< std::is_scalar< T >::value>::type*>
void inline minval(
    const T                                     &x,
    T                                           &min_value
) {
 
// ================================================================================== //
// template <typename T,                                                              //
//           typename std::enable_if< std::is_scalar< T >::value>::type* = nullptr>   //
// void inline minval(                                                                //
//     const T              &x,                                                       //
//     T                    &min_value)                                               //
//                                                                                    //
// Overloading of minval function for scalar types.                                   //
// ================================================================================== //
// INPUT                                                                              //
// ================================================================================== //
// - x         : T, input scalar                                                      //
// - min_value : T, on output stores the input value.                                 //
// ================================================================================== //
// OUTPUT                                                                             //
// ================================================================================== //
// - none                                                                             //
// ================================================================================== //
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
// none
 
// Counters
// none
 
// ================================================================================== //
// COMPUTE THE MIN VALUE                                                              //
// ================================================================================== //
min_value = x;
 
return; };
 
// ---------------------------------------------------------------------------------- //
template <class T, class T1>
void minval(
    const std::vector<T>                        &x,
    T1                                          &min_value
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
T1               value;
 
// Counters
int              i, n;
 
// ================================================================================== //
// FIND THE MIN-VALUE                                                                 //
// ================================================================================== //
n = x.size();
if (n > 0) {
    minval(x[0], min_value);
    for (i = 1; i < n; i++) {
        minval(x[i], value);
        if (value < min_value) {
            min_value = value;
        }
    } //next i
}
 
return; };
 
// Operator "max" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector<T> max(
    const std::vector< T >                      &x,
    const std::vector< T >                      &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
int            n = std::min(x.size(), y.size());
int            m = std::max(x.size(), y.size());
std::vector< T >    z(m);
 
// ================================================================================== //
// COMPARE VECTORS                                                                    //
// ================================================================================== //
for (int i = 0; i < n; i++) {
    z[i] = max(x[i], y[i]);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< T > max(
    const std::vector< T >                      &x,
    const T                                     &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
int            n = x.size();
std::vector<T>      z(n);
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
for (int i = 0; i < n; i++) {
    z[i] = max(x[i], y);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< T > max(
    const T                                     &x,
    const std::vector< T >                      &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::vector< T >       z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
z = max(y, x);
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< std::vector< T > > max(
    const std::vector< std::vector< T > >       &x,
    const T                                     &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
int                        n = x.size();
std::vector< std::vector< T > >      z(n);
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
for (int i = 0; i < n; i++) {
    z[i] = max(x[i], y);
} //next i
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< std::vector< T > > max(
    const T                                     &x,
    const std::vector< std::vector< T > >       &y
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
std::vector< std::vector< T > >       z;
 
// ================================================================================== //
// COMPARE VECTOR AND SCALAR                                                          //
// ================================================================================== //
z = max(y, x);
 
return(z); };
 
// Operator "maxval" ================================================================ //
 
// ---------------------------------------------------------------------------------- //
template <typename T, typename std::enable_if< std::is_scalar< T >::value>::type*>
void inline maxval(
    const T                                     &x,
    T                                           &max_value
) {
 
// ================================================================================== //
// template <typename T,                                                              //
//           typename std::enable_if< std::is_scalar< T >::value>::type* = nullptr>   //
// void inline maxval(                                                                //
//     const T              &x,                                                       //
//     T                    &max_value)                                               //
//                                                                                    //
// Overloading of maxval function for scalar types.                                   //
// ================================================================================== //
// INPUT                                                                              //
// ================================================================================== //
// - x         : T, input scalar                                                      //
// - max_value : T, on output stores the input value                                  //
// ================================================================================== //
// OUTPUT                                                                             //
// ================================================================================== //
// - none                                                                             //
// ================================================================================== //
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
// none
 
// Counters
// none
 
// ================================================================================== //
// COMPUTE THE MIN VALUE                                                              //
// ================================================================================== //
max_value = x;
 
return; };
 
// ---------------------------------------------------------------------------------- //
template <class T, class T1>
void maxval(
    const std::vector<T>                        &x,
    T1                                          &max_value
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
T1               value;
 
// Counters
int              i, n;
 
// ================================================================================== //
// FIND THE MIN-VALUE                                                                 //
// ================================================================================== //
n = x.size();
if (n > 0) {
    maxval(x[0], max_value);
    for (i = 1; i < n; i++) {
        maxval(x[i], value);
        if (value > max_value) {
            max_value = value;
        }
    } //next i
}
 
return; };
 
// Operator "sum" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template<class T, typename std::enable_if< std::is_scalar< T >::value>::type*>
void inline sum(
    const T                                     &x,
    T                                           &s
) {
 
// ================================================================================== //
// template<class T,                                                                  //
//          typename std::enable_if< std::is_scalar< T >::value>::type* = nullptr>    //
// void inline sum(                                                                   //
//     const T              &x,                                                       //
//     T                    &s)                                                       //
//                                                                                    //
// Overloading of sum() function for scalar type.                                     //
// ================================================================================== //
// INPUT                                                                              //
// ================================================================================== //
// - x         : int, dummy input                                                     //
// - s         : int, dummy output                                                    //
// ================================================================================== //
// OUTPUT                                                                             //
// ================================================================================== //
// - none                                                                             //
// ================================================================================== //
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
// none
 
// Counters
// none
 
// ================================================================================== //
// PERFORM DUMMY SUM                                                                  //
// ================================================================================== //
s = x;
 
return; };
 
// ---------------------------------------------------------------------------------- //
template <class T, class T1>
void sum(
    const std::vector< T >                      &x,
    T1                                          &s
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
int           n = x.size();
T1            value;
 
// Counters
int           i;
 
// ================================================================================== //
// PERFORM SUMMATION                                                                  //
// ================================================================================== //
if (n > 0) {
    sum(x[0], s);
    for (i = 1; i < n; i++) {
        sum(x[i], value);
        s += value;
    } //next i
}
 
return; };
 
// Operator "abs" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector<T> abs(
    const std::vector< T >                      &x
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
int                     n = x.size();
std::vector< T >             z;
 
// Counters
int                     i;
 
// ================================================================================== //
// COMPUTE THE ABSOLUTE VALUE OF A VECTOR                                             //
// ================================================================================== //
if (n > 0) {
    z.resize(n);
    for (i = 0; i < n; i++) {
        z[i] = abs(x[i]);
    } //next i
}
 
return(z); };
 
// Operator "pow" =================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector< T > pow(
    std::vector< T >                            &x,
    double                                       p
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
int              n;
std::vector< T >      y;
 
// Counters
int              i;
 
// ================================================================================== //
// COMPUTE ELEMENT-WISE POWER                                                         //
// ================================================================================== //
n = x.size();
y.resize(n);
for (i = 0; i < n; i++) {
    y[i] = static_cast<T>(pow(x[i], p));
} //next i
 
return(y); };
 
// Operator "norm" ================================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T>
double norm1(
    const std::vector< T >                      &x
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
int             n = x.size();
double          z = 0.0;
 
// Counters
int             i;
 
// ================================================================================== //
// COMPUTE THE P-NORM                                                                 //
// ================================================================================== //
if (n > 0) {
    for (i = 0; i < n; i++) {
        z += abs(x[i]);
    } //next i
}
 
return(z); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
double norm2(
    const std::vector< T >                      &x
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
int             n = x.size();
double          z = 0.0;
 
// Counters
int             i;
 
// ================================================================================== //
// COMPUTE THE P-NORM                                                                 //
// ================================================================================== //
if (n > 0) {
    for (i = 0; i < n; i++) {
        z += x[i]*x[i];
    } //next i
}
 
return(sqrt(z)); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
double norm(
    const std::vector< T >                      &x,
    int                                          p
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
int             n = x.size();
double          z = 0.0;
double          t, y;
 
// Counters
int             i, j;
 
// ================================================================================== //
// COMPUTE THE P-NORM                                                                 //
// ================================================================================== //
if (p == 1) { return(norm1(x)); }
if (p == 2) { return(norm2(x)); }
 
if (n > 0) {
    for (i = 0; i < n; i++) {
        y = 1.0;
        t = x[i];
        for (j = 1; j <= p; j++) {
            y = y*t;
        } //next j
        z += abs(y);
    } //next i
}
 
return(std::exp(std::log(std::max(z, 1.0e-307))/((double) p))); };
 
// ---------------------------------------------------------------------------------- //
template <class T>
double normInf(
    const std::vector< T >                      &x
) {
 
using namespace std;
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
int             n = x.size();
double          z = 0.0, y;
 
// Counters
int             i;
 
// ================================================================================== //
// COMPUTE THE inf-NORM                                                               //
// ================================================================================== //
if (n > 0) {
    z = abs(x[0]);
    for (i = 1; i < n; i++) {
        y = abs(x[i]);
        if (y > z) {
            z = y;
        }
    } //next i
}
 
return(z); };
 
// Operator "dotProduct" =========================================================== //
 
// ---------------------------------------------------------------------------------- //
template <class T>
T dotProduct(
    const std::vector< T >                      &x,
    const std::vector< T >                      &y
) {
 
// ================================================================================== //
// VARIABLES DECLARATION                                                              //
// ================================================================================== //
 
// Local variables
int                 n = x.size(), m = y.size();
T                   dp = ((T) 0.0);
 
// Counters
int                 i;
 
// ================================================================================== //
// COMPUTE THE DOT PRODUCT                                                            //
// ================================================================================== //
if ((n > 0) && (n == m)) {
    for (i = 0; i < n; i++) {
        dp += x[i]*y[i];
    } //next i
}
 
return(dp); };
 
// Operator "crossProduct" ========================================================= //
 
// ---------------------------------------------------------------------------------- //
template <class T>
std::vector<T> crossProduct(
    const std::vector<T>                        &x,
    const std::vector<T>                        &y
) {
 
// =================================================================================== //
// VARIABLES DECLARATION                                                               //
// =================================================================================== //
std::vector<T>      z(3, T{});
 
// =================================================================================== //
// 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);}