floatmatrix.h

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2013   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
 */

/*
 * The original idea for this class comes from
 * Dubois-Pelerin, Y.: "Object-Oriented  Finite Elements: Programming concepts and Implementation",
 * PhD Thesis, EPFL, Lausanne, 1992.
 */

#ifndef flotmtrx_h
#define flotmtrx_h

#include "oofemcfg.h"
#include "contextioresulttype.h"
#include "contextmode.h"

#include <vector>
#include <iosfwd>
#include <initializer_list>
#include <algorithm>

#ifdef BOOST_PYTHON
namespace boost {
namespace python {
namespace api {
class object;
};
};
};
#endif

namespace oofem {
class FloatArray;
class IntArray;
class DataStream;

class OOFEM_EXPORT FloatMatrix
{
protected:
    int nRows;
    int nColumns;
    std :: vector< double >values;

public:
    FloatMatrix(int n, int m) : nRows(n), nColumns(m), values(n * m) {}
    FloatMatrix() : nRows(0), nColumns(0), values() {}
    FloatMatrix(const FloatArray &vector, bool transpose = false);
    FloatMatrix(const FloatMatrix &mat) : nRows(mat.nRows), nColumns(mat.nColumns), values(mat.values) {}
    FloatMatrix(FloatMatrix && mat) : nRows(mat.nRows), nColumns(mat.nColumns), values( std :: move(mat.values) ) {}
    FloatMatrix(std :: initializer_list< std :: initializer_list< double > >mat);
    FloatMatrix &operator=(std :: initializer_list< std :: initializer_list< double > >mat);
    FloatMatrix &operator=(std :: initializer_list< FloatArray >mat);
    FloatMatrix &operator=(const FloatMatrix &mat) {
        nRows = mat.nRows;
        nColumns = mat.nColumns;
        values = mat.values;
        return * this;
    }
    FloatMatrix &operator=(FloatMatrix && mat) {
        nRows = std :: move(mat.nRows);
        nColumns = std :: move(mat.nColumns);
        values = std :: move(mat.values);
        return * this;
    }
    ~FloatMatrix() {}

    void checkBounds(int i, int j) const;
    inline int giveNumberOfRows() const { return nRows; }
    inline int giveNumberOfColumns() const { return nColumns; }
    inline bool isSquare() const { return nRows == nColumns; }
    inline bool isNotEmpty() const { return nRows > 0 && nColumns > 0; }

    bool isFinite() const;

#ifdef DEBUG
    double at(int i, int j) const;
#else
    inline double at(int i, int j) const { return values [ ( j - 1 ) * nRows + i - 1 ]; }
#endif

#ifdef DEBUG
    double &at(int i, int j);
#else
    inline double &at(int i, int j) { return values [ ( j - 1 ) * nRows + i - 1 ]; }
#endif

#ifdef DEBUG
    double &operator()(int i, int j);
#else
    inline double &operator()(int i, int j) { return values [ j * nRows + i ]; }
#endif

#ifdef DEBUG
    double operator()(int i, int j) const;
#else
    inline double operator()(int i, int j) const { return values [ j * nRows + i ]; }
#endif

    void assemble(const FloatMatrix &src, const IntArray &loc);
    void assemble(const FloatMatrix &src, const IntArray &rowind, const IntArray &colind);
    void assemble(const FloatMatrix &src, const int *rowind, const int *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, int topRow, int bottomRow, int topCol, int 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;
    void beInverseOf(const FloatMatrix &src);
    bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose = false);
    void 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 add(const FloatMatrix &a);
    void add(double s, const FloatMatrix &a);
    void subtract(const FloatMatrix &a);
    void times(double f);
    void negated();
    void initFromVector(const FloatArray &vector, bool transposed);
    void symmetrized();
    void rotatedWith(const FloatMatrix &r, char mode = 'n');
    void resize(int rows, int cols);
    void resizeWithData(int, int);
    void hardResize(int r, int c);
    void clear() {
        this->nRows = 0;
        this->nColumns = 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 values.data(); }
    inline double *givePointer() { return values.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);


#ifdef BOOST_PYTHON
    void __setitem__(boost :: python :: api :: object t, double val);
    double __getitem__(boost :: python :: api :: object t);
    void beCopyOf(const FloatMatrix &src) { this->operator=(src); }
#endif
};
} // end namespace oofem
#endif // flotmtrx_h