Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/math/full_matrix_inline_impl.hh
1e97b3b33259master

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