9 #include "GaudiKernel/StatusCode.h"
25 #include <TMatrixDSparse.h>
29 #include <Eigen/IterativeLinearSolvers>
30 #include <Eigen/SparseCholesky>
34 charAddress(
auto &
v){
35 return reinterpret_cast<char *
>(&
v);
94 if(
size() !=
m.size()) {
95 throw std::range_error(
"AlSpaMat::copy: size do not match!" );
105 if(
size() !=
m.size()) {
106 throw std::range_error(
"AlSpaMat::copy: size do not match!" );
113 for (
int j = 0; j<=
i; j++)
114 if (
m.elemc(
i,j) != 0.) {
115 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j),
m.elemc(
i,j)));
124 if(
size() !=
m.nrow() ||
size() !=
m.ncol() ) {
125 throw std::range_error(
"AlSpaMat::copy: size do not match!" );
134 for (
int j = 0; j<=
i; j++)
135 if (
m.elemc(
i,j) != 0.) {
136 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j),
m.elemc(
i,j)));
147 throw std::out_of_range(
"AlSpaMat::elemr: Index 1 < zero! " );
150 throw std::out_of_range(
"AlSpaMat::elemr: Index 1 too large! " );
153 throw std::out_of_range(
"AlSpaMat::elemr: Index 2 < zero! " );
156 throw std::out_of_range(
"AlSpaMat::elemr: Index 2 too large! " );;
161 indices key = j <
i ? std::make_pair(
i,j) : std::make_pair(j,
i);
168 return (
m_ptr_map.insert(std::make_pair(
key, 0.))).first->second;
183 throw std::out_of_range(
"AlSpaMat::elemc: Index 1 < zero! " );
186 throw std::out_of_range(
"AlSpaMat::elemc: Index 1 too large! " );
189 throw std::out_of_range(
"AlSpaMat::elemc: Index 2 < zero! " );
192 throw std::out_of_range(
"AlSpaMat::elemc: Index 2 too large! " );
197 indices key = j <
i ? std::make_pair(
i,j) : std::make_pair(j,
i);
214 throw std::out_of_range(
"AlSpaMat::elem: Index 1 < zero! " );
217 throw std::out_of_range(
"AlSpaMat::elem: Index 1 too large! " );
220 throw std::out_of_range(
"AlSpaMat::elem: Index 2 < zero! " );
223 throw std::out_of_range(
"AlSpaMat::elem: Index 2 too large! " );
227 indices key = j <
i ? std::make_pair(
i,j) : std::make_pair(j,
i);
254 if(
m.nrow() !=
m.ncol() ) {
255 throw std::range_error(
"AlSpaMat::=operator: allowed for square matrices only!" );
275 if(
size() !=
m.size()) {
276 throw std::range_error(
"AlSpaMat::operator+: size do not match!" );
280 for (
const datamap::value_type&
p :
m_ptr_map)
281 b.m_ptr_map[
p.first] +=
p.second;
282 b.m_nele =
b.m_ptr_map.size();
290 if(
size() !=
m.size()) {
291 throw std::range_error(
"AlSpaMat::operator+=: size do not match!" );
294 for (
const datamap::value_type&
p :
m.m_ptr_map)
295 (*this).m_ptr_map[
p.first] +=
p.second;
304 if(
size() !=
m.size()) {
305 throw std::range_error(
"AlSpaMat::operator-: size do not match!" );
309 for (
const datamap::value_type&
p :
m_ptr_map)
310 b.m_ptr_map[
p.first] -=
p.second;
311 b.m_nele =
b.m_ptr_map.size();
319 if(
size() !=
m.size()) {
320 throw std::range_error(
"AlSpaMat::operator-=: size do not match!" );
323 for (
const datamap::value_type&
p :
m.m_ptr_map)
324 (*this).m_ptr_map[
p.first] -=
p.second;
333 if(
size() !=
m.size() ) {
334 throw std::range_error(
"AlSpaMat::operator*: size do not match!" );
337 long int isiz(
size());
340 for(
long int i=0;
i<isiz;
i++ )
341 for(
long int j=0; j<isiz; j++ )
342 for(
int k=0;
k<isiz;
k++)
350 if(
size() !=
m.nrow() ) {
351 throw std::range_error(
"AlSpaMat::operator*: size do not match!" );
354 long int isiz(
size());
355 long int icol(
m.ncol());
358 for(
long int i=0;
i<isiz;
i++ )
359 for(
long int j=0; j<
icol; j++ )
360 for(
long int k=0;
k<isiz;
k++)
369 if(
size() !=
v.size() ) {
370 throw std::range_error(
"AlSpaMat::operator*: size do not match! " );
373 long int isiz(
size());
376 for(
long int i=0;
i<isiz;
i++ )
377 for(
long int j=0; j<isiz; j++ )
396 for (
const datamap::value_type&
p :
m_ptr_map)
397 a.m_ptr_map.emplace(
p.first,
p.second*
d);
405 throw std::invalid_argument(
"AlSpaMat::determinant: not implemented!" );
411 throw std::range_error(
"AlSpaMat::SolveWithEigen vector size is incorrect" );
414 Eigen::VectorXd eigenBigVector( RHS.
size() );
415 for(
int i(0);
i< RHS.
size(); ++
i ){
416 eigenBigVector[
i] = RHS[
i];
419 Eigen::SparseMatrix<double> eigenBigMatrix(
m_size,
m_size );
421 using Triplet = Eigen::Triplet<double>;
422 std::vector<Triplet> tripletList;
423 tripletList.reserve(
m_nele);
425 for (
const datamap::value_type&
p :
m_ptr_map) {
427 tripletList.emplace_back(
i,j,
p.second);
428 if(
i!=j) tripletList.emplace_back(j,
i,
p.second);
430 eigenBigMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
438 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
440 solver.compute(eigenBigMatrix);
441 if(solver.info()!=Eigen::Success) {
443 throw std::domain_error(
"AlSpaMat::SolveWithEigen: failed to compute: -- is the input matrix singular?" );
447 Eigen::VectorXd
x = solver.solve( eigenBigVector );
448 if(solver.info()!=Eigen::Success) {
450 throw std::domain_error(
"AlSpaMat::SolveWithEigen: Failed to solve: -- is your matrix singular? ");
455 for(
int i(0);
i< RHS.
size(); ++
i ){
460 Eigen::VectorXd
residual = eigenBigMatrix *
x - eigenBigVector;
461 double sumresidual = 0;
467 if( sumresidual > 1
e-3 ){
468 throw std::overflow_error(
"AlSpaMat::SolveWithEigen: your solution is no good! ");
483 throw std::invalid_argument(
"AlSpaMat::diagonalize: not implemented!" );
489 throw std::invalid_argument(
"AlSpaMat::invert: not implemented!" );
538 throw std::invalid_argument(
"AlSpaMat::RemoveCollsRows: Vector of indices larger than matrix size." );
544 for (
int i=0;
i<
n-1;
i++)
545 for (
int j=
i+1; j<
n; j++)
553 for (
int i=0;
i<
n;
i++) {
555 throw std::out_of_range(
"AlSpaMat::RemoveCollsRows: Index goes beyond matrix." );
582 for(
int col=0; col<
row+1; col++) {
600 if (counterCol==5-control) {
607 if (counterRow==5-control) {
621 const int ind=shift*
index;
636 for (
const datamap::value_type&
p :
m.m_ptr_map) {
637 m.elem(
p.first,
i, j);
664 std::ofstream outmat;
666 int32_t msizz = 1000000+
size();
672 return StatusCode::FAILURE;
675 outmat.write(charAddress(msizz),
sizeof (msizz));
677 outmat.write(charAddress(nelem),
sizeof (nelem));
682 return StatusCode::FAILURE;
683 outmat.setf(std::ios::fixed);
684 outmat.setf(std::ios::showpoint);
686 outmat <<
"msizz: " << std::setw(6) << msizz << std::endl;
687 outmat <<
"nelem: " << std::setw(6) << nelem << std::endl;
688 outmat <<
"AlSpaMat version: " << std::setw(6) <<
version << std::endl;
695 for (
const datamap::value_type&
p :
m_ptr_map) {
697 elem(
p.first, ii, jj);
i=ii; j=jj;
699 outmat.write(charAddress((
i)),
sizeof (
i));
700 outmat.write(charAddress((j)),
sizeof (j));
701 outmat.write(charAddress((melem)),
sizeof (melem));
704 outmat << std::setw(6) <<
i << std::setw(6) << j << std::setw(18) << melem << std::endl;
708 return StatusCode::SUCCESS;
714 std::ifstream inmat((
filename).c_str(), std::ios::binary);
716 return StatusCode::FAILURE;
721 inmat.read(charAddress(msiz),
sizeof (msiz));
730 return StatusCode::SUCCESS;
736 bool stdUnits =
true;
738 return StatusCode::FAILURE;
742 return StatusCode::FAILURE;
748 inmat.read(charAddress(msiz),
sizeof (msiz));
763 inmat.read(charAddress(nelem),
sizeof (nelem));
765 for(
int k=0;
k<nelem;
k++) {
766 inmat.read(charAddress(
i),
sizeof (
i));
767 inmat.read(charAddress(j),
sizeof (j));
768 inmat.read(charAddress(melem),
sizeof (melem));
769 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
774 for(int32_t
i=0;
i<msiz;
i++) {
775 for(int32_t j=0; j<msiz; j++) {
776 inmat.read(charAddress(melem),
sizeof (melem));
777 if(
i>=j && melem!=0. )
778 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
785 for( int32_t
i=0;
i<msiz;
i++) {
786 for( int32_t j=0; j<=
i; j++) {
787 inmat.read(charAddress(melem),
sizeof (melem));
789 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
796 return StatusCode::SUCCESS;
805 return StatusCode::FAILURE;
811 inmat.read(charAddress(msiz),
sizeof (msiz));
815 dofs = std::abs(msiz);
825 inmat.read(charAddress(nelem),
sizeof (nelem));
827 for(
int k=0;
k<nelem;
k++) {
828 inmat.read(charAddress(
i),
sizeof (
i));
829 inmat.read(charAddress(j),
sizeof (j));
830 inmat.read(charAddress(melem),
sizeof (melem));
831 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
836 for(int32_t
i=0;
i<msiz;
i++) {
837 for(int32_t j=0; j<msiz; j++) {
838 inmat.read(charAddress(melem),
sizeof (melem));
839 if(
i>=j && melem!=0. )
840 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
847 for( int32_t
i=0;
i<msiz;
i++) {
848 for( int32_t j=0; j<=
i; j++) {
849 inmat.read(charAddress(melem),
sizeof (melem));
851 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
858 return StatusCode::SUCCESS;
864 long int nonZeroElements =
m_ptr_map.size();
865 int *irow =
new int[nonZeroElements];
866 int *
icol =
new int[nonZeroElements];
867 double *
val =
new double[nonZeroElements];
871 for (
const datamap::value_type&
p :
m_ptr_map) {
881 myTMatrix->SetMatrixArray(nonZeroElements,irow,
icol,
val);