matrix_utilities.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/>.
 *
\*---------------------------------------------------------------------------*/
 
// ========================================================================== //
//                         LINEAR ALGEBRA PACKAGE                             //
//                                                                            //
// Functions for basic linear algebra computations.                           //
// ========================================================================== //
// INFO                                                                       //
// ========================================================================== //
// Author            : Alessandro Alaia                                       //
// Data              : Sept 26, 2014                                          //
// Version           : v2.0                                                   //
//                                                                            //
// All rights reserved.                                                       //
// ========================================================================== //
 
namespace bitpit{
 
namespace linearalgebra{
// -------------------------------------------------------------------------- //
template<class T>
void display(
    std::ostream                                &out,
    std::vector< std::vector< T > >             &A
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
int             m = A.size();
 
// Counters
int             i;
 
// ========================================================================== //
// DISPLAY MATRIX CONTENT                                                     //
// ========================================================================== //
for (i = 0; i < m; ++i) {
    out << A[i] << std::endl;
} //next i
 
return; };
 
// -------------------------------------------------------------------------- //
template<class T, size_t m, size_t n>
void display(
    std::ostream                                &out,
    std::array<std::array<T, n>, m>             &A
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
// none
 
// Counters
size_t          i;
 
// ========================================================================== //
// DISPLAY MATRIX CONTENT                                                     //
// ========================================================================== //
for (i = 0; i < m; ++i) {
    out << A[i] << std::endl;
} //next i
 
return; };
 
// -------------------------------------------------------------------------- //
template<class T>
void displayColMajor(
    std::ostream &out,
    T *A,
    int nRows,
    int nCols
) {
    for (int i = 0; i < nRows; ++i) {
        for (int j = 0; j < nCols; ++j) {
            out << A[linearIndexColMajor(i, j, nRows, nCols)] << "  ";
        }
        out << std::endl;
    }
}
 
// -------------------------------------------------------------------------- //
template<class T>
void displayRowMajor(
    std::ostream &out,
    T *A,
    int nRows,
    int nCols
) {
    for (int i = 0; i < nRows; ++i) {
        for (int j = 0; j < nCols; ++j) {
            out << A[linearIndexRowMajor(i, j, nRows, nCols)] << "  ";
        }
        out << std::endl;
    }
}
 
// -------------------------------------------------------------------------- //
template<class T>
void displayColMajorSymmetric(
    std::ostream &out, T *A,
    int nRows,
    int nCols,
    char uplo
) {
    assert(uplo=='L' || uplo=='U');
 
    for (int i = 0; i < nRows; ++i) {
        for (int j = 0; j < nCols; ++j) {
            out << A[linearIndexColMajorSymmetric(i, j, nRows, nCols, uplo)] << "  ";
        }
        out << std::endl;
    }
}
 
// -------------------------------------------------------------------------- //
template<class T>
void displayRowMajorSymmetric(
    std::ostream &out,
    T *A,
    int nRows,
    int nCols,
    char uplo
) {
    assert(uplo=='L' || uplo=='U');
 
    for (int i = 0; i < nRows; ++i) {
        for (int j = 0; j < nCols; ++j) {
            out << A[linearIndexRowMajorSymmetric(i, j, nRows, nCols, uplo)] << "  ";
        }
        out << std::endl;
    }
}
 
// -------------------------------------------------------------------------- //
template <class T>
void complement(
    int                                          i,
    int                                          j,
    std::vector< std::vector< T > >             &A,
    std::vector< std::vector< T > >             &B
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
int            m, n;
 
// Counters
int            l, k;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
m = A.size();
if (m == 0) { return; }
n = A[0].size();
if (n == 0) { return; }
if ((i >= m) || (i < 0)) {
    return;
}
if ((j >= n) || (j < 0)) {
    return;
}
 
// ========================================================================== //
// EXTRACT COMPLEMENT                                                         //
// ========================================================================== //
 
// Resize output variables
B.resize(m-1);
for (int p = 0; p < m-1; ++p) {
    B[p].resize(n-1, T{});
} //next p
 
// Extract complement
for (l = 0; l < i; l++) {
    for (k = 0; k < j; k++) {
        B[l][k] = A[l][k];
    } //next k
    for (k = j+1; k < n; k++) {
        B[l][k-1] = A[l][k];
    } //next k
} //next l
for (l = i+1; l < m; l++) {
    for (k = 0; k < j; k++) {
        B[l-1][k] = A[l][k];
    } //next k
    for (k = j+1; k < n; k++) {
        B[l-1][k-1] = A[l][k];
    } //next k
} //next l
 
return; };
 
// -------------------------------------------------------------------------- //
template <class T, size_t m, size_t n>
void complement(
    int                                          i,
    int                                          j,
    std::array< std::array< T, n >, m >         &A,
    std::array< std::array<T, n-1>, m-1>        &B
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
// none
 
// Counters
int            l, k;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
if (m == 0) { return; }
if (n == 0) { return; }
if ((i >= (long) m) || (i < 0)) {
    return;
}
if ((j >= (long) n) || (j < 0)) {
    return;
}
 
// ========================================================================== //
// EXTRACT COMPLEMENT                                                         //
// ========================================================================== //
for (l = 0; l < i; l++) {
    for (k = 0; k < j; k++) {
        B[l][k] = A[l][k];
    } //next k
    for (k = j+1; k < (long) n; k++) {
        B[l][k-1] = A[l][k];
    } //next k
} //next l
for (l = i+1; l < (long) m; l++) {
    for (k = 0; k < j; k++) {
        B[l-1][k] = A[l][k];
    } //next k
    for (k = j+1; k < (long) n; k++) {
        B[l-1][k-1] = A[l][k];
    } //next k
} //next l
 
return; };
 
// -------------------------------------------------------------------------- //
template <class T>
void zeros(
    std::vector< std::vector < T > >            &A,
    int                                          m,
    int                                          n
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
// none
 
// Counters
int      i, j;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
if ((m == 0) || (n == 0)) {
    return;
}
 
// ========================================================================== //
// CREATE MATRIX                                                              //
// ========================================================================== //
A.resize(m);
for (i = 0; i < m; i++) {
    A[i].resize(n, (T) 0.0);
    for (j = 0; j < n; j++) {
        A[i][j] = (T) 0.0;
    } //next j
} //next i
 
return; };
 
// -------------------------------------------------------------------------- //
template <class T, size_t m, size_t n>
void zeros(
    std::array< std::array < T, n >, m >        &A
) {
 
// ========================================================================== //
// template <class T, size_t m, size_t n>                                     //
// void zeros(                                                                //
//     array< array < T, n >, m >  &A)                                        //
//                                                                            //
// Initialize a m-by-n matrix of zeros.                                       //
// ========================================================================== //
// INPUT                                                                      //
// ========================================================================== //
// - A      : array< array< T, n >, m >, with m-by-n matrix of zeros          //
// - none                                                                     //
// ========================================================================== //
// OUTPUT                                                                     //
// ========================================================================== //
// - none                                                                     //
// ========================================================================== //
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
// none
 
// Counters
size_t   i;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
if ((m == 0) || (n == 0)) {
    return;
}
 
// ========================================================================== //
// CREATE MATRIX                                                              //
// ========================================================================== //
for (i = 0; i < m; i++) {
    A[i].fill(0.0);
} //next i
 
return; };
 
// -------------------------------------------------------------------------- //
template <class T>
void ones(
    std::vector< std::vector < T > >            &A,
    int                                          m,
    int                                          n
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
// none
 
// Counters
int      i, j;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
if ((m == 0) || (n == 0)) {
    std::cout << "ERROR: number of rows (columns) must be > 0!!" << std::endl;
}
 
// ========================================================================== //
// CREATE MATRIX                                                              //
// ========================================================================== //
A.resize(m);
for (i = 0; i < m; i++) {
    A[i].resize(n, (T) 1.0);
    for (j = 0; j < n; j++) {
        A[i][j] = (T) 1.0;
    } //next j
} //next i
 
return; };
 
// -------------------------------------------------------------------------- //
template <class T, size_t m, size_t n>
void ones(
    std::array< std::array < T, n >, m >        &A
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
// none
 
// Counters
size_t   i;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
if ((m == 0) || (n == 0)) {
    return;
}
 
// ========================================================================== //
// CREATE MATRIX                                                              //
// ========================================================================== //
for (i = 0; i < m; i++) {
    A[i].fill(1.0);
} //next i
 
return; };
 
// -------------------------------------------------------------------------- //
template <class T>
void eye(
    std::vector< std::vector < T > >            &A,
    int                                          m,
    int                                          n
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
int      s = std::min(m, n);
 
// Counters
int      i;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
if ((m == 0) || (n == 0)) {
    return;
}
 
// ========================================================================== //
// CREATE MATRIX                                                              //
// ========================================================================== //
zeros(A, m, n);
for (i = 0; i < s; i++) {
    A[i][i] = (T) 1.0;
} //next i
 
return; };
 
// -------------------------------------------------------------------------- //
template <class T, size_t m, size_t n>
void eye(
    std::array< std::array < T, n >, m >        &A
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
int      s = std::min(m, n);
 
// Counters
int      i;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
if ((m == 0) || (n == 0)) {
    return;
}
 
// ========================================================================== //
// CREATE MATRIX                                                              //
// ========================================================================== //
zeros(A);
for (i = 0; i < s; i++) {
    A[i][i] = (T) 1.0;
} //next i
 
return; };
 
// Matrix determinant ======================================================= //
 
// -------------------------------------------------------------------------- //
template <class T>
T det(
    std::array< std::array< T, 1 >, 1 >         &A
) {
 
return(A[0][0]); };
 
// -------------------------------------------------------------------------- //
template <class T>
T det(
    std::vector< std::vector < T > >            &A
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
int                     m, n;
T                       d = (T) 1.0e+18;
std::vector< std::vector < T > >  C;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
m = A.size();
if (m == 0) { return(d); };
n = A[0].size();
if (n == 0) { return(d); };
if (m != n) {
    return(d);
}
 
// ========================================================================== //
// COMPUTE DETERMINANT                                                        //
// ========================================================================== //
if (m == 1) {
    d = A[0][0];
}
else if (m == 2) {
    d = A[0][0]*A[1][1] - A[0][1]*A[1][0];
}
else if (m == 3) {
    d = A[0][0]*A[1][1]*A[2][2] - A[0][0]*A[1][2]*A[2][1]
        + A[0][1]*A[1][2]*A[2][0] - A[0][1]*A[1][0]*A[2][2]
        + A[0][2]*A[1][0]*A[2][1] - A[0][2]*A[1][1]*A[2][0];
}
else {
    d = (T) 0.0;
    for (int i = 0; i < m; i++) {
        complement(0, i, A, C);
        d += pow(-1.0, i) * A[0][i] * det(C);
    } //next i
}
 
return(d); };
 
// -------------------------------------------------------------------------- //
template <class T, size_t m, size_t n>
T det(
    std::array< std::array < T, n >, m >        &A
) {
 
// ========================================================================== //
// VARIABLES DECLARATION                                                      //
// ========================================================================== //
 
// Local variables
T                       d = (T) 0.0;
 
// Counters
int                     i;
 
// ========================================================================== //
// CHECK INPUT                                                                //
// ========================================================================== //
if (m == 0) { return(d); };
if (n == 0) { return(d); };
if (m != n) {
    return(d);
}
 
// ========================================================================== //
// COMPUTE DETERMINANT                                                        //
// ========================================================================== //
if (m == 1) {
    d = A[0][0];
}
else if (m == 2) {
    d = A[0][0]*A[1][1] - A[0][1]*A[1][0];
}
else if (m == 3) {
    d = A[0][0]*A[1][1]*A[2][2] - A[0][0]*A[1][2]*A[2][1]
        + A[0][1]*A[1][2]*A[2][0] - A[0][1]*A[1][0]*A[2][2]
        + A[0][2]*A[1][0]*A[2][1] - A[0][2]*A[1][1]*A[2][0];
}
else {
    std::array< std::array < double, m-1 >, m-1 >     C;
    for (i = 0; i < m; i++) {
        complement(0, i, A, C);
        d += pow(-1.0, i) * A[0][i] * det(C);
    } //next i
}
 
return(d); };
 
 
}
 
}