Homec4science
Diffusion LibMultiScale (features/new_parralelism) Browse

libmultiscale/src/math/full_matrix_inline_impl.hh
b05d5f2aba78features/new_parralelism

Last Change
Raw File

full_matrix_inline_impl.hh
View Options

/**
* @file full_matrix_inline_impl.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Jul 11 15:47:44 2014
*
* @brief This is a wrapper to a dense matrix
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale 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 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale 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 LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
extern "C" {
#include "matrix_blaswrapper.h"
// libmultiscale::UInt * factoLUFull(libmultiscale::Real * A,libmultiscale::UInt n);
// libmultiscale::UInt solveFullMatrix(libmultiscale::Real * A,libmultiscale::Real * B,libmultiscale::UInt n,libmultiscale::UInt * pivot);
// libmultiscale::UInt * Multiplication(libmultiscale::Real * A,libmultiscale::UInt m1,libmultiscale::UInt n1,libmultiscale::Real * B,libmultiscale::UInt m2,libmultiscale::UInt n2,libmultiscale::Real * C);
// libmultiscale::UInt inverseFull(libmultiscale::Real * A, libmultiscale::UInt * pivot,libmultiscale::UInt n);
}
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
namespace math{
inline void Matrix::transpose(){
std::vector<Real> temp(_m*_n);
for (UInt i = 0 ; i < _n ; ++i)
for (UInt j = 0 ; j < _m ; ++j){
temp[i+_n*j] = matrix[j+_m*i];
}
UInt sav = _n;
_n = _m;
_m = sav;
matrix = temp;
}
/* -------------------------------------------------------------------------- */
// return inverseFull(&matrix[0], &pivot[0],_m);
//int inverseFull(double * A, int * pivot,int n){
inline void Matrix::inverse(){
int ok;
//int ispec = 1;
int ldwork = 100*_m;
std::vector<Real> work(ldwork);
int n = _m;
dgetri_(&n,&matrix[0],&n,&pivot[0],&work[0],&ldwork,&ok);
}
/* -------------------------------------------------------------------------- */
inline void Matrix::factoLU() {
if (_m!=_n){LM_FATAL("FactoLU possible :) que pour matrice carree");}
int n = _m;
double moy_err;
int ok,c=1,i;
char trans = 'N';
double unf = 1.0f,zero= 0.0f;
std::vector<Real> B(n);
B.assign(n,1.);
std::vector<Real> sav_A(n*n);
sav_A = matrix;
Real * A = &matrix[0];
std::vector<Real> C(n);
C.assign(n,0.);
pivot.resize(n);
pivot.assign(n,0.);
dgesv_(&n, &c, A, &n, &pivot[0], &B[0], &n,&ok);
for(i=0;i<n;++i)
{
B[i]=1.0f;
}
dgetrs_(&trans,&n,&c, A,&n,&pivot[0],&B[0],&n,&ok);
dgemv_(&trans, &n, &n,&unf, &sav_A[0], &n, &B[0], &c,&zero,&C[0], &c );
for(moy_err=0,i=0;i<n;++i)
{
moy_err +=fabs(1.0f-C[i]);
}
moy_err /= n;
if (moy_err > 1e-12)
LM_FATAL("FactoLU resoud avec erreur a " << moy_err);
}
/* -------------------------------------------------------------------------- */
inline void Matrix::solve(Array * r){
if (_m != _n){
LM_FATAL("Solve for square matrix");
}
int ok,c=1;
char trans = 'N';
int n = _m;
Real * B = &r->getArray()[0];
dgetrs_(&trans,&n,&c, &matrix[0],&n,&pivot[0],B,&n,&ok);
}
/* -------------------------------------------------------------------------- */
template <typename _Array>
inline void Matrix::rightMultiply(_Array & r, _Array & output){
if (_n != r.size()){
LM_FATAL("multiplying vector size does not match " << _n << " " << r.size());
}
if (_m != output.size()){
LM_FATAL("result vector size does not match " << _m << " " << output.size());
}
for (UInt i=0;i < _m ; ++i) {
output[i] = 0.;
for (UInt j=0;j < _n ; ++j)
output[i] += r[j]*(*this)(i,j);
}
}
/* -------------------------------------------------------------------------- */
inline void Matrix::zero(){
matrix.assign(_m*_n,0.);
}
/* -------------------------------------------------------------------------- */
// Multiplication(f.matrix,f.m(),f.n(),g.matrix,g.m(),g.n(),this->matrix);
//int * Multiplication(double * A,int m1,int n1,double * B,int m2,int n2,double * C)
inline void Matrix::matProduct(math::Matrix & f, math::Matrix & g){
char trans = 'N';
int n1 = f.n();
int n2 = g.n();
int m2 = g.m();
int m1 = f.m();
int k = n1;
Real unf = 1.0;
Real zero = 0.;
if (n1 != m2)exit(-1);
dgemm_(&trans,&trans, &m1, &n2,&k,&unf,
&f.getArray()[0], &m1,
&g.getArray()[0], &k,
&zero,&this->matrix[0], &m1);
}
/* -------------------------------------------------------------------------- */
}
__END_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
c4science · Help