Commit c04103bc authored by Benoît LÉ's avatar Benoît LÉ

[xTLS] Added printMatrixMarket function to coefWarperSparse

May be useful to debug.
parent 4e1d4427
/*
/*
This file is a part of eXlibris C++ Library
under the GNU Lesser General Public License.
See the NOTICE.md & LICENSE.md files for terms
See the NOTICE.md & LICENSE.md files for terms
and conditions.
*/
......@@ -9,150 +9,157 @@
#include "TLSFastMarchingTools.h"
// solverBase
#include "xGraphMatrix.h"
#include "xGenericSparseMatrix.h"
#include "xGraphMatrix.h"
namespace xtls
{
// ------------------------------------------------------------------------------------------------------
// uniteMetric implementation
double uniteMetric(const AOMD::mVertex &v){ return 1.; }
double uniteMetric(const AOMD::mVertex &v) { return 1.; }
// ------------------------------------------------------------------------------------------------------
// default function to filter column: take all
bool filterMatrixColumn(const AOMD::mEntity *e) {return true; }
bool filterMatrixColumn(const AOMD::mEntity *e) { return true; }
// ------------------------------------------------------------------------------------------------------
// coefWarperSparse implementation
coefWarperSparse::coefWarperSparse() :
coef_matrix(nullptr),graph(nullptr),n(0),m(0),ptr(nullptr),rowidx(nullptr),dat(nullptr) {}
coefWarperSparse::coefWarperSparse()
: coef_matrix(nullptr), graph(nullptr), n(0), m(0), ptr(nullptr), rowidx(nullptr), dat(nullptr)
{
}
coefWarperSparse::~coefWarperSparse()
{
if (coef_matrix) delete coef_matrix;
if (graph) delete graph;
if (coef_matrix) delete coef_matrix;
if (graph) delete graph;
}
void coefWarperSparse::reset(int n_,int m_)
void coefWarperSparse::reset(int n_, int m_)
{
if (graph) delete graph;
graph = nullptr;
n = n_;
m = m_;
graph = new xlinalg::xGraphMatrix(n,m);
if (graph) delete graph;
graph = nullptr;
n = n_;
m = m_;
graph = new xlinalg::xGraphMatrix(n, m);
}
void coefWarperSparse::addGraph(int i, int j)
{
assert(graph);
graph->add(i,j);
assert(graph);
graph->add(i, j);
}
void coefWarperSparse::finalizeGraph()
{
assert(graph);
graph->countNNZ();
if (coef_matrix) delete coef_matrix;
coef_matrix = nullptr;
coef_matrix = new coef_t(*graph);
int nc,mc,nnz;
xTraitMatrixSparceCSC storage;
xTraitMatrixCindex indexing;
coef_matrix->getMemoryAccess(nc,mc,nnz,&ptr,&rowidx,&dat,storage,indexing);
assert(ptr);
assert(n==nc);
assert(m==mc);
assert(rowidx);
graph->clear();
assert(graph);
graph->countNNZ();
if (coef_matrix) delete coef_matrix;
coef_matrix = nullptr;
coef_matrix = new coef_t(*graph);
int nc, mc, nnz;
xTraitMatrixSparceCSC storage;
xTraitMatrixCindex indexing;
coef_matrix->getMemoryAccess(nc, mc, nnz, &ptr, &rowidx, &dat, storage, indexing);
assert(ptr);
assert(n == nc);
assert(m == mc);
assert(rowidx);
graph->clear();
}
void coefWarperSparse::addMatrix(int i, int j, double val)
{
assert(coef_matrix);
int *pbeg = &rowidx[ptr[j]];
if (i < *pbeg)
throw -1;
int *pend = &rowidx[ptr[j+1]];
int *pos = std::lower_bound(pbeg, pend,i);
if (pos != pend && *pos == i)
dat[std::distance(rowidx, pos)] += val;
else
throw -1;
assert(coef_matrix);
int *pbeg = &rowidx[ptr[j]];
if (i < *pbeg) throw -1;
int *pend = &rowidx[ptr[j + 1]];
int *pos = std::lower_bound(pbeg, pend, i);
if (pos != pend && *pos == i)
dat[std::distance(rowidx, pos)] += val;
else
throw -1;
}
int coefWarperSparse::nline() const
{
assert(coef_matrix);
return coef_matrix->getN();
assert(coef_matrix);
return coef_matrix->getN();
}
int coefWarperSparse::ncolumn() const
{
assert(coef_matrix);
return coef_matrix->getM();
assert(coef_matrix);
return coef_matrix->getM();
}
double coefWarperSparse::operator ()( int i, int j ) const
double coefWarperSparse::operator()(int i, int j) const
{
assert(coef_matrix);
int *pbeg = &rowidx[ptr[j]];
if (i < *pbeg)
return 0.;
int *pend = &rowidx[ptr[j+1]];
int *pos = std::lower_bound(pbeg, pend,i);
if (pos != pend && *pos == i)
return dat[std::distance(rowidx, pos)];
else
return 0.;
assert(coef_matrix);
int *pbeg = &rowidx[ptr[j]];
if (i < *pbeg) return 0.;
int *pend = &rowidx[ptr[j + 1]];
int *pos = std::lower_bound(pbeg, pend, i);
if (pos != pend && *pos == i)
return dat[std::distance(rowidx, pos)];
else
return 0.;
}
void coefWarperSparse::gemv(int size_X,int size_Y, double alpha, double beta,double *X, double*Y,bool sparse_X, bool transpose )
void coefWarperSparse::gemv(int size_X, int size_Y, double alpha, double beta, double *X, double *Y, bool sparse_X,
bool transpose)
{
assert(coef_matrix);
coef_matrix->gemv(size_X,size_Y,alpha,beta,X,Y,sparse_X,transpose);
assert(coef_matrix);
coef_matrix->gemv(size_X, size_Y, alpha, beta, X, Y, sparse_X, transpose);
}
void coefWarperSparse::gemm(char * trans_B,int nbr_OPB, int nbc_OPB,int nbr_C, int ldb, int ldc, double alpha, double beta,double * B, double * C)
void coefWarperSparse::gemm(char *trans_B, int nbr_OPB, int nbc_OPB, int nbr_C, int ldb, int ldc, double alpha, double beta,
double *B, double *C)
{
coef_matrix->gemm(trans_B,nbr_OPB,nbc_OPB,nbr_C,ldb,ldc,alpha,beta,B,C);
coef_matrix->gemm(trans_B, nbr_OPB, nbc_OPB, nbr_C, ldb, ldc, alpha, beta, B, C);
}
std::unique_ptr < int[] > coefWarperSparse::giveNonNullRow()
std::unique_ptr<int[]> coefWarperSparse::giveNonNullRow()
{
assert(coef_matrix);
std::unique_ptr < int[] > non_null (new int[n]);
for (int idx = 0; idx < n; ++idx) non_null[idx] = 0;
assert(coef_matrix);
std::unique_ptr<int[]> non_null(new int[n]);
for (int idx = 0; idx < n; ++idx) non_null[idx] = 0;
for (int idx = 0,idxe = ptr[m]; idx < idxe; ++idx)
{
++non_null[rowidx[idx]];
}
return non_null;
for (int idx = 0, idxe = ptr[m]; idx < idxe; ++idx)
{
++non_null[rowidx[idx]];
}
return non_null;
}
void coefWarperSparse::mergeRow(int *sel_n)
{
assert(sel_n);
assert(coef_matrix);
coef_t *tmp;
tmp = new coef_t(*coef_matrix,sel_n,nullptr,true);
delete coef_matrix;
coef_matrix = tmp;
int nnz;
xTraitMatrixSparceCSC storage;
xTraitMatrixCindex indexing;
coef_matrix->getMemoryAccess(n,m,nnz,&ptr,&rowidx,&dat,storage,indexing);
assert(ptr);
assert(rowidx);
}
void coefWarperSparse::eliminateRow(P_t &P)
{
xlinalg::xGraphMatrix g_tmp(P.getN(),m);
coef_t *tmp;
tmp = new coef_t(g_tmp);
assert(coef_matrix);
P.gemm(1.,0.,*coef_matrix,*tmp);
delete coef_matrix;
coef_matrix = tmp;
int nnz;
xTraitMatrixSparceCSC storage;
xTraitMatrixCindex indexing;
coef_matrix->getMemoryAccess(n,m,nnz,&ptr,&rowidx,&dat,storage,indexing);
assert(ptr);
assert(rowidx);
}
void coefWarperSparse::toDense(int nbr, int nbc, double * dense_matrix, int* ids_r, int *ids_c, std::function < int (int,int) > denseIndex)
{
assert(coef_matrix);
coef_matrix->toDense(nbr,nbc,dense_matrix,ids_r,ids_c,denseIndex);
}
} //end namespace xtls
assert(sel_n);
assert(coef_matrix);
coef_t *tmp;
tmp = new coef_t(*coef_matrix, sel_n, nullptr, true);
delete coef_matrix;
coef_matrix = tmp;
int nnz;
xTraitMatrixSparceCSC storage;
xTraitMatrixCindex indexing;
coef_matrix->getMemoryAccess(n, m, nnz, &ptr, &rowidx, &dat, storage, indexing);
assert(ptr);
assert(rowidx);
}
void coefWarperSparse::eliminateRow(P_t &P)
{
xlinalg::xGraphMatrix g_tmp(P.getN(), m);
coef_t *tmp;
tmp = new coef_t(g_tmp);
assert(coef_matrix);
P.gemm(1., 0., *coef_matrix, *tmp);
delete coef_matrix;
coef_matrix = tmp;
int nnz;
xTraitMatrixSparceCSC storage;
xTraitMatrixCindex indexing;
coef_matrix->getMemoryAccess(n, m, nnz, &ptr, &rowidx, &dat, storage, indexing);
assert(ptr);
assert(rowidx);
}
void coefWarperSparse::toDense(int nbr, int nbc, double *dense_matrix, int *ids_r, int *ids_c,
std::function<int(int, int)> denseIndex)
{
assert(coef_matrix);
coef_matrix->toDense(nbr, nbc, dense_matrix, ids_r, ids_c, denseIndex);
}
void coefWarperSparse::printMatrixMarket(std::string filename)
{
std::ofstream oss(filename);
coef_matrix->printMatrixMarket(oss);
}
} // end namespace xtls
/*
/*
This file is a part of eXlibris C++ Library
under the GNU Lesser General Public License.
See the NOTICE.md & LICENSE.md files for terms
See the NOTICE.md & LICENSE.md files for terms
and conditions.
*/
#ifndef TLSFASTMARCHINGTOOLS_H
#define TLSFASTMARCHINGTOOLS_H
//std
#include <iostream>
#include <fstream>
// std
#include <algorithm>
#include <cassert>
#include <fstream>
#include <functional>
#include <iostream>
#include <vector>
#include <cassert>
// eXlibris_tools
#include "xIteratorTools.h"
#include "xDataExchanger.h"
#include "xIteratorTools.h"
// xfem
#include "xRegion.h"
#include "xLevelSet.h"
#include "xRegion.h"
//xlinalg
#include "xTraitsMatrix.h"
// xlinalg
#include "xSparseVectorPackUnPack.h"
#include "xTraitsMatrix.h"
#ifdef HAVE_FASTMARCHING
// Fast Marching
#include "meshinterfacexRegion.h"
#include "FM.h"
template < typename T >
using datamanagerFM_t = xinterface::aomd::xAttachedDataManagerAOMDGeneric < AOMD::mVertex, T > ;
#include "meshinterfacexRegion.h"
template <typename T>
using datamanagerFM_t = xinterface::aomd::xAttachedDataManagerAOMDGeneric<AOMD::mVertex, T>;
#ifdef HAVE_XGRAPH
// xgraph
#include "nodeAndConnectedEdgeDist.h"
#include "nodeAndConnectedEdge.h"
#include "nodeAndConnectedEdgeDist.h"
#endif
#endif
namespace xlinalg
{
template < typename T, typename DEFINED, typename PATTERN, typename STORAGE_TYPE, typename INDEXING >
template <typename T, typename DEFINED, typename PATTERN, typename STORAGE_TYPE, typename INDEXING>
class xGenericSparseMatrix;
class xGraphMatrix;
}
} // namespace xlinalg
namespace xtls
{
/// function used as default metric for FM
double uniteMetric(const AOMD::mVertex &v);
......@@ -60,86 +59,74 @@ bool filterMatrixColumn(const AOMD::mEntity *e);
/// warper class to store coefficient of modal matrix in a standalone container that fit to requirement of TLSAlgo for now
class coefWarperSparse
{
public:
typedef xlinalg::xGenericSparseMatrix < double, xTraitMatrixNonSingular, xTraitMatrixUnSym, xTraitMatrixSparceCSC, xTraitMatrixCindex > P_t;
coefWarperSparse();
~coefWarperSparse();
void reset(int n,int m);
void addGraph(int i, int j);
void finalizeGraph();
void addMatrix(int i, int j, double val);
int nline() const;
int ncolumn() const;
double operator () (int i, int j) const;
void gemv(int size_X,int size_Y, double alpha, double beta,double *X, double*Y,bool sparse_X = false, bool transpose = false);
void gemm(char * trans_B,int nbr_OPB, int nbc_OPB,int nbr_C, int ldb, int ldc, double alpha, double beta,double * B, double * C);
std::unique_ptr < int[] > giveNonNullRow();
void mergeRow(int *sel_n);
void eliminateRow(P_t & P);
void toDense(int nbr,int nbc, double * dense_matrix,int* ids_r,int *ids_c,std::function < int (int,int) > denseIndex);
private:
typedef xlinalg::xGenericSparseMatrix < double, xTraitMatrixNonSingular, xTraitMatrixUnSym, xTraitMatrixSparceCSC, xTraitMatrixCindex > coef_t;
coef_t *coef_matrix;
xlinalg::xGraphMatrix * graph;
int n,m;
int *ptr;
int *rowidx;
double *dat;
public:
typedef xlinalg::xGenericSparseMatrix<double, xTraitMatrixNonSingular, xTraitMatrixUnSym, xTraitMatrixSparceCSC,
xTraitMatrixCindex>
P_t;
coefWarperSparse();
~coefWarperSparse();
void reset(int n, int m);
void addGraph(int i, int j);
void finalizeGraph();
void addMatrix(int i, int j, double val);
int nline() const;
int ncolumn() const;
double operator()(int i, int j) const;
void gemv(int size_X, int size_Y, double alpha, double beta, double *X, double *Y, bool sparse_X = false,
bool transpose = false);
void gemm(char *trans_B, int nbr_OPB, int nbc_OPB, int nbr_C, int ldb, int ldc, double alpha, double beta, double *B,
double *C);
std::unique_ptr<int[]> giveNonNullRow();
void mergeRow(int *sel_n);
void eliminateRow(P_t &P);
void toDense(int nbr, int nbc, double *dense_matrix, int *ids_r, int *ids_c, std::function<int(int, int)> denseIndex);
void printMatrixMarket(std::string filename);
private:
typedef xlinalg::xGenericSparseMatrix<double, xTraitMatrixNonSingular, xTraitMatrixUnSym, xTraitMatrixSparceCSC,
xTraitMatrixCindex>
coef_t;
coef_t *coef_matrix;
xlinalg::xGraphMatrix *graph;
int n, m;
int *ptr;
int *rowidx;
double *dat;
};
#ifdef HAVE_FASTMARCHING
/// Compute updated level set values, gradiant and coeficient matrix of modes
template < typename ITERKNOWN,
typename ITERKNOWNV,
typename ITERTOFIND,
typename ITERNEWKNOWN,
typename ITERMATNOD,
typename MATRIX >
void TLSFastMarchingValGradMode(std::function < double (AOMD::mVertex *) > get_val,
std::function < void (AOMD::mVertex *, double) > set_val,
std::function < void (AOMD::mVertex *, const xtensor::xVector<> &, double) > set_grad,
const xfem::xRegion & region,
xtool::xRange < ITERKNOWN > known,
xtool::xRange < ITERKNOWNV > known_volatil,
xtool::xRange < ITERTOFIND > to_find,
ITERNEWKNOWN new_known,
ITERMATNOD mat_node_it,
MATRIX * coeffs,
double epsilon_ratio,
const std::function < bool (const AOMD::mEntity *) > & filter_func = filterMatrixColumn,
const std::function < double (const AOMD::mVertex &) > & f_func = uniteMetric);
} //end namespace xtls
template <typename ITERKNOWN, typename ITERKNOWNV, typename ITERTOFIND, typename ITERNEWKNOWN, typename ITERMATNOD,
typename MATRIX>
void TLSFastMarchingValGradMode(std::function<double(AOMD::mVertex *)> get_val,
std::function<void(AOMD::mVertex *, double)> set_val,
std::function<void(AOMD::mVertex *, const xtensor::xVector<> &, double)> set_grad,
const xfem::xRegion &region, xtool::xRange<ITERKNOWN> known,
xtool::xRange<ITERKNOWNV> known_volatil, xtool::xRange<ITERTOFIND> to_find,
ITERNEWKNOWN new_known, ITERMATNOD mat_node_it, MATRIX *coeffs, double epsilon_ratio,
const std::function<bool(const AOMD::mEntity *)> &filter_func = filterMatrixColumn,
const std::function<double(const AOMD::mVertex &)> &f_func = uniteMetric);
} // end namespace xtls
#include "TLSFastMarchingTools_imp.h"
#else
// fake
template < typename ITERKNOWN,
typename ITERKNOWNV,
typename ITERTOFIND,
typename ITERNEWKNOWN,
typename ITERMATNOD,
typename MATRIX >
void TLSFastMarchingValGradMode(std::function < double (AOMD::mVertex *) > get_val,
std::function < void (AOMD::mVertex *, double) > set_val,
std::function < void (AOMD::mVertex *, const xtensor::xVector<> &, double) > set_grad,
const xfem::xRegion & region,
xtool::xRange < ITERKNOWN > known,
xtool::xRange < ITERKNOWNV > known_volatil,
xtool::xRange < ITERTOFIND > to_find,
ITERNEWKNOWN new_known,
ITERMATNOD mat_node_it,
MATRIX * coeffs,
double epsilon_ratio,
const std::function < bool (const AOMD::mEntity *) > & filter_func = filterMatrixColumn,
const std::function < double (const AOMD::mVertex &) > & f_func = uniteMetric)
template <typename ITERKNOWN, typename ITERKNOWNV, typename ITERTOFIND, typename ITERNEWKNOWN, typename ITERMATNOD,
typename MATRIX>
void TLSFastMarchingValGradMode(std::function<double(AOMD::mVertex *)> get_val,
std::function<void(AOMD::mVertex *, double)> set_val,
std::function<void(AOMD::mVertex *, const xtensor::xVector<> &, double)> set_grad,
const xfem::xRegion &region, xtool::xRange<ITERKNOWN> known,
xtool::xRange<ITERKNOWNV> known_volatil, xtool::xRange<ITERTOFIND> to_find,
ITERNEWKNOWN new_known, ITERMATNOD mat_node_it, MATRIX *coeffs, double epsilon_ratio,
const std::function<bool(const AOMD::mEntity *)> &filter_func = filterMatrixColumn,
const std::function<double(const AOMD::mVertex &)> &f_func = uniteMetric)
{
throw -1;
throw -1;
}
} //end namespace xtls
} // end namespace xtls
#endif
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment