floatmatrix.C

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
 */
 
#include "floatmatrix.h"
#include "floatarray.h"
#include "intarray.h"
#include "mathfem.h"
#include "error.h"
#include "datastream.h"
 
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ostream>
#include <iostream>
#include <fstream>
#include <iomanip>
 
 
#ifdef _USE_EIGEN
    #include<Eigen/Eigenvalues>
#endif
 
 
 
namespace oofem {
 
 
#define _LOOP_MATRIX(r,c) for(Index c=0; c<cols(); c++) for(Index r=0; r<rows(); r++)
 
#ifdef _USE_EIGEN
    void FloatMatrix::_resize_internal(int nr, int nc){
        MatrixXXd::resize(nr,nc);
    }
#else
    void FloatMatrix::_resize_internal(int nr, int nc){
        this->nRows = nr; this->nColumns = nc;
        std::size_t nsize = this->nRows * this->nColumns;
        if ( nsize < this->values.size() ) {
            this->values.resize(nsize);
        } else if ( nsize > this->values.size() ) {
            this->values.assign(nsize, 0.);
        }
    };
 
 
    void FloatMatrix :: resize(Index rows, Index columns)
    //
    // resizes receiver, all data will be lost
    //
    {
        this->nRows = rows;
        this->nColumns = columns;
        this->values.assign(rows * columns, 0.);
    }
 
#endif
 
void FloatMatrix :: resizeWithData(Index rows, Index columns)
//
// resizes receiver, all data kept
//
{
    // Check of resize if necessary at all.
    if ( rows == this->rows() && columns == this->cols() ) {
        return;
    }
 
 
 
    #ifdef _USE_EIGEN
        FloatMatrix old(*this); // inefficient, but for testing okay
        resize(rows,columns); // FIXME: useless zeroing
    #else
        FloatMatrix old( std :: move(* this) );
        this->nRows = rows;
        this->nColumns = columns;
        this->values.resize(rows * columns);
    #endif
 
    Index ii = min( rows, (Index)old.rows() );
    Index jj = min( columns, (Index)old.cols() );
    // copy old values if possible
    for (Index i = 1; i <= ii; i++ ) {
        for (Index j = 1; j <= jj; j++ ) {
            this->at(i, j) = old.at(i, j);
        }
    }
}
 
 
 
 
FloatMatrix FloatMatrix::fromIniList(std::initializer_list<std::initializer_list<double>> ini){
    FloatMatrix ret;
    ret._resize_internal( ( int ) ini.begin()->size(), ( int ) ini.size() );
    Index c=0;
    for (auto col : ini) {
        #if DEBUG
            if ( ret.rows() != (Index) col.size() ) {
                    OOFEM_ERROR("Initializer list has inconsistent column sizes.");
            }
        #endif
        Index r=0;
        for(const double& x: col) ret(r++,c)=x;
        c++;
    }
    return ret;
}
 
 
FloatMatrix FloatMatrix :: fromCols ( std :: initializer_list< FloatArray >mat )
{
    FloatMatrix ret( mat.begin()->giveSize(), ( int ) mat.size() );
    Index c=0;
    for ( auto col : mat ) {
        #if DEBUG
            if ( ret.rows() != (Index) col.size() ) {
                OOFEM_ERROR("Columns have inconsistent column sizes.");
            }
        #endif
        for(Index r=0; r<col.size(); r++) ret(r,c)=col(r);
        c++;
    }
 
    return ret;
}
 
void FloatMatrix :: checkBounds(Index i, Index j) const
// Checks that the receiver includes a position (i,j).
{
    if ( i <= 0 ) {
        OOFEM_ERROR("matrix error on rows : %d < 0", (int)i);
    }
    if ( j <= 0 ) {
        OOFEM_ERROR("matrix error on columns : %d < 0", (int)j);
    }
 
    if ( i > rows() ) {
        OOFEM_ERROR("matrix error on rows : %d > %d", (int)i, (int)rows());
    }
 
    if ( j > cols() ) {
        OOFEM_ERROR("matrix error on columns : %d > %d", (int)j, (int)cols());
    }
}
 
 
void FloatMatrix :: assemble(const FloatMatrix &src, const IntArray &loc)
{
    std::size_t ii, jj, size = src.rows();
 
#ifndef NDEBUG
    if ( size != loc.size() ) {
        OOFEM_ERROR("dimensions of 'src' and 'loc' mismatch");
    }
 
    if ( !src.isSquare() ) {
        OOFEM_ERROR("'src' is not sqaure matrix");
    }
#endif
 
    for ( std::size_t i = 1; i <= size; i++ ) {
        if ( ( ii = loc.at(i) ) ) {
            for ( std::size_t j = 1; j <= size; j++ ) {
                if ( ( jj = loc.at(j) ) ) {
                    this->at(ii, jj) += src.at(i, j);
                }
            }
        }
    }
}
 
#ifdef _USE_EIGEN___BROKEN
    void FloatMatrix :: assemble(const FloatMatrix &src, const IntArray &rowind, const IntArray &colind){
        /* rowind may contain 0s, which mean that the element is skipped */
        Eigen::Map<const Eigen::ArrayXi> rr(rowind.givePointer(),rowind.size());
        Eigen::Map<const Eigen::ArrayXi> cc(rowind.givePointer(),rowind.size());
        // auto rOk=(rr>0), cOk=(cc>0);
        Eigen::VectorXi R=Eigen::VectorXi::LinSpaced(rr.size(), 0, rr.size()-1).array();
        Eigen::VectorXi C=Eigen::VectorXi::LinSpaced(cc.size(), 0, cc.size()-1).array();
        // Eigen::VectorXi rThis=rr[rr>0];
        (*this)(rr[rr>0]-1,cc[cc>0]-1)=src(R,C);
        // auto rOk=(rr>0), cOk=(cc>0);
        //(*this)((rr-1)(rr!=0),(cc-1)(cOk))+=src(Eigen::seq(0,src.rows()(rOk),cOk);
        // FloatMatrix src2=src(rowind>0
        //(*this)(rowind.minusOne(),colind.minusOne())+=src;
    }
#else
void FloatMatrix :: assemble(const FloatMatrix &src, const IntArray &rowind, const IntArray &colind)
{
    int ii, jj;
    int nr = src.giveNumberOfRows();
    int nc = src.giveNumberOfColumns();
 
#ifndef NDEBUG
    if ( nr != rowind.giveSize() ) {
        OOFEM_ERROR("row dimensions of 'src' and 'rowind' mismatch");
    }
 
    if ( nc != colind.giveSize() ) {
        OOFEM_ERROR("column dimensions of 'src' and 'colind' mismatch");
    }
#endif
 
    for ( int i = 1; i <= nr; i++ ) {
        if ( ( ii = rowind.at(i) ) ) {
            for ( int j = 1; j <= nc; j++ ) {
                if ( ( jj = colind.at(j) ) ) {
                    this->at(ii, jj) += src.at(i, j);
                }
            }
        }
    }
}
#endif
 
#ifdef _USE_EIGEN___BROKEN
void FloatMatrix :: assembleT(const FloatMatrix &src, const IntArray &rowind, const IntArray &colind){
    (*this)(rowind.minusOne(),colind.minusOne())+=src.transpose();
}
#else
void FloatMatrix :: assembleT(const FloatMatrix &src, const IntArray &rowind, const IntArray &colind)
{
    int ii, jj;
    int nr = src.giveNumberOfRows();
    int nc = src.giveNumberOfColumns();
 
#ifndef NDEBUG
    if ( nr != rowind.giveSize() ) {
        OOFEM_ERROR("row dimensions of 'src' and 'rowind' mismatch");
    }
 
    if ( nc != colind.giveSize() ) {
        OOFEM_ERROR("column dimensions of 'src' and 'colind' mismatch");
    }
#endif
 
    for ( int i = 1; i <= nr; i++ ) {
        if ( ( ii = rowind.at(i) ) ) {
            for ( int j = 1; j <= nc; j++ ) {
                if ( ( jj = colind.at(j) ) ) {
                    this->at(jj, ii) += src.at(i, j);
                }
            }
        }
    }
}
#endif
 
 
#ifdef _USE_EIGEN
// helper routines
void _zeroedIfEmpty(FloatMatrix& A, Eigen::Index rows, Eigen::Index cols){ if(A.isNotEmpty()) return; A.setZero(rows,cols); }
void _ensureRowsCols(FloatMatrix& A, Eigen::Index minRows, Eigen::Index minCols){ if(A.rows()>=minRows && A.cols()>=minCols) return; A.resizeWithData(max(A.rows(),minRows),max(A.cols(),minCols)); }
 
 
bool FloatMatrix::isAllFinite() const { return this->array().isFinite().all(); }
void FloatMatrix :: beTranspositionOf(const FloatMatrix &src) { *this=src.transpose(); }
void FloatMatrix :: beProductOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix){ *this=aMatrix*bMatrix; }
void FloatMatrix :: beTProductOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix){ *this=aMatrix.transpose()*bMatrix; }
void FloatMatrix :: beProductTOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix){ *this=aMatrix*bMatrix.transpose(); }
void FloatMatrix :: addProductOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix){ *this+=aMatrix*bMatrix; }
void FloatMatrix :: addTProductOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix){ *this+=aMatrix.transpose()*bMatrix; }
void FloatMatrix :: beDyadicProductOf(const FloatArray &vec1, const FloatArray &vec2){ *this=vec1*vec2.transpose(); }
void FloatMatrix :: symmetrized() { *this = this->selfadjointView<Eigen::Upper>().toDenseMatrix(); }
void FloatMatrix :: times(double factor) { (*this)*=factor; }
void FloatMatrix :: negated() { (*this)*=-1; }
double FloatMatrix :: computeFrobeniusNorm() const { return this->norm(); }
double FloatMatrix :: computeNorm(char p) const {
    if(p=='1') return this->lpNorm<1>();
    OOFEM_ERROR("p == %d not implemented.\n", p);
}
void FloatMatrix::setSubMatrix(const FloatMatrix &src, int sr, int sc){ this->block(sr-1,sc-1,src.rows(),src.cols())=src; }
void FloatMatrix::setTSubMatrix(const FloatMatrix &src, int sr, int sc){ this->block(sr-1,sc-1,src.cols(),src.rows())=src.transpose(); }
void FloatMatrix :: addSubVectorRow(const FloatArray &src, int sr, int sc){ Index sr0=sr-1, sc0=sc-1; _ensureRowsCols(*this,sr0+1,sc0+src.size()+1); row(sr0).segment(sc0,src.size())=src; }
void FloatMatrix :: addSubVectorCol(const FloatArray &src, int sr, int sc){ Index sr0=sr-1, sc0=sc-1; _ensureRowsCols(*this,sr0+src.size()+1,sc0+1); col(sc0).segment(sr0,src.size())=src; }
void FloatMatrix :: setColumn(const FloatArray &src, int c){ this->col(c-1)=src; }
void FloatMatrix :: copyColumn(FloatArray &dest, int c) const { dest=this->col(c-1); }
// void FloatMatrix :: plusProductSymmUpper(const FloatMatrix &a, const FloatMatrix &b, double dV){ _zeroedIfEmpty(*this,a.cols(),b.cols()); /*...*/ } // not sure who to write this nicely with Eigen
void FloatMatrix :: plusDyadSymmUpper(const FloatArray &a, double dV)                             { _zeroedIfEmpty(*this,a.size(),a.size()); this->selfadjointView<Eigen::Upper>().rankUpdate(a,dV); }
void FloatMatrix :: plusProductUnsym(const FloatMatrix &a, const FloatMatrix &b, double dV)       { _zeroedIfEmpty(*this,a.cols(),b.cols()); *this+=a.transpose()*b*dV; }
void FloatMatrix :: plusDyadUnsym(const FloatArray &a, const FloatArray &b, double dV)            { _zeroedIfEmpty(*this,a.size(),b.size()); *this+=a*b.transpose()*dV; }
void FloatMatrix :: plus_Nt_a_otimes_b_B(const FloatMatrix &N, const FloatArray &a, const FloatArray &b, const FloatMatrix &B, double dV){ _zeroedIfEmpty(*this,N.cols(),B.cols()); *this=dV*((N*a)*(b.transpose()*B).transpose()); }
bool FloatMatrix :: beInverseOf(const FloatMatrix &src){
    Eigen::FullPivLU<Eigen::MatrixXd> lu(src);
    // lu.setThreshold(1e-30); // this is what the original code uses
    if(!lu.isInvertible()) return false;
    *this=lu.inverse();
    return true;
}
void FloatMatrix :: beSubMatrixOf(const FloatMatrix &src, Index topRow, Index bottomRow, Index topCol, Index bottomCol){ *this=src.block(topRow-1,topCol-1,bottomRow-topRow+1,bottomCol-topCol+1); }
void FloatMatrix :: beSubMatrixOf(const FloatMatrix &src, const IntArray &indxRow, const IntArray &indxCol)            { *this=src(indxRow.minusOne(),indxCol.minusOne()); }
void FloatMatrix::add(const FloatMatrix& a)          { if(!a.isNotEmpty()) return; if(!isNotEmpty()){ *this=a; return; }   *this+=a;   }
void FloatMatrix::add(double s, const FloatMatrix& a){ if(!a.isNotEmpty()) return; if(!isNotEmpty()){ *this=s*a; return; } *this+=s*a; }
void FloatMatrix :: subtract(const FloatMatrix &a)   { if(!a.isNotEmpty()) return; if(!isNotEmpty()){ *this=-a; return; }  *this-=a;   }
FloatMatrix FloatMatrix :: fromArray(const FloatArray &vector, bool transposed){ if(transposed) return vector.transpose(); return vector; }
void FloatMatrix :: zero() { this->setZero(); }
void FloatMatrix :: beUnitMatrix(){ if(!this->isSquare()) OOFEM_ERROR("cannot make unit matrix of %ld by %ld matrix", rows(), cols()); *this=Eigen::MatrixXd::setIdentity(); }
double FloatMatrix :: giveDeterminant() const { return this->determinant(); }
void FloatMatrix :: beDiagonal(const FloatArray &diag) { *this=diag.asDiagonal().toDenseMatrix(); }
double FloatMatrix :: giveTrace() const { return this->trace(); }
#endif
 
#ifndef _USE_EIGEN
 
bool FloatMatrix :: isAllFinite() const
{
    for(Index r=0; r<rows(); r++) for(Index c=0; c<cols(); c++){
        if ( !std::isfinite((*this)(r,c)) ) {
            return false;
        }
    }
 
    return true;
}
 
 
void FloatMatrix :: beTranspositionOf(const FloatMatrix &src)
{
    // receiver becomes a transposition of src
    int nrows = src.giveNumberOfColumns(), ncols = src.giveNumberOfRows();
    _resize_internal(nrows, ncols);
 
    for ( int i = 1; i <= nrows; i++ ) {
        for ( int j = 1; j <= ncols; j++ ) {
            this->at(i, j) = src.at(j, i);
        }
    }
}
 
 
void FloatMatrix :: beProductOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix)
// Receiver = aMatrix * bMatrix
{
#  ifndef NDEBUG
    if ( aMatrix.cols() != bMatrix.rows() ) {
        OOFEM_ERROR("error in product A*B : dimensions do not match, A(*,%d), B(%d,*)", aMatrix.cols(), bMatrix.rows());
    }
#  endif
    _resize_internal(aMatrix.rows(), bMatrix.cols());
#  ifdef __LAPACK_MODULE
    const int this_nColumns=cols(), this_nRows=rows();
    const int aMatrix_nColumns=aMatrix.cols(), aMatrix_nRows=aMatrix.rows();
    const int bMatrix_nColumns=bMatrix.cols(), bMatrix_nRows=bMatrix.rows();
    double alpha = 1., beta = 0.;
    dgemm_("n", "n", & this_nRows, & this_nColumns, & aMatrix_nColumns,
           & alpha, aMatrix.givePointer(), & aMatrix_nRows, bMatrix.givePointer(), & bMatrix_nRows,
           & beta, this->givePointer(), & this_nRows,
           aMatrix_nColumns, bMatrix_nColumns, this_nColumns);
#  else
    for ( Index i = 1; i <= aMatrix.rows(); i++ ) {
        for ( Index j = 1; j <= bMatrix.cols(); j++ ) {
            double coeff = 0.;
            for ( Index k = 1; k <= aMatrix.cols(); k++ ) {
                coeff += aMatrix.at(i, k) * bMatrix.at(k, j);
            }
 
            this->at(i, j) = coeff;
        }
    }
#  endif
}
 
 
void FloatMatrix :: beTProductOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix)
// Receiver = aMatrix^T * bMatrix
{
#  ifndef NDEBUG
    if ( aMatrix.rows() != bMatrix.rows() ) {
        OOFEM_ERROR("error in product A*B : dimensions do not match");
    }
#  endif
    _resize_internal(aMatrix.cols(), bMatrix.cols());
#  ifdef __LAPACK_MODULE
    double alpha = 1., beta = 0.;
    const int this_nColumns=cols(), this_nRows=rows();
    const int aMatrix_nColumns=aMatrix.cols(), aMatrix_nRows=aMatrix.rows();
    const int bMatrix_nColumns=bMatrix.cols(), bMatrix_nRows=bMatrix.rows();
    dgemm_("t", "n", & this_nRows, & this_nColumns, & aMatrix_nRows,
           & alpha, aMatrix.givePointer(), & aMatrix_nRows, bMatrix.givePointer(), & bMatrix_nRows,
           & beta, this->givePointer(), & this_nRows,
           aMatrix_nColumns, bMatrix_nColumns, this_nColumns);
#  else
    for ( Index i = 1; i <= aMatrix.cols(); i++ ) {
        for (Index j = 1; j <= bMatrix.cols(); j++ ) {
            double coeff = 0.;
            for (Index k = 1; k <= aMatrix.rows(); k++ ) {
                coeff += aMatrix.at(k, i) * bMatrix.at(k, j);
            }
 
            this->at(i, j) = coeff;
        }
    }
#endif
}
 
 
void FloatMatrix :: beProductTOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix)
// Receiver = aMatrix * bMatrix^T
{
#  ifndef NDEBUG
    if ( aMatrix.cols() != bMatrix.cols() ) {
        OOFEM_ERROR("error in product A*B : dimensions do not match");
    }
#  endif
    _resize_internal(aMatrix.rows(), bMatrix.rows());
#  ifdef __LAPACK_MODULE
    const int this_nColumns=cols(), this_nRows=rows();
    const int aMatrix_nColumns=aMatrix.cols(), aMatrix_nRows=aMatrix.rows();
    const int bMatrix_nColumns=bMatrix.cols(), bMatrix_nRows=bMatrix.rows();
    double alpha = 1., beta = 0.;
    dgemm_("n", "t", & this_nRows, & this_nColumns, & aMatrix_nColumns,
           & alpha, aMatrix.givePointer(), & aMatrix_nRows, bMatrix.givePointer(), & bMatrix_nRows,
           & beta, this->givePointer(), & this_nRows,
           aMatrix_nColumns, bMatrix_nColumns, this_nColumns);
#  else
    for (Index i = 1; i <= aMatrix.rows(); i++ ) {
        for (Index j = 1; j <= bMatrix.rows(); j++ ) {
            double coeff = 0.;
            for (Index k = 1; k <= aMatrix.cols(); k++ ) {
                coeff += aMatrix.at(i, k) * bMatrix.at(j, k);
            }
 
            this->at(i, j) = coeff;
        }
    }
#  endif
}
 
 
void FloatMatrix :: addProductOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix)
// Receiver = aMatrix * bMatrix
{
#  ifndef NDEBUG
    if ( aMatrix.cols() != bMatrix.rows() ) {
        OOFEM_ERROR("error in product A*B : dimensions do not match");
    }
    if ( aMatrix.rows() != this->rows() || bMatrix.cols() != this->cols() ) {
        OOFEM_ERROR("error in product receiver : dimensions do not match");
    }
#  endif
 
#  ifdef __LAPACK_MODULE
    const int this_nColumns=cols(), this_nRows=rows();
    const int aMatrix_nColumns=aMatrix.cols(), aMatrix_nRows=aMatrix.rows();
    const int bMatrix_nColumns=bMatrix.cols(), bMatrix_nRows=bMatrix.rows();
    double alpha = 1., beta = 1.;
    dgemm_("n", "n", & this_nRows, & this_nColumns, & aMatrix_nColumns,
           & alpha, aMatrix.givePointer(), & aMatrix_nRows, bMatrix.givePointer(), & bMatrix_nRows,
           & beta, this->givePointer(), & this_nRows,
           aMatrix_nColumns, bMatrix_nColumns, this_nColumns);
#  else
    for (Index i = 1; i <= aMatrix.rows(); i++ ) {
        for (Index j = 1; j <= bMatrix.cols(); j++ ) {
            double coeff = 0.;
            for (Index k = 1; k <= aMatrix.cols(); k++ ) {
                coeff += aMatrix.at(i, k) * bMatrix.at(k, j);
            }
 
            this->at(i, j) += coeff;
        }
    }
#  endif
}
 
void FloatMatrix :: addTProductOf(const FloatMatrix &aMatrix, const FloatMatrix &bMatrix)
// Receiver += aMatrix^T * bMatrix
{
#  ifndef NDEBUG
    if ( aMatrix.rows() != bMatrix.rows() ) {
        OOFEM_ERROR("error in product A*B : dimensions do not match");
    }
    if ( aMatrix.cols() != this->cols() || bMatrix.cols() != this->rows() ) {
        OOFEM_ERROR("error in product receiver : dimensions do not match");
    }
#  endif
 
#  ifdef __LAPACK_MODULE
    double alpha = 1., beta = 1.;
    const int this_nColumns=cols(), this_nRows=rows();
    const int aMatrix_nColumns=aMatrix.cols(), aMatrix_nRows=aMatrix.rows();
    const int bMatrix_nColumns=bMatrix.cols(), bMatrix_nRows=bMatrix.rows();
    dgemm_("t", "n",  & this_nRows, & this_nColumns, & aMatrix_nRows,
           & alpha, aMatrix.givePointer(), & aMatrix_nRows, bMatrix.givePointer(), & bMatrix_nRows,
           & beta, this->givePointer(), & this_nRows,
           aMatrix_nColumns, bMatrix_nColumns, this_nColumns);
#  else
    for (Index i = 1; i <= aMatrix.cols(); i++ ) {
        for (Index j = 1; j <= bMatrix.cols(); j++ ) {
            double coeff = 0.;
            for (Index k = 1; k <= aMatrix.rows(); k++ ) {
                coeff += aMatrix.at(k, i) * bMatrix.at(k, j);
            }
 
            this->at(i, j) += coeff;
        }
    }
#  endif
}
 
 
void FloatMatrix :: beDyadicProductOf(const FloatArray &vec1, const FloatArray &vec2)
// Receiver = vec1 * vec2^T
{
    int n1 = vec1.giveSize();
    int n2 = vec2.giveSize();
    _resize_internal(n1, n2);
    for ( int j = 1; j <= n2; j++ ) {
        for ( int i = 1; i <= n1; i++ ) {
            this->at(i, j) = vec1.at(i) * vec2.at(j);
        }
    }
}
 
 
void FloatMatrix :: setSubMatrix(const FloatMatrix &src, int sr, int sc)
{
    sr--;
    sc--;
 
    int srcRows = src.giveNumberOfRows(), srcCols = src.giveNumberOfColumns();
#ifndef NDEBUG
    int nr = sr + srcRows;
    int nc = sc + srcCols;
 
    if ( ( this->giveNumberOfRows() < nr ) || ( this->giveNumberOfColumns() < nc ) ) {
        OOFEM_ERROR("Sub matrix doesn't fit inside allocated space.");
    }
#endif
 
    // add sub-matrix
    for ( int j = 0; j < srcCols; j++ ) {
        for ( int i = 0; i < srcRows; i++ ) {
            ( * this )( sr + i, sc + j ) = src(i, j);
        }
    }
    //memcpy( &(*this)(sr + 0, sc + j), &src(0,j), srcRows * sizeof(int));
}
 
 
void FloatMatrix :: setTSubMatrix(const FloatMatrix &src, int sr, int sc)
{
    sr--;
    sc--;
 
    int srcRows = src.giveNumberOfRows(), srcCols = src.giveNumberOfColumns();
#ifndef NDEBUG
    int nr = sr + srcCols;
    int nc = sc + srcRows;
 
    if ( ( this->giveNumberOfRows() < nr ) || ( this->giveNumberOfColumns() < nc ) ) {
        OOFEM_ERROR("Sub matrix doesn't fit inside allocated space");
    }
#endif
 
    // add sub-matrix
    for ( int i = 0; i < srcCols; i++ ) {
        for ( int j = 0; j < srcRows; j++ ) {
            ( * this )( sr + i, sc + j ) = src(j, i);
        }
    }
}
 
 
 
 
void FloatMatrix :: addSubVectorRow(const FloatArray &src, int sr, int sc)
{
    sc--;
 
    int srcCols = src.giveSize();
 
    int nr = sr;
    int nc = sc + srcCols;
 
    if ( ( this->giveNumberOfRows() < nr ) || ( this->giveNumberOfColumns() < nc ) ) {
        this->resizeWithData( max(this->giveNumberOfRows(), nr), max(this->giveNumberOfColumns(), nc) );
    }
 
    // add sub-matrix
    for ( int j = 1; j <= srcCols; j++ ) {
        this->at(sr, sc + j) += src.at(j);
    }
}
 
 
 
 
void FloatMatrix :: addSubVectorCol(const FloatArray &src, int sr, int sc)
{
    sr--;
 
    int srcRows = src.giveSize();
 
    int nr = sr + srcRows;
    int nc = sc;
 
    if ( ( this->giveNumberOfRows() < nr ) || ( this->giveNumberOfColumns() < nc ) ) {
        this->resizeWithData( max(this->giveNumberOfRows(), nr), max(this->giveNumberOfColumns(), nc) );
    }
 
    // add sub-matrix
    for ( int j = 1; j <= srcRows; j++ ) {
        this->at(sr + j, sc) += src.at(j);
    }
}
 
void FloatMatrix :: setColumn(const FloatArray &src, int c)
{
    int nr = src.giveSize();
#ifndef NDEBUG
    if ( this->giveNumberOfRows() != nr || c < 1 || c > this->giveNumberOfColumns() ) {
        OOFEM_ERROR("Size mismatch");
    }
#endif
 
    for(int r=0; r<nr; r++) (*this)(r,c-1)=src(r);
}
 
 
void FloatMatrix :: copyColumn(FloatArray &dest, int c) const
{
    int nr = this->giveNumberOfRows();
#ifndef NDEBUG
    if ( c < 1 || c > this->giveNumberOfColumns() ) {
        OOFEM_ERROR("Column outside range (%d)", c);
    }
#endif
 
    dest.resize(nr);
    for(int r=0; r<nr; r++) dest(r)=(*this)(r,c-1);
}
 
void FloatMatrix :: plusDyadSymmUpper(const FloatArray &a, double dV)
{
    if ( !this->isNotEmpty() ) {
        this->resize(a.giveSize(),a.giveSize());
    }
#ifdef __LAPACK_MODULE
    int inc = 1;
    int sizeA = a.giveSize();
    const int this_nColumns=cols(), this_nRows=rows();
    dsyr_("u", & sizeA, & dV, a.givePointer(), & inc,
          this->givePointer(), & this_nRows,
          sizeA, this_nColumns);
#else
    for (Index i = 1; i <= rows(); i++ ) {
        for (Index j = i; j <= cols(); j++ ) {
            this->at(i, j) += a.at(i) * a.at(j) * dV;
        }
    }
#endif
}
 
 
 
 
void FloatMatrix :: plusProductUnsym(const FloatMatrix &a, const FloatMatrix &b, double dV)
// Adds to the receiver the product  a(transposed).b dV .
// If the receiver has a null size, it is expanded.
// Advantage : does not compute the transposition of matrix a.
{
    if ( !this->isNotEmpty() ) {
        this->resize(a.cols(),b.cols());
    }
#ifdef __LAPACK_MODULE
    const int this_nColumns=cols(), this_nRows=rows();
    const int a_nColumns=a.cols(), a_nRows=a.rows();
    const int b_nColumns=b.cols(), b_nRows=b.rows();
    double beta = 1.;
    dgemm_("t", "n", & this_nRows, & this_nColumns, & a_nRows,
           & dV, a.givePointer(), & a_nRows, b.givePointer(), & b_nRows,
           & beta, this->givePointer(), & this_nRows,
           a_nColumns, b_nColumns, this_nColumns);
#else
    for (Index i = 1; i <= rows(); i++ ) {
        for (Index j = 1; j <= cols(); j++ ) {
            double summ = 0.;
            for (Index k = 1; k <= a.rows(); k++ ) {
                summ += a.at(k, i) * b.at(k, j);
            }
 
            this->at(i, j) += summ * dV;
        }
    }
#endif
}
 
 
void FloatMatrix :: plusDyadUnsym(const FloatArray &a, const FloatArray &b, double dV)
{
    if ( !this->isNotEmpty() ) {
        this->resize(a.giveSize(),b.giveSize());
    }
#ifdef __LAPACK_MODULE
    int inc = 1;
    int sizeA = a.giveSize();
    int sizeB = b.giveSize();
    const int this_nColumns=cols(), this_nRows=rows();
    dger_(& sizeA, & sizeB, & dV, a.givePointer(), & inc,
          b.givePointer(), & inc, this->givePointer(), & this_nRows,
          sizeA, sizeB, this_nColumns);
#else
    for (Index i = 1; i <= rows(); i++ ) {
        for (Index j = 1; j <= cols(); j++ ) {
            this->at(i, j) += a.at(i) * b.at(j) * dV;
        }
    }
#endif
}
 
void FloatMatrix :: plus_Nt_a_otimes_b_B(const FloatMatrix &N, const FloatArray &a, const FloatArray &b, const FloatMatrix &B, double dV)
{
    if ( !this->isNotEmpty() ) {
      this->nRows = N.nColumns;
      this->nColumns = B.nColumns;
      this->values.assign(this->nRows * this->nColumns, 0.);
    }
    auto a_size = a.giveSize();
    auto b_size = b.giveSize();
    for (std::size_t i = 1; i <= nRows; i++ ) {
      for (std::size_t j = 1; j <= nColumns; j++ ) {
    for (int k = 1; k <= a_size; k++ ) {
      for (int l = 1; l <= b_size; l++ ) {
        this->at(i, j) += N.at(k,i) * a.at(k) * b.at(l) * B.at(l,j) * dV;
      }
    }
      }
    }
}
 
 
 
 
bool FloatMatrix :: beInverseOf(const FloatMatrix &src)
// Receiver becomes inverse of given parameter src. If necessary, size is adjusted.
{
    double det;
 
#  ifndef NDEBUG
    if ( !src.isSquare() ) {
        OOFEM_ERROR("cannot inverse a %d by %d matrix", src.rows(), src.cols());
    }
#  endif
 
    _resize_internal(src.rows(), src.cols());
    int nRows=this->rows();
 
    if ( nRows == 1 ) {
        if (fabs(src.at(1,1)) > 1.e-30) {
            this->at(1, 1) = 1. / src.at(1, 1);
            return true;
        } else {
            return false;
        }
    } else if ( nRows == 2 ) {
        det = src.at(1, 1) * src.at(2, 2) - src.at(1, 2) * src.at(2, 1);
        if (fabs(det)>1.e-30) {
            this->at(1, 1) =  src.at(2, 2) / det;
            this->at(2, 1) = -src.at(2, 1) / det;
            this->at(1, 2) = -src.at(1, 2) / det;
            this->at(2, 2) =  src.at(1, 1) / det;
            return true;
        } else {
            return false;
        }
    } else if ( nRows == 3 ) {
        det = src.at(1, 1) * src.at(2, 2) * src.at(3, 3) + src.at(1, 2) * src.at(2, 3) * src.at(3, 1) +
              src.at(1, 3) * src.at(2, 1) * src.at(3, 2) - src.at(1, 3) * src.at(2, 2) * src.at(3, 1) -
              src.at(2, 3) * src.at(3, 2) * src.at(1, 1) - src.at(3, 3) * src.at(1, 2) * src.at(2, 1);
        if (fabs(det)>1.e-30) {
            this->at(1, 1) = ( src.at(2, 2) * src.at(3, 3) - src.at(2, 3) * src.at(3, 2) ) / det;
            this->at(2, 1) = ( src.at(2, 3) * src.at(3, 1) - src.at(2, 1) * src.at(3, 3) ) / det;
            this->at(3, 1) = ( src.at(2, 1) * src.at(3, 2) - src.at(2, 2) * src.at(3, 1) ) / det;
            this->at(1, 2) = ( src.at(1, 3) * src.at(3, 2) - src.at(1, 2) * src.at(3, 3) ) / det;
            this->at(2, 2) = ( src.at(1, 1) * src.at(3, 3) - src.at(1, 3) * src.at(3, 1) ) / det;
            this->at(3, 2) = ( src.at(1, 2) * src.at(3, 1) - src.at(1, 1) * src.at(3, 2) ) / det;
            this->at(1, 3) = ( src.at(1, 2) * src.at(2, 3) - src.at(1, 3) * src.at(2, 2) ) / det;
            this->at(2, 3) = ( src.at(1, 3) * src.at(2, 1) - src.at(1, 1) * src.at(2, 3) ) / det;
            this->at(3, 3) = ( src.at(1, 1) * src.at(2, 2) - src.at(1, 2) * src.at(2, 1) ) / det;
            return true;
        } else {
            return false;
        }
    } else {
        #ifdef __LAPACK_MODULE
            const int this_nRows=rows();
            IntArray ipiv(this_nRows);
            int lwork, info;
            * this = src;
 
            // LU-factorization
            dgetrf_(& this_nRows, & this_nRows, this->givePointer(), & this_nRows, ipiv.givePointer(), & info);
            if ( info != 0 ) {
                OOFEM_WARNING("dgetrf error %d", info);
                return false;
            }
 
            // Inverse
            lwork = this_nRows * this_nRows;
            FloatArray work(lwork);
            dgetri_(& this_nRows, this->givePointer(), & this_nRows, ipiv.givePointer(), work.givePointer(), & lwork, & info);
            if ( info > 0 ) {
                OOFEM_WARNING("Singular at %d", info);
                return false;
            } else if ( info < 0 ) {
                OOFEM_ERROR("Error on input %d", info);
            }
        #else
            // size >3 ... gaussian elimination - slow but safe
            //
            double piv, linkomb;
            FloatMatrix tmp = src;
            // initialize answer to be unity matrix;
            this->zero();
            for (Index i = 1; i <= nRows; i++ ) {
                this->at(i, i) = 1.0;
            }
 
            // lower triangle elimination by columns
            for (Index i = 1; i < nRows; i++ ) {
                piv = tmp.at(i, i);
                if ( fabs(piv) < 1.e-30 ) {
                    OOFEM_WARNING("pivot (%d,%d) to close to small (< 1.e-20)", i, i);
                    return false;
                }
 
                for (Index j = i + 1; j <= nRows; j++ ) {
                    linkomb = tmp.at(j, i) / tmp.at(i, i);
                    for (Index k = i; k <= nRows; k++ ) {
                        tmp.at(j, k) -= tmp.at(i, k) * linkomb;
                    }
 
                    for (Index k = 1; k <= nRows; k++ ) {
                        this->at(j, k) -= this->at(i, k) * linkomb;
                    }
                }
            }
 
            // upper triangle elimination by columns
            for ( Index i = nRows; i > 1; i-- ) {
                piv = tmp.at(i, i);
                for ( Index j = i - 1; j > 0; j-- ) {
                    linkomb = tmp.at(j, i) / piv;
                    for ( Index k = i; k > 0; k-- ) {
                        tmp.at(j, k) -= tmp.at(i, k) * linkomb;
                    }
 
                    for ( Index k = nRows; k > 0; k-- ) {
                        // tmp -> at(j,k)-= tmp  ->at(i,k)*linkomb;
                        this->at(j, k) -= this->at(i, k) * linkomb;
                    }
                }
            }
 
            // diagonal scaling
            for ( Index i = 1; i <= nRows; i++ ) {
                for ( Index j = 1; j <= nRows; j++ ) {
                    this->at(i, j) /= tmp.at(i, i);
                }
            }
        #endif
        return true;
    }
}
 
 
void FloatMatrix :: beSubMatrixOf(const FloatMatrix &src,
                                  Index topRow, Index bottomRow, Index topCol, Index bottomCol)
/*
 * modifies receiver to be  submatrix of the src matrix
 * size of receiver  submatrix is determined from
 * input parameters
 */
{
#ifndef NDEBUG
    if ( ( topRow < 1 ) || ( bottomRow < 1 ) || ( topCol < 1 ) || ( bottomCol < 1 ) ) {
        OOFEM_ERROR("subindexes size mismatch");
    }
 
    if ( ( src.rows() < bottomRow ) || ( src.cols() < bottomCol ) || ( ( bottomRow - topRow ) > src.rows() ) ||
         ( ( bottomCol - topCol ) > src.cols() ) ) {
        OOFEM_ERROR("subindexes size mismatch");
    }
#endif
 
 
    Index topRm1, topCm1;
    topRm1 = topRow - 1;
    topCm1 = topCol - 1;
 
    // allocate return value
    this->resize(bottomRow - topRm1, bottomCol - topCm1);
    for (Index i = topRow; i <= bottomRow; i++ ) {
        for (Index j = topCol; j <= bottomCol; j++ ) {
            this->at(i - topRm1, j - topCm1) = src.at(i, j);
        }
    }
}
 
 
 
void
FloatMatrix :: beSubMatrixOf(const FloatMatrix &src, const IntArray &indxRow, const IntArray &indxCol)
/*
 * Modifies receiver to be a sub-matrix of the src matrix.
 * sub-matrix has size(indxRow) x size(indxCol) with values given as
 * this(i,j) = src( indxRow(i), indxCol(j) )
 */
{
#  ifndef NDEBUG
    if ( indxRow.maximum() > src.giveNumberOfRows()  ||  indxCol.maximum() > src.giveNumberOfColumns()  ||
         indxRow.minimum() < 1  ||  indxCol.minimum() < 1 ) {
        OOFEM_ERROR("index exceeds source dimensions");
    }
# endif
 
    int szRow = indxRow.giveSize();
    int szCol = indxCol.giveSize();
    this->resize(szRow, szCol);
 
    for ( int i = 1; i <= szRow; i++ ) {
        for ( int j = 1; j <= szCol; j++ ) {
            this->at(i, j) = src.at( indxRow.at(i), indxCol.at(j) );
        }
    }
}
 
void FloatMatrix :: add(const FloatMatrix &aMatrix)
// Adds aMatrix to the receiver. If the receiver has a null size,
// adjusts its size to that of aMatrix. Returns the modified receiver.
{
    if ( aMatrix.rows() == 0 || aMatrix.cols() == 0 ) {
        return;
    }
 
    if ( !this->isNotEmpty() ) {
        this->operator = ( aMatrix );
        return;
    }
#     ifndef NDEBUG
    if ( ( aMatrix.rows() != rows() || aMatrix.cols() != cols() ) && aMatrix.isNotEmpty() ) {
        OOFEM_ERROR("dimensions mismatch : (r1,c1)+(r2,c2) : (%d,%d)+(%d,%d)", rows(), cols(), aMatrix.rows(), aMatrix.cols());
    }
#     endif
 
#ifdef __LAPACK_MODULE
    int aSize = aMatrix.rows() * aMatrix.cols();
    int inc = 1;
    double s = 1.;
    daxpy_(& aSize, & s, aMatrix.givePointer(), & inc, this->givePointer(), & inc, aSize, aSize);
#else
    _LOOP_MATRIX(r,c){ (*this)(r,c)+=aMatrix(r,c); }
#endif
}
 
 
void FloatMatrix :: add(double s, const FloatMatrix &aMatrix)
// Adds aMatrix to the receiver. If the receiver has a null size,
// adjusts its size to that of aMatrix. Returns the modified receiver.
{
    if ( aMatrix.rows() == 0 || aMatrix.cols() == 0 ) {
        return;
    }
 
    if ( !this->isNotEmpty() ) {
        this->operator = ( aMatrix );
        this->times(s);
        return;
    }
#     ifndef NDEBUG
    if ( ( aMatrix.rows() != rows() || aMatrix.cols() != cols() ) && aMatrix.isNotEmpty() ) {
        OOFEM_ERROR("dimensions mismatch : (r1,c1)+(r2,c2) : (%d,%d)+(%d,%d)", rows(), cols(), aMatrix.rows(), aMatrix.cols());
    }
#     endif
 
#ifdef __LAPACK_MODULE
    int aSize = aMatrix.rows() * aMatrix.cols();
    int inc = 1;
    daxpy_(& aSize, & s, aMatrix.givePointer(), & inc, this->givePointer(), & inc, aSize, aSize);
#else
    _LOOP_MATRIX(r,c){ (*this)(r,c)+=s*aMatrix(r,c); }
#endif
}
 
 
void FloatMatrix :: subtract(const FloatMatrix &aMatrix)
// Adds aMatrix to the receiver. If the receiver has a null size,
// adjusts its size to that of aMatrix. Returns the modified receiver.
{
    if ( !this->isNotEmpty() ) {
        this->operator = ( aMatrix );
        this->negated();
        return;
    }
#     ifndef NDEBUG
    if ( ( aMatrix.rows() != rows() || aMatrix.cols() != cols() ) && aMatrix.isNotEmpty() ) {
        OOFEM_ERROR("dimensions mismatch : (r1,c1)-(r2,c2) : (%d,%d)-(%d,%d)", rows(), cols(), aMatrix.rows(), aMatrix.cols());
    }
#     endif
 
#ifdef __LAPACK_MODULE
    int aSize = aMatrix.rows() * aMatrix.cols();
    int inc = 1;
    double s = -1.;
    daxpy_(& aSize, & s, aMatrix.givePointer(), & inc, this->givePointer(), & inc, aSize, aSize);
#else
    _LOOP_MATRIX(r,c){ (*this)(r,c)-=aMatrix(r,c); }
#endif
}
 
 
FloatMatrix FloatMatrix :: fromArray(const FloatArray &vector, bool transposed)
//
// constructor : creates (vector->giveSize(),1) FloatMatrix
// if transpose = 1 creates (1,vector->giveSize()) FloatMatrix
//
{
    FloatMatrix ret;
    if ( transposed ) {
        ret.resize(1,vector.size());
        for(int c=0; c<ret.cols(); c++) ret(0,c)=vector(c);
    } else {
        ret.resize(vector.size(),1);
        for(int r=0; r<ret.rows(); r++) ret(r,0)=vector(r);
    }
    return ret;
}
 
 
void FloatMatrix :: zero()
{
    _LOOP_MATRIX(r,c) (*this)(r,c)=0.;
}
 
 
void FloatMatrix :: beUnitMatrix()
{
    #ifndef NDEBUG
    if ( !this->isSquare() ) {
        OOFEM_ERROR("cannot make unit matrix of %d by %d matrix", rows(), cols());
    }
    #endif
 
    this->zero();
    for (Index i = 1; i <= rows(); i++ ) {
        this->at(i, i) = 1.0;
    }
}
 
 
double FloatMatrix :: giveDeterminant() const
// Returns the determinant of the receiver.
{
    #  ifndef NDEBUG
    if ( !this->isSquare() ) {
        OOFEM_ERROR("cannot compute the determinant of a non-square %d by %d matrix", rows(), cols());
    }
    #  endif
    const FloatMatrix& m(*this);
    if ( rows() == 1 ) {
        return m(0,0);
    } else if ( rows() == 2 ) {
        // return ( values [ 0 ] * values [ 3 ] - values [ 1 ] * values [ 2 ] );
        return m(0,0)*m(1,1)-m(1,0)*m(0,1);
        // return ( values [ 0 ] * values [ 3 ] - values [ 1 ] * values [ 2 ] );
    } else if ( rows() == 3 ) {
        return
        m(0,0)*m(1,1)*m(2,2)+m(1,0)*m(2,1)*m(0,2)+m(2,0)*m(0,1)*m(1,2)
        -m(0,2)*m(1,1)*m(2,0)-m(1,2)*m(2,1)*m(0,0)-m(2,2)*m(0,1)*m(1,0);
    }
 
    OOFEM_ERROR("sorry, cannot compute the determinant of a matrix larger than 3x3");
    // return 0.;
}
 
 
void FloatMatrix :: beDiagonal(const FloatArray &diag)
{
    int n = diag.giveSize();
    this->resize(n, n);
    for ( int i = 0; i < n; ++i ) {
        (*this)(i, i) = diag[i];
    }
}
 
 
double FloatMatrix :: giveTrace() const
{
    #  ifndef NDEBUG
    if ( !this->isSquare() ) {
        OOFEM_ERROR("cannot compute the trace of a non-square %d by %d matrix", rows(), cols());
    }
    #  endif
    double answer = 0.;
    for (Index k = 0; k < rows(); k++ ) {
        answer += (*this)(k,k);;
    }
    return answer;
}
 
 
 
void FloatMatrix :: symmetrized()
// Initializes the lower half of the receiver to the upper half.
{
    #  ifndef NDEBUG
    if ( rows() != cols() ) {
        OOFEM_ERROR("cannot symmetrize a non-square matrix");
    }
 
    #   endif
 
    for (Index i = 2; i <= rows(); i++ ) {
        for (Index j = 1; j < i; j++ ) {
            this->at(i, j) = this->at(j, i);
        }
    }
}
 
 
void FloatMatrix :: times(double factor)
// Multiplies every coefficient of the receiver by factor. Answers the
// modified receiver.
{
    _LOOP_MATRIX(r,c) (*this)(r,c)*=factor;
    // dscal_ seemed to be slower for typical usage of this function.
}
 
 
void FloatMatrix :: negated()
{
    _LOOP_MATRIX(r,c) (*this)(r,c)=-(*this)(r,c);
}
 
 
double FloatMatrix :: computeFrobeniusNorm() const
{
    // return sqrt( std :: inner_product(this->values.begin(), this->values.end(), this->values.begin(), 0.) );
    double ret=0;
    _LOOP_MATRIX(r,c) ret+=(*this)(r,c)*(*this)(r,c);
    return std::sqrt(ret);
}
 
 
double FloatMatrix :: computeNorm(char p) const
{
    #  ifdef __LAPACK_MODULE
    const int this_nRows=rows(), this_nColumns=cols();
    FloatArray work( this_nRows );
    int lda = max(this_nRows, 1);
    double norm = dlange_(& p, & this_nRows, & this_nColumns, this->givePointer(), & lda, work.givePointer(), 1);
    return norm;
 
    #  else
    if ( p == '1' ) { // Maximum absolute column sum.
        double col_sum, max_col = 0.0;
        for (Index j = 1; j <= this->cols(); j++ ) {
            col_sum  = 0.0;
            for (Index i = 1; i <= this->rows(); i++ ) {
                col_sum += fabs( this->at(i, j) );
            }
            if ( col_sum > max_col ) {
                max_col = col_sum;
            }
        }
        return max_col;
    }
    /*else if (p == '2') {
     *  double lambda_max;
     *  FloatMatrix AtA;
     *  FloatArray eigs;
     *  AtA.beTProductOf(*this,this);
     *  Ata.eigenValues(eigs, 1);
     *  return sqrt(eigs(0));
     * } */
    else {
        OOFEM_ERROR("p == %d not implemented.\n", p);
        // return 0.0;
    }
    #  endif
}
 
#endif
 
 
 
void FloatMatrix :: copySubVectorRow(const FloatArray &src, int sr, int sc)
{
    sc--;
 
    int srcCols = src.giveSize();
 
    int nr = sr;
    int nc = sc + srcCols;
 
    if ( ( this->giveNumberOfRows() < nr ) || ( this->giveNumberOfColumns() < nc ) ) {
        this->resizeWithData( max(this->giveNumberOfRows(), nr), max(this->giveNumberOfColumns(), nc) );
    }
 
    // add sub-matrix
    for ( int j = 1; j <= srcCols; j++ ) {
        this->at(sr, sc + j) = src.at(j);
    }
}
 
 
 
void FloatMatrix :: plusProductSymmUpper(const FloatMatrix &a, const FloatMatrix &b, double dV)
// Adds to the receiver the product  a(transposed).b dV .
// The receiver size is adjusted, if necessary.
// This method assumes that both the receiver and the product above are
// symmetric matrices, and therefore computes only the upper half of the
// receiver ; the lower half is not modified. Other advantage : it does
// not compute the transposition of matrix a.
{
    if ( !this->isNotEmpty() ) {
        this->resize(a.cols(),b.cols());
    }
 
    #ifdef __LAPACK_MODULE
    const int this_nColumns=cols(), this_nRows=rows();
    const int a_nColumns=a.cols(), a_nRows=a.rows();
    const int b_nColumns=b.cols(), b_nRows=b.rows();
    double beta = 1.;
    if ( this_nRows < 20 ) {
        // Split the matrix into 2 columns, s1 + s2 = n ( = nRows = nColumns ).
        int s1 = this_nRows / 2;
        int s2 = this_nRows - s1;
        // First column block, we only take the first s rows by only taking the first s columns in the matrix a.
        dgemm_("t", "n", & s1, & s1, & a_nRows, & dV, a.givePointer(), & a_nRows, b.givePointer(), & b_nRows, & beta, this->givePointer(), & this_nRows, a_nColumns, b_nColumns, this_nColumns);
        // Second column block starting a memory position c * nRows
        dgemm_("t", "n", & this_nRows, & s2, & a_nRows, & dV, a.givePointer(), & a_nRows, & b.givePointer() [ s1 * b_nRows ], & b_nRows, & beta, & this->givePointer() [ s1 * this_nRows ], & this_nRows, a_nColumns, b_nColumns, this_nColumns);
    } else {
        // Get suitable blocksize. Around 10 rows should be suitable (slightly adjusted to minimize number of blocks):
        // Smaller blocks than ~10 didn't show any performance gains in my benchmarks. / Mikael
        int block = ( this_nRows - 1 ) / ( this_nRows / 10 ) + 1;
        int start = 0;
        int end = block;
        while ( start < this_nRows ) {
            int s = end - start;
            dgemm_("t", "n", & end, & s, & a_nRows, & dV, a.givePointer(), & a_nRows, & b.givePointer() [ start * b_nRows ], & b_nRows, & beta, & this->givePointer() [ start * this_nRows ],
                   & this_nRows, a_nColumns, b_nColumns, this_nColumns);
            start = end;
            end += block;
            if ( end > this_nRows ) {
                end = this_nRows;
            }
        }
    }
    #else
    for (Index i = 1; i <= rows(); i++ ) {
        for (Index j = i; j <= cols(); j++ ) {
            double summ = 0.;
            for (Index k = 1; k <= a.rows(); k++ ) {
                summ += a.at(k, i) * b.at(k, j);
            }
 
            this->at(i, j) += summ * dV;
        }
    }
    #endif
}
 
 
 
bool FloatMatrix :: solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose)
// solves equation b = this * x
{
#  ifndef NDEBUG
    if ( !this->isSquare() ) {
        OOFEM_ERROR("cannot solve a %d by %d matrix", rows(), cols());
    }
 
    if ( rows() != b.size() ) {
        OOFEM_ERROR("dimension mismatch");
    }
#  endif
 
#ifdef __LAPACK_MODULE
    int info, nrhs = 1;
    const int this_nRows=rows();
    IntArray ipiv(this_nRows);
    answer = b;
    //dgesv_( &this_nRows, &nrhs, this->givePointer(), &this_nRows, ipiv.givePointer(), answer.givePointer(), &this_nRows, &info );
    dgetrf_(& this_nRows, & this_nRows, this->givePointer(), & this_nRows, ipiv.givePointer(), & info);
    if ( info == 0 ) {
        dgetrs_(transpose ? "Transpose" : "No transpose", & this_nRows, & nrhs, this->givePointer(), & this_nRows, ipiv.givePointer(), answer.givePointer(), & this_nRows, & info);
    }
    if ( info != 0 ) {
        return false;
    }
#else
    Index pivRow;
    double piv, linkomb, help;
    FloatMatrix *mtrx, trans;
    if ( transpose ) {
        trans.beTranspositionOf(* this);
        mtrx = & trans;
    } else {
        mtrx = this;
    }
 
    answer = b;
 
    // initialize answer to be unity matrix;
    // lower triangle elimination by columns
    int nRows=this->rows();
    for ( Index i = 1; i < nRows; i++ ) {
        // find the suitable row and pivot
        piv = fabs( mtrx->at(i, i) );
        pivRow = i;
        for (Index j = i + 1; j <= nRows; j++ ) {
            if ( fabs( mtrx->at(j, i) ) > piv ) {
                pivRow = j;
                piv = fabs( mtrx->at(j, i) );
            }
        }
 
        if ( piv < 1.e-20 ) {
            return false;
        }
 
        // exchange rows
        if ( pivRow != i ) {
            for (Index j = i; j <= nRows; j++ ) {
                help = mtrx->at(i, j);
                mtrx->at(i, j) = mtrx->at(pivRow, j);
                mtrx->at(pivRow, j) = help;
            }
            help = answer.at(i);
            answer.at(i) = answer.at(pivRow);
            answer.at(pivRow) = help;
        }
 
        for (Index j = i + 1; j <= nRows; j++ ) {
            linkomb = mtrx->at(j, i) / mtrx->at(i, i);
            for (Index k = i; k <= nRows; k++ ) {
                mtrx->at(j, k) -= mtrx->at(i, k) * linkomb;
            }
 
            answer.at(j) -= answer.at(i) * linkomb;
        }
    }
 
    // back substitution
    for ( Index i = nRows; i >= 1; i-- ) {
        help = 0.;
        for ( Index j = i + 1; j <= nRows; j++ ) {
            help += mtrx->at(i, j) * answer.at(j);
        }
 
        answer.at(i) = ( answer.at(i) - help ) / mtrx->at(i, i);
    }
#endif
 
    return true;
}
 
 
bool FloatMatrix :: solveForRhs(const FloatMatrix &b, FloatMatrix &answer, bool transpose)
// solves equation b = this * x
// returns x. this and b are kept untouched
//
// gaussian elimination - slow but safe
//
{
#  ifndef NDEBUG
    if ( !this->isSquare() ) {
        OOFEM_ERROR("cannot solve a %d by %d matrix", rows(), cols());
    }
 
    if ( rows() != b.rows() ) {
        OOFEM_ERROR("dimension mismatch");
    }
#  endif
 
#ifdef __LAPACK_MODULE
    int info;
    const int this_nRows=rows();
    IntArray ipiv(this_nRows);
    answer = b;
    const int answer_nColumns=answer.cols();
    dgetrf_(& this_nRows, & this_nRows, this->givePointer(), & this_nRows, ipiv.givePointer(), & info);
    if ( info == 0 ) {
        dgetrs_(transpose ? "t" : "n", & this_nRows, & answer_nColumns, this->givePointer(), & this_nRows, ipiv.givePointer(), answer.givePointer(), & this_nRows, & info);
    }
    if ( info != 0 ) {
        return false; //OOFEM_ERROR("error %d", info);
    }
#else
    Index pivRow, nPs;
    double piv, linkomb, help;
    FloatMatrix *mtrx, trans;
    if ( transpose ) {
        trans.beTranspositionOf(* this);
        mtrx = & trans;
    } else {
        mtrx = this;
    }
 
    nPs = b.giveNumberOfColumns();
    answer = b;
    int nRows=this->rows();
    // initialize answer to be unity matrix;
    // lower triangle elimination by columns
    for (Index i = 1; i < nRows; i++ ) {
        // find the suitable row and pivot
        piv = fabs( mtrx->at(i, i) );
        pivRow = i;
        for (Index j = i + 1; j <= nRows; j++ ) {
            if ( fabs( mtrx->at(j, i) ) > piv ) {
                pivRow = j;
                piv = fabs( mtrx->at(j, i) );
            }
        }
 
        if ( fabs(piv) < 1.e-20 ) {
            return false; //OOFEM_ERROR("pivot too small, cannot solve %d by %d matrix", nRows, nColumns);
        }
 
        // exchange rows
        if ( pivRow != i ) {
            for (Index j = i; j <= nRows; j++ ) {
                help = mtrx->at(i, j);
                mtrx->at(i, j) = mtrx->at(pivRow, j);
                mtrx->at(pivRow, j) = help;
            }
 
            for (Index j = 1; j <= nPs; j++ ) {
                help = answer.at(i, j);
                answer.at(i, j) = answer.at(pivRow, j);
                answer.at(pivRow, j) = help;
            }
        }
 
        for (Index j = i + 1; j <= nRows; j++ ) {
            linkomb = mtrx->at(j, i) / mtrx->at(i, i);
            for (Index k = i; k <= nRows; k++ ) {
                mtrx->at(j, k) -= mtrx->at(i, k) * linkomb;
            }
 
            for (Index k = 1; k <= nPs; k++ ) {
                answer.at(j, k) -= answer.at(i, k) * linkomb;
            }
        }
    }
 
    // back substitution
    for ( Index i = nRows; i >= 1; i-- ) {
        for (Index k = 1; k <= nPs; k++ ) {
            help = 0.;
            for ( Index j = i + 1; j <= nRows; j++ ) {
                help += mtrx->at(i, j) * answer.at(j, k);
            }
 
            answer.at(i, k) = ( answer.at(i, k) - help ) / mtrx->at(i, i);
        }
    }
#endif
    return true;
}
 
 
void FloatMatrix :: printYourself() const
// Prints the receiver on screen.
{
    printf("FloatMatrix with dimensions : %ld %ld\n",
           (long)rows(), (long)cols());
    if ( rows() <= 250 && cols() <= 250 ) {
        for (Index i = 1; i <= rows(); ++i ) {
            for (Index j = 1; j <= cols() && j <= 100; ++j ) {
                printf( "%10.3e  ", this->at(i, j) );
            }
 
            printf("\n");
        }
    } else {
        printf("   large matrix : coefficients not printed \n");
    }
}
 
 
void FloatMatrix :: printYourselfToFile(const std::string filename, const bool showDimensions) const
// Prints the receiver to file.
{
    std :: ofstream matrixfile (filename);
    if (matrixfile.is_open()) {
        if (showDimensions)
            matrixfile << "FloatMatrix with dimensions : " << rows() << ", " << cols() << "\n";
        matrixfile << std::scientific << std::right << std::setprecision(3);
        for (Index i = 1; i <= rows(); ++i ) {
            for (Index j = 1; j <= cols(); ++j ) {
                matrixfile << std::setw(10) << this->at(i, j) << "\t";
            }
 
            matrixfile << "\n";
        }
        matrixfile.close();
    } else {
        OOFEM_ERROR("Failed to write to file");
    }
}
 
 
void FloatMatrix :: printYourself(const std::string &name) const
// Prints the receiver on screen.
{
    printf("%s (%ld x %ld): \n", name.c_str(), (long)rows(), (long)cols());
    if ( rows() <= 250 && cols() <= 250 ) {
        for (Index i = 1; i <= rows(); ++i ) {
            for (Index j = 1; j <= cols() && j <= 100; ++j ) {
                printf( "%10.3e  ", this->at(i, j) );
            }
 
            printf("\n");
        }
    } else {
        for (Index i = 1; i <= rows() && i <= 20; ++i ) {
            for (Index j = 1; j <= cols() && j <= 10; ++j ) {
                printf( "%10.3e  ", this->at(i, j) );
            }
            if ( cols() > 10 ) printf(" ...");
            printf("\n");
        }
        if ( rows() > 20 )  printf(" ...\n");
    }
}
 
 
void FloatMatrix :: pY() const
// Prints the receiver on screen with higher accuracy than printYourself.
{
    printf("[");
    for (Index i = 1; i <= rows(); ++i ) {
        for (Index j = 1; j <= cols(); ++j ) {
            printf( "%20.15e", this->at(i, j) );
            if ( j < cols() ) {
                printf(",");
            } else {
                printf(";");
            }
        }
    }
 
    printf("];\n");
}
 
 
void FloatMatrix :: writeCSV(const std :: string &name) const
{
    FILE *file = fopen(name.c_str(), "w");
    for (Index i = 1; i <= rows(); ++i ) {
        for (Index j = 1; j <= cols(); ++j ) {
            fprintf(file, "%10.3e, ", this->at(i, j) );
        }
 
        fprintf(file, "\n");
    }
    fclose(file);
}
 
 
void FloatMatrix :: bePinvID()
// this matrix is the product of the 6x6 deviatoric projection matrix ID
// and the inverse scaling matrix Pinv
{
    FloatMatrix& m(*this);
    m.resize(6, 6);
    m(0,0)=m(1,1)=m(2,2)=2/3.;
    m(0,1)=m(0,2)=m(1,0)=m(2,0)=m(1,2)=m(2,1)=-1/3.;
    m(3,3)=m(4,4)=m(5,5)=1/2.;
}
 
 
void FloatMatrix :: rotatedWith(const FloatMatrix &r, char mode)
// Returns the receiver 'a' rotated according the change-of-base matrix r.
// The method performs the operation  a = r^T . a . r . or the inverse
{
    FloatMatrix rta;
 
    if ( mode == 'n' ) {
        rta.beTProductOf(r, * this);     //  r^T . a
        this->beProductOf(rta, r);       //  r^T . a . r
    } else if ( mode == 't' ) {
        rta.beProductOf(r, * this);      //  r . a
        this->beProductTOf(rta, r);      //  r . a . r^T
    } else {
        OOFEM_ERROR("unsupported mode");
    }
}
 
 
 
 
 
void FloatMatrix :: beNMatrixOf(const FloatArray &n, int nsd)
{
    this->resize(nsd, n.giveSize() * nsd);
    for ( int i = 0; i < n.giveSize(); ++i ) {
        for ( int j = 0; j < nsd; ++j ) {
            ( * this )( j, i * nsd + j ) = n(i);
        }
    }
}
 
 
 
 
void FloatMatrix :: beLocalCoordSys(const FloatArray &normal)
{ //normal should be at the first position, easier for interface material models
    if ( normal.giveSize() == 1 ) {
        this->resize(1, 1);
        this->at(1, 1) = normal(0);
    } else if ( normal.giveSize() == 2 ) {
        this->resize(2, 2);
        this->at(1, 1) = normal(0);
        this->at(1, 2) = normal(1);
 
        this->at(2, 1) = normal(1);
        this->at(2, 2) = -normal(0);
 
    } else if ( normal.giveSize() == 3 ) {
        // Create a permutated vector of n, *always* length 1 and significantly different from n.
        FloatArray b, t = {
            normal(1), -normal(2), normal(0)
        };                                                    // binormal and tangent
 
        // Construct orthogonal vector
        double npn = t.dotProduct(normal);
        t.add(-npn, normal);
        t.normalize();
        b.beVectorProductOf(t, normal);
 
        this->resize(3, 3);
        this->at(1, 1) = normal.at(1);
        this->at(1, 2) = normal.at(2);
        this->at(1, 3) = normal.at(3);
 
        this->at(2, 1) = b.at(1);
        this->at(2, 2) = b.at(2);
        this->at(2, 3) = b.at(3);
 
        this->at(3, 1) = t.at(1);
        this->at(3, 2) = t.at(2);
        this->at(3, 3) = t.at(3);
    } else {
        OOFEM_ERROR("Normal needs 1 to 3 components.");
    }
}
 
 
void FloatMatrix :: beMatrixForm(const FloatArray &aArray)
{
    // Revrites the vector on matrix form (symmetrized matrix used if size is 6),
    // order: 11, 22, 33, 23, 13, 12
    // order: 11, 22, 33, 23, 13, 12, 32, 31, 21
#  ifndef NDEBUG
    if ( aArray.giveSize() != 6 && aArray.giveSize() != 9 ) {
        OOFEM_ERROR("matrix dimension is not 3x3");
    }
#  endif
    this->resize(3, 3);
    if ( aArray.giveSize() == 9 ) {
        this->at(1, 1) = aArray.at(1);
        this->at(2, 2) = aArray.at(2);
        this->at(3, 3) = aArray.at(3);
        this->at(2, 3) = aArray.at(4);
        this->at(1, 3) = aArray.at(5);
        this->at(1, 2) = aArray.at(6);
        this->at(3, 2) = aArray.at(7);
        this->at(3, 1) = aArray.at(8);
        this->at(2, 1) = aArray.at(9);
    } else if ( aArray.giveSize() == 6 ) {
        this->at(1, 1) = aArray.at(1);
        this->at(2, 2) = aArray.at(2);
        this->at(3, 3) = aArray.at(3);
        this->at(2, 3) = aArray.at(4);
        this->at(1, 3) = aArray.at(5);
        this->at(1, 2) = aArray.at(6);
        this->at(3, 2) = aArray.at(4);
        this->at(3, 1) = aArray.at(5);
        this->at(2, 1) = aArray.at(6);
    }
}
 
void FloatMatrix :: changeComponentOrder()
{
    // Changes index order between abaqus <-> OOFEM
//#  ifndef NDEBUG
//    if ( nRows != 6 || cols() != 6 ) {
//        OOFEM_ERROR("matrix dimension is not 6x6");
//    }
//#  endif
 
    if ( rows() == 6 && cols() == 6 ) {
        // This could probably be done more beautifully + efficiently.
 
        std :: swap( this->at(4, 1), this->at(6, 1) );
 
        std :: swap( this->at(4, 2), this->at(6, 2) );
 
        std :: swap( this->at(4, 3), this->at(6, 3) );
 
        std :: swap( this->at(1, 4), this->at(1, 6) );
        std :: swap( this->at(2, 4), this->at(2, 6) );
        std :: swap( this->at(3, 4), this->at(3, 6) );
        std :: swap( this->at(4, 4), this->at(6, 6) );
        std :: swap( this->at(5, 4), this->at(5, 6) );
        std :: swap( this->at(6, 4), this->at(4, 6) );
 
        std :: swap( this->at(4, 5), this->at(6, 5) );
    } else if ( rows() == 9 && cols() == 9 ) {
        // OOFEM:           11, 22, 33, 23, 13, 12, 32, 31, 21
        // UMAT:            11, 22, 33, 12, 13, 23, 32, 21, 31
        const int abq2oo [ 9 ] = {
            1,  2,  3,  6,  5,  4,  7,  9,  8
        };
 
        FloatMatrix tmp(9, 9);
        for ( int i = 1; i <= 9; i++ ) {
            for ( int j = 1; j <= 9; j++ ) {
                tmp.at(i, j) = this->at(abq2oo [ i - 1 ], abq2oo [ j - 1 ]);
            }
        }
 
        * this = tmp;
    }
}
 
 
 
double FloatMatrix :: computeReciprocalCondition(char p) const
{
#  ifndef NDEBUG
    if ( !this->isSquare() ) {
        OOFEM_ERROR("receiver must be square (is %d by %d)", this->rows(), this->cols());
    }
#  endif
    double anorm = this->computeNorm(p);
 
#  ifdef __LAPACK_MODULE
    const int n = this->rows();
    FloatArray work(4 *n);
    IntArray iwork(n);
    int info;
    double rcond;
 
    if ( n > 3 ) { // Only use this routine for larger matrices. (not sure if it's suitable for n = 3) // Mikael
        FloatMatrix a_cpy = * this;
        dgetrf_(& n, & n, a_cpy.givePointer(), & n, iwork.givePointer(), & info);
        if ( info < 0 ) {
            OOFEM_ERROR("dgetfr error %d\n", info);
        }
        dgecon_(& p, & n, a_cpy.givePointer(), & n, & anorm, & rcond, work.givePointer(), iwork.givePointer(), & info, 1);
        if ( info < 0 ) {
            OOFEM_ERROR("dgecon error %d\n", info);
        }
        return rcond;
    }
#  endif
    if ( this->giveDeterminant() <= 1e-6 * anorm ) {
        return 0.0;
    }
    FloatMatrix inv;
    inv.beInverseOf(* this);
    return 1.0 / ( inv.computeNorm(p) * anorm );
}
 
void FloatMatrix :: beMatrixFormOfStress(const FloatArray &aArray)
{
    // Revrites the  matrix on vector form (symmetrized matrix used), order: 11, 22, 33, 23, 13, 12
#  ifndef NDEBUG
    if ( aArray.giveSize() != 6 && aArray.giveSize() != 9 ) {
        OOFEM_ERROR("matrix dimension is not 3x3");
    }
#  endif
    this->resize(3, 3);
    if ( aArray.giveSize() == 9 ) {
        this->at(1, 1) = aArray.at(1);
        this->at(2, 2) = aArray.at(2);
        this->at(3, 3) = aArray.at(3);
        this->at(2, 3) = aArray.at(4);
        this->at(1, 3) = aArray.at(5);
        this->at(1, 2) = aArray.at(6);
        this->at(3, 2) = aArray.at(7);
        this->at(3, 1) = aArray.at(8);
        this->at(2, 1) = aArray.at(9);
    } else if ( aArray.giveSize() == 6 ) {
        this->at(1, 1) = aArray.at(1);
        this->at(2, 2) = aArray.at(2);
        this->at(3, 3) = aArray.at(3);
        this->at(2, 3) = aArray.at(4);
        this->at(1, 3) = aArray.at(5);
        this->at(1, 2) = aArray.at(6);
        this->at(3, 2) = aArray.at(4);
        this->at(3, 1) = aArray.at(5);
        this->at(2, 1) = aArray.at(6);
    }
}
 
#if 0
bool FloatMatrix :: computeEigenValuesSymmetric(FloatArray &lambda, FloatMatrix &v, int neigs) const
{
 #ifdef __LAPACK_MODULE
    double abstol = 1.0; 
    int lda, n, ldz, info, found, lwork;
    n = this->rows();
    lda = n;
    ldz = n;
    FloatMatrix a;
    a = * this;
    if ( neigs == 0 ) {
        neigs = n;
    }
 
    lambda.resize(neigs);
    v.resize(n, neigs);
 
    IntArray ifail(n), iwork(5 * n);
    FloatArray work( ( n + 3 ) *n ); // maximum block size should be less than n (?)
    lwork = ( n + 3 ) * n;
    if ( neigs > 0 ) {
        int one = 1;
        dsyevx_("N", "I", "U",
                & n, a.givePointer(), & n,
                NULL, NULL, & one, & neigs,
                & abstol, & found, lambda.givePointer(), v.givePointer(), & ldz,
                work.givePointer(), & lwork, iwork.givePointer(), ifail.givePointer(), & info, 1, 1, 1);
    } else {
        dsyevx_("N", "A", "U",
                & n, a.givePointer(), & n,
                NULL, NULL, NULL, NULL,
                & abstol, & found, lambda.givePointer(), v.givePointer(), & ldz,
                work.givePointer(), & lwork, iwork.givePointer(), ifail.givePointer(), & info, 1, 1, 1);
    }
 
    return info == 0;
 
 #else
    OOFEM_ERROR("Requires LAPACK");
    return false;
 
 #endif
}
#endif
 
contextIOResultType FloatMatrix :: storeYourself(DataStream &stream) const
// writes receiver's binary image into stream
// use id to distinguish some instances
// return value >0 success
//              =0 file i/o error
{
    // write size
    if ( !stream.write(rows()) ) {
        return ( CIO_IOERR );
    }
 
    if ( !stream.write(cols()) ) {
        return ( CIO_IOERR );
    }
 
    // write raw data
    if ( !stream.write(this->givePointer(), rows() * cols()) ) {
        return ( CIO_IOERR );
    }
 
    // return result back
    return CIO_OK;
}
 
 
contextIOResultType FloatMatrix :: restoreYourself(DataStream &stream)
// reads receiver from stream
// warning - overwrites existing data!
// returns 0 if file i/o error
//        -1 if id of class id is not correct
{
    int r,c;
    // read size
    if ( !stream.read(r) ) {
        return ( CIO_IOERR );
    }
 
    if ( !stream.read(c) ) {
        return ( CIO_IOERR );
    }
 
    this->resize(r,c); // FIXME: uselessly fills with zeros
 
    // read raw data
    if ( !stream.read(this->givePointer(), r*c) ) {
        return ( CIO_IOERR );
    }
 
    // return result back
    return CIO_OK;
}
 
 
int
FloatMatrix :: givePackSize(DataStream &buff) const
{
    return buff.givePackSizeOfSizet(1) + buff.givePackSizeOfSizet(1) +
           buff.givePackSizeOfDouble(rows() * cols());
}
 
#ifdef _USE_EIGEN__BROKEN
bool FloatMatrix :: jaco_(FloatArray &eval, FloatMatrix &v, int nf){
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(*this,nf);
    if(eig.info()!=Eigen::Success) return false;
    v=eig.eigenvectors();
    eval=eig.eigenvalues();
    return true;
}
#else
bool FloatMatrix :: jaco_(FloatArray &eval, FloatMatrix &v, int nf)
{
    using std::pow, std::sin, std::cos, std::atan2;
    /*
     * Solves the eigenvalues and eigenvectors of real
     * symmetric matrix by jacobi method.
     *  Written by bp. Inspired by ED WILSON jaco_ procedure.
     *
     * Parameters (input):
     * nf - number of significant figures
     *
     * Output params:
     * eval - eigen values (not sorted)
     * v    - eigenvectors (stored columvise)
     */
 
 
    /* Local variables */
    int neq = this->giveNumberOfRows();
 
    double c_b2 = .10;
    //double c_b27 = .01;
 
    /* Function Body */
#ifndef NDEBUG
    if ( !isSquare() ) {
        OOFEM_ERROR("Not square matrix");
    }
    // check for symmetry
    for ( int i = 1; i <= neq; i++ ) {
        for ( int j = i + 1; j <= neq; j++ ) {
            //if ( this->at(i, j) != this->at(j, i) ) {
            if ( fabs( this->at(i, j) - this->at(j, i) ) > 1.0e-6 ) {
                OOFEM_ERROR("Not Symmetric matrix");
            }
        }
    }
 
#endif
 
    eval.resize(neq);
    v.resize(neq, neq);
 
    for ( int i = 1; i <= neq; i++ ) {
        eval.at(i) = this->at(i, i);
    }
 
    double tol = pow(c_b2, nf);
    double sum = 0.0;
    for ( int i = 1; i <= neq; ++i ) {
        for ( int j = 1; j <= neq; ++j ) {
            sum += fabs( this->at(i, j) );
            v.at(i, j) = 0.0;
        }
 
        v.at(i, i) = 1.0;
    }
 
    if ( sum <= 0.0 ) {
        return 0;
    }
 
 
    /* ---- REDUCE MATRIX TO DIAGONAL ---------------- */
    int ite = 0;
    double ssum;
    do {
        ssum = 0.0;
        for ( int j = 2; j <= neq; ++j ) {
            int ih = j - 1;
            for ( int i = 1; i <= ih; ++i ) {
                if ( ( fabs( this->at(i, j) ) / sum ) > tol ) {
                    ssum += fabs( this->at(i, j) );
                    /* ---- CALCULATE ROTATION ANGLE ----------------- */
                    double aa = atan2( this->at(i, j) * 2.0, eval.at(i) - eval.at(j) ) /  2.0;
                    double si = sin(aa);
                    double co = cos(aa);
                    /*
                     *   // ---- MODIFY "I" AND "J" COLUMNS OF "A" AND "V"
                     *   for (k = 1; k <= neq; ++k) {
                     *    tt = this->at(k, i);
                     *    this->at(k, i) = co * tt + si * this->at(k, j);
                     *    this->at(k, j) = -si * tt + co * this->at(k, j);
                     *    tt = v.at(k, i);
                     *    v.at(k, i) = co * tt + si * v.at(k, j);
                     *    // L500:
                     *    v.at(k, j) = -si * tt + co * v.at(k, j);
                     *   }
                     *   // ---- MODIFY DIAGONAL TERMS --------------------
                     *   this->at(i, i) = co * this->at(i, i) + si * this->at(j, i);
                     *   this->at(j, j) = -si * this->at(i, j) + co * this->at(j, j);
                     *   this->at(i, j) = 0.0;
                     *   // ---- MAKE "A" MATRIX SYMMETRICAL --------------
                     *   for (k = 1; k <= neq; ++k) {
                     *    this->at(i, k) = this->at(k, i);
                     *    this->at(j, k) = this->at(k, j);
                     *    // L600:
                     *   }
                     */
                    // ---- MODIFY "I" AND "J" COLUMNS OF "A" AND "V"
                    for ( int k = 1; k < i; ++k ) {
                        double tt = this->at(k, i);
                        this->at(k, i) = co * tt + si * this->at(k, j);
                        this->at(k, j) = -si * tt + co * this->at(k, j);
                        tt = v.at(k, i);
                        v.at(k, i) = co * tt + si *v.at(k, j);
                        v.at(k, j) = -si * tt + co *v.at(k, j);
                    }
 
                    // diagonal term (i,i)
                    double tt = eval.at(i);
                    eval.at(i) = co * tt + si * this->at(i, j);
                    double aij = -si * tt + co *this->at(i, j);
                    tt = v.at(i, i);
                    v.at(i, i) = co * tt + si *v.at(i, j);
                    v.at(i, j) = -si * tt + co *v.at(i, j);
 
                    for ( int k = i + 1; k < j; ++k ) {
                        double tt = this->at(i, k);
                        this->at(i, k) = co * tt + si * this->at(k, j);
                        this->at(k, j) = -si * tt + co * this->at(k, j);
                        tt = v.at(k, i);
                        v.at(k, i) = co * tt + si *v.at(k, j);
                        v.at(k, j) = -si * tt + co *v.at(k, j);
                    }
 
                    // diagonal term (j,j)
                    tt = this->at(i, j);
                    double aji = co * tt + si *eval.at(j);
                    eval.at(j) = -si * tt + co *eval.at(j);
 
                    tt = v.at(j, i);
                    v.at(j, i) = co * tt + si *v.at(j, j);
                    v.at(j, j) = -si * tt + co *v.at(j, j);
                    //
                    for ( int k = j + 1; k <= neq; ++k ) {
                        double tt = this->at(i, k);
                        this->at(i, k) = co * tt + si *this->at(j, k);
                        this->at(j, k) = -si * tt + co *this->at(j, k);
                        tt = v.at(k, i);
                        v.at(k, i) = co * tt + si *v.at(k, j);
                        v.at(k, j) = -si * tt + co *v.at(k, j);
                    }
 
                    // ---- MODIFY DIAGONAL TERMS --------------------
                    eval.at(i) = co * eval.at(i) + si * aji;
                    eval.at(j) = -si * aij + co *eval.at(j);
                    this->at(i, j) = 0.0;
                } else {
                    /* ---- A(I,J) MADE ZERO BY ROTATION ------------- */
                    ;
                }
            }
        }
 
        /* ---- CHECK FOR CONVERGENCE -------------------- */
        if ( ++ite > 50 ) {
            OOFEM_ERROR("too many iterations");
        }
    } while ( fabs(ssum) / sum > tol );
 
    // restore original matrix
    for ( int i = 1; i <= neq; i++ ) {
        for ( int j = i; j <= neq; j++ ) {
            this->at(i, j) = this->at(j, i);
        }
    }
 
    return 0;
} /* jaco_ */
#endif
 
 
std :: ostream &operator << ( std :: ostream & out, const FloatMatrix & x )
{
    out << x.rows() << " " << x.cols() << " {";
    for (Index i = 0; i < x.rows(); ++i ) {
        for (Index j = 0; j < x.cols(); ++j ) {
            out << " " << x(i, j);
        }
        out << ";";
    }
    out << "}";
    return out;
}
 
FloatMatrix &operator *= ( FloatMatrix & x, const double & a ) {x.times(a); return x;}
FloatMatrix operator * ( const FloatMatrix &x, const double & a ) {FloatMatrix ans(x); ans.times(a); return ans;}
FloatMatrix operator * ( const double & a, const FloatMatrix &x ) {FloatMatrix ans(x); ans.times(a); return ans;}
FloatMatrix operator *( const FloatMatrix & a, const FloatMatrix & b ) {FloatMatrix ans; ans.beProductOf (a,b); return ans;}
FloatArray operator *( const FloatMatrix & a, const FloatArray & b ) {FloatArray ans; ans.beProductOf (a,b); return ans;}
FloatMatrix operator +( const FloatMatrix & a, const FloatMatrix & b ) {FloatMatrix ans(a); ans.add(b); return ans;}
FloatMatrix operator -( const FloatMatrix & a, const FloatMatrix & b ) {FloatMatrix ans(a); ans.subtract(b); return ans;}
FloatMatrix &operator += ( FloatMatrix & a, const FloatMatrix & b ) {a.add(b); return a;}
FloatMatrix &operator -= ( FloatMatrix & a, const FloatMatrix & b ) {a.subtract(b); return a;}
 
} // end namespace oofem