floatmatrix.h

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2025   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  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 program 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */
 
#ifndef flotmtrx_h
#define flotmtrx_h
 
#include "oofemenv.h"
#include "contextioresulttype.h"
#include "contextmode.h"
#include "numerics.h"
 
#include <vector>
#include <iosfwd>
#include <initializer_list>
#include <algorithm>
#include <string>
 
namespace oofem {
class FloatArray;
template<std::size_t N> class FloatArrayF;
template<std::size_t M, std::size_t N> class FloatMatrixF;
class IntArray;
class DataStream;

 
class OOFEM_EXPORT FloatMatrix
#ifdef _USE_EIGEN
    : public MatrixXXd
#endif
{
#ifdef _USE_EIGEN
    public:
        OOFEM_EIGEN_DERIVED(FloatMatrix,MatrixXXd);
        FloatMatrix(int r, int c): MatrixXXd(r,c) {}
#else
protected:
    std::size_t nRows;
    std::size_t nColumns;
    std :: vector< double >values;
 
public:
    constexpr static int Dim = 2;
    typedef double Scalar;
    FloatMatrix(std::size_t  n, std::size_t m) : nRows(n), nColumns(m), values(n * m) {}
    FloatMatrix() : nRows(0), nColumns(0), values() {}
    FloatMatrix(const FloatMatrix &mat) : nRows(mat.nRows), nColumns(mat.nColumns), values(mat.values) {}
    FloatMatrix(FloatMatrix && mat) noexcept : nRows(mat.nRows), nColumns(mat.nColumns), values( std :: move(mat.values) ) {}
    FloatMatrix &operator=(const FloatMatrix &mat) {
        nRows = mat.nRows;
        nColumns = mat.nColumns;
        values = mat.values;
        return * this;
    }
    FloatMatrix &operator=(FloatMatrix && mat) noexcept {
        nRows = std :: move(mat.nRows);
        nColumns = std :: move(mat.nColumns);
        values = std :: move(mat.values);
        return * this;
    }
    inline int rows() const { return (int)nRows; }
    inline int cols() const { return (int)nColumns; }

    inline double &operator()(std::size_t i, std::size_t j)
    {
        #ifndef NDEBUG
            this->checkBounds(i + 1, j + 1);
        #endif
        return values [ j * nRows + i ];
    }

    inline double operator()(std::size_t i, std::size_t j) const
    {
        #ifndef NDEBUG
        this->checkBounds(i + 1, j + 1);
        #endif
        return values [ j * nRows + i ];
    }
    inline const double *data() const { return values.data(); }
    inline double *data() { return values.data(); }
#endif
 
    template<std::size_t N, std::size_t M>
    void assignFloatMatrixF(const FloatMatrixF<N,M> &mat){
        resize(N,M); // FIXME: useless 0-write
        for(Index c=0; c<cols(); c++) for(Index r=0; r<rows(); r++) (*this)(r,c)=mat(r,c);
    }
 
    void _resize_internal(int nr, int nc);
    FloatMatrix(std :: initializer_list< std :: initializer_list< double > > ini){ (*this)=FloatMatrix::fromIniList(ini); }
    FloatMatrix &operator=(std :: initializer_list< std :: initializer_list< double > >mat){ (*this)=fromIniList(mat); return *this; }
    static FloatMatrix fromIniList(std :: initializer_list< std :: initializer_list< double > >);
    static FloatMatrix fromCols(std :: initializer_list< FloatArray >mat);
    template<std::size_t M, std::size_t N>
    FloatMatrix(const FloatMatrixF<N,M> &src) { assignFloatMatrixF(src); }

    static FloatMatrix fromArray(const FloatArray &vector, bool transpose = false);
 
    template<std::size_t N, std::size_t M>
    FloatMatrix &operator=(const FloatMatrixF<N,M> &mat) {
        this->assignFloatMatrixF(mat);
        return *this;
    }

    void checkBounds(Index i, Index j) const;
    //[[deprecated("use rows() instead")]]
    inline int giveNumberOfRows() const { return rows(); }
    //[[deprecated("use cols() instead")]]
    inline int giveNumberOfColumns() const { return cols(); }
 
 

    inline bool isSquare() const { return rows() == cols(); }
    inline bool isNotEmpty() const { return rows() > 0 && cols() > 0; }

    bool isAllFinite() const;

    inline double at(std::size_t i, std::size_t j) const
    {
        return (*this)(i-1,j-1);
    }

    inline double &at(std::size_t i, std::size_t j)
    {
        return (*this)(i-1,j-1);
    }

    void assemble(const FloatMatrix &src, const IntArray &loc);
    void assemble(const FloatMatrix &src, const IntArray &rowind, const IntArray &colind);
    void assembleT(const FloatMatrix &src, const IntArray &rowind, const IntArray &colind);

    double computeFrobeniusNorm() const;
    double computeNorm(char p) const;
    double computeReciprocalCondition(char p = '1') const;
    /*
     * Computes the eigenvalues of a symmetric matrix.
     * The receiver must be square and symmetric.
     * @param lambda Eigenvalues.
     * @param v Eigenvectors (stored column wise).
     * @param neigs If set, only neigs largest eigenvalues are computed.
     * @return True if successful.
     */
    //bool computeEigenValuesSymmetric(FloatArray &lambda, FloatMatrix &v, int neigs = 0) const;
    void beDiagonal(const FloatArray &diag);
    double giveTrace() const;
    double giveDeterminant() const;

    void zero();
    void beUnitMatrix();
    void bePinvID();
    void beTranspositionOf(const FloatMatrix &src);
    void beProductOf(const FloatMatrix &a, const FloatMatrix &b);
    void addProductOf(const FloatMatrix &a, const FloatMatrix &b);
    void addTProductOf(const FloatMatrix &a, const FloatMatrix &b);
    void beTProductOf(const FloatMatrix &a, const FloatMatrix &b);
    void beProductTOf(const FloatMatrix &a, const FloatMatrix &b);
    void beDyadicProductOf(const FloatArray &vec1, const FloatArray &vec2);
    void beNMatrixOf(const FloatArray &n, int nsd);
    void beLocalCoordSys(const FloatArray &normal);
    void setSubMatrix(const FloatMatrix &src, int sr, int sc);
    void setTSubMatrix(const FloatMatrix &src, int sr, int sc);
    void beSubMatrixOf(const FloatMatrix &src, Index topRow, Index bottomRow, Index topCol, Index bottomCol);
    void beSubMatrixOf(const FloatMatrix &src, const IntArray &indxRow, const IntArray &indxCol);
    void addSubVectorRow(const FloatArray &src, int sr, int sc);
    void addSubVectorCol(const FloatArray &src, int sr, int sc);

    void copySubVectorRow(const FloatArray &src, int sr, int sc);
    void setColumn(const FloatArray &src, int c);
    void copyColumn(FloatArray &dest, int c) const;
    bool beInverseOf(const FloatMatrix &src);
    bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose = false);
    bool solveForRhs(const FloatMatrix &B, FloatMatrix &answer, bool transpose = false);
    void plusProductSymmUpper(const FloatMatrix &a, const FloatMatrix &b, double dV);
    void plusDyadSymmUpper(const FloatArray &a, double dV);
    void plusProductUnsym(const FloatMatrix &a, const FloatMatrix &b, double dV);
    void plusDyadUnsym(const FloatArray &a, const FloatArray &b, double dV);

    void plus_Nt_a_otimes_b_B(const FloatMatrix &N, const FloatArray &a, const FloatArray &b, const FloatMatrix &B, double dV = 1);

    void add(const FloatMatrix &a);
    void add(double s, const FloatMatrix &a);
    void subtract(const FloatMatrix &a);
    void times(double f);
    void negated();
    void symmetrized();
    void rotatedWith(const FloatMatrix &r, char mode = 'n');
    #ifndef _USE_EIGEN
        void resize(Index rows, Index cols);
    #endif
    void resizeWithData(Index, Index);

    void clear() { resize(0,0); }
    bool jaco_(FloatArray &eval, FloatMatrix &v, int nf);

    void printYourself() const;
    void printYourself(const std::string &name) const;
    
    void printYourselfToFile(const std::string filename, const bool showDimensions=true) const;
    
    void pY() const;

    void writeCSV(const std :: string &name) const;

    inline const double *givePointer() const { return data(); }
    inline double *givePointer() { return data(); }

    void beMatrixFormOfStress(const FloatArray &aArray);
    void beMatrixForm(const FloatArray &aArray);

    void changeComponentOrder();
 
    // Overloaded methods:
    contextIOResultType storeYourself(DataStream &stream) const;
    contextIOResultType restoreYourself(DataStream &stream);
    int givePackSize(DataStream &buff) const;
 
    friend std :: ostream &operator<<(std :: ostream &out, const FloatMatrix &r);
 
}; // class FloatMatrix
 
#ifndef _USE_EIGEN
    //@name operators

    OOFEM_EXPORT FloatMatrix &operator *= ( FloatMatrix & x, const double & a );
    OOFEM_EXPORT FloatMatrix operator * ( const FloatMatrix & x, const double & a );
    OOFEM_EXPORT FloatMatrix operator * ( const double & a, const FloatMatrix & x);
    OOFEM_EXPORT FloatMatrix operator *( const FloatMatrix & a, const FloatMatrix & b ) ;
    OOFEM_EXPORT FloatArray operator *( const FloatMatrix & a, const FloatArray & b ) ;
    OOFEM_EXPORT FloatMatrix operator +( const FloatMatrix & a, const FloatMatrix & b ) ;
    OOFEM_EXPORT FloatMatrix operator -( const FloatMatrix & a, const FloatMatrix & b ) ;
    OOFEM_EXPORT FloatMatrix &operator += ( FloatMatrix & a, const FloatMatrix & b ) ;
    OOFEM_EXPORT FloatMatrix &operator -= ( FloatMatrix & a, const FloatMatrix & b ) ;
#endif
 
 
} // end namespace oofem
#endif // flotmtrx_h