9 #include "GaudiKernel/StatusCode.h"
25 #include <TMatrixDSparse.h>
29 #include <Eigen/IterativeLinearSolvers>
30 #include <Eigen/SparseCholesky>
87 if(
size() !=
m.size()) {
88 throw std::range_error(
"AlSpaMat::copy: size do not match!" );
98 if(
size() !=
m.size()) {
99 throw std::range_error(
"AlSpaMat::copy: size do not match!" );
106 for (
int j = 0; j<=
i; j++)
107 if (
m.elemc(
i,j) != 0.) {
108 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j),
m.elemc(
i,j)));
117 if(
size() !=
m.nrow() ||
size() !=
m.ncol() ) {
118 throw std::range_error(
"AlSpaMat::copy: size do not match!" );
127 for (
int j = 0; j<=
i; j++)
128 if (
m.elemc(
i,j) != 0.) {
129 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j),
m.elemc(
i,j)));
140 throw std::out_of_range(
"AlSpaMat::elemr: Index 1 < zero! " );
143 throw std::out_of_range(
"AlSpaMat::elemr: Index 1 too large! " );
146 throw std::out_of_range(
"AlSpaMat::elemr: Index 2 < zero! " );
149 throw std::out_of_range(
"AlSpaMat::elemr: Index 2 too large! " );;
154 indices key = j <
i ? std::make_pair(
i,j) : std::make_pair(j,
i);
161 return (
m_ptr_map.insert(std::make_pair(
key, 0.))).first->second;
176 throw std::out_of_range(
"AlSpaMat::elemc: Index 1 < zero! " );
179 throw std::out_of_range(
"AlSpaMat::elemc: Index 1 too large! " );
182 throw std::out_of_range(
"AlSpaMat::elemc: Index 2 < zero! " );
185 throw std::out_of_range(
"AlSpaMat::elemc: Index 2 too large! " );
190 indices key = j <
i ? std::make_pair(
i,j) : std::make_pair(j,
i);
207 throw std::out_of_range(
"AlSpaMat::elem: Index 1 < zero! " );
210 throw std::out_of_range(
"AlSpaMat::elem: Index 1 too large! " );
213 throw std::out_of_range(
"AlSpaMat::elem: Index 2 < zero! " );
216 throw std::out_of_range(
"AlSpaMat::elem: Index 2 too large! " );
220 indices key = j <
i ? std::make_pair(
i,j) : std::make_pair(j,
i);
247 if(
m.nrow() !=
m.ncol() ) {
248 throw std::range_error(
"AlSpaMat::=operator: allowed for square matrices only!" );
268 if(
size() !=
m.size()) {
269 throw std::range_error(
"AlSpaMat::operator+: size do not match!" );
273 for (
const datamap::value_type&
p :
m_ptr_map)
274 b.m_ptr_map[
p.first] +=
p.second;
275 b.m_nele =
b.m_ptr_map.size();
283 if(
size() !=
m.size()) {
284 throw std::range_error(
"AlSpaMat::operator+=: size do not match!" );
287 for (
const datamap::value_type&
p :
m.m_ptr_map)
288 (*this).m_ptr_map[
p.first] +=
p.second;
297 if(
size() !=
m.size()) {
298 throw std::range_error(
"AlSpaMat::operator-: size do not match!" );
302 for (
const datamap::value_type&
p :
m_ptr_map)
303 b.m_ptr_map[
p.first] -=
p.second;
304 b.m_nele =
b.m_ptr_map.size();
312 if(
size() !=
m.size()) {
313 throw std::range_error(
"AlSpaMat::operator-=: size do not match!" );
316 for (
const datamap::value_type&
p :
m.m_ptr_map)
317 (*this).m_ptr_map[
p.first] -=
p.second;
326 if(
size() !=
m.size() ) {
327 throw std::range_error(
"AlSpaMat::operator*: size do not match!" );
330 long int isiz(
size());
333 for(
long int i=0;
i<isiz;
i++ )
334 for(
long int j=0; j<isiz; j++ )
335 for(
int k=0;
k<isiz;
k++)
343 if(
size() !=
m.nrow() ) {
344 throw std::range_error(
"AlSpaMat::operator*: size do not match!" );
347 long int isiz(
size());
348 long int icol(
m.ncol());
351 for(
long int i=0;
i<isiz;
i++ )
352 for(
long int j=0; j<
icol; j++ )
353 for(
long int k=0;
k<isiz;
k++)
362 if(
size() !=
v.size() ) {
363 throw std::range_error(
"AlSpaMat::operator*: size do not match! " );
366 long int isiz(
size());
369 for(
long int i=0;
i<isiz;
i++ )
370 for(
long int j=0; j<isiz; j++ )
389 for (
const datamap::value_type&
p :
m_ptr_map)
390 a.m_ptr_map.emplace(
p.first,
p.second*
d);
398 throw std::invalid_argument(
"AlSpaMat::determinant: not implemented!" );
404 throw std::range_error(
"AlSpaMat::SolveWithEigen vector size is incorrect" );
407 Eigen::VectorXd eigenBigVector( RHS.
size() );
408 for(
int i(0);
i< RHS.
size(); ++
i ){
409 eigenBigVector[
i] = RHS[
i];
412 Eigen::SparseMatrix<double> eigenBigMatrix(
m_size,
m_size );
414 using Triplet = Eigen::Triplet<double>;
415 std::vector<Triplet> tripletList;
416 tripletList.reserve(
m_nele);
418 for (
const datamap::value_type&
p :
m_ptr_map) {
420 tripletList.emplace_back(
i,j,
p.second);
421 if(
i!=j) tripletList.emplace_back(j,
i,
p.second);
423 eigenBigMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
431 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
433 solver.compute(eigenBigMatrix);
434 if(solver.info()!=Eigen::Success) {
436 throw std::domain_error(
"AlSpaMat::SolveWithEigen: failed to compute: -- is the input matrix singular?" );
440 Eigen::VectorXd
x = solver.solve( eigenBigVector );
441 if(solver.info()!=Eigen::Success) {
443 throw std::domain_error(
"AlSpaMat::SolveWithEigen: Failed to solve: -- is your matrix singular? ");
448 for(
int i(0);
i< RHS.
size(); ++
i ){
453 Eigen::VectorXd
residual = eigenBigMatrix *
x - eigenBigVector;
454 double sumresidual = 0;
460 if( sumresidual > 1
e-3 ){
461 throw std::overflow_error(
"AlSpaMat::SolveWithEigen: your solution is no good! ");
476 throw std::invalid_argument(
"AlSpaMat::diagonalize: not implemented!" );
482 throw std::invalid_argument(
"AlSpaMat::invert: not implemented!" );
531 throw std::invalid_argument(
"AlSpaMat::RemoveCollsRows: Vector of indices larger than matrix size." );
537 for (
int i=0;
i<
n-1;
i++)
538 for (
int j=
i+1; j<
n; j++)
546 for (
int i=0;
i<
n;
i++) {
548 throw std::out_of_range(
"AlSpaMat::RemoveCollsRows: Index goes beyond matrix." );
593 if (counterCol==5-control) {
600 if (counterRow==5-control) {
629 for (
const datamap::value_type&
p :
m.m_ptr_map) {
630 m.elem(
p.first,
i, j);
657 std::ofstream outmat;
659 int32_t msizz = 1000000+
size();
665 return StatusCode::FAILURE;
668 outmat.write((
char*)&msizz,
sizeof (msizz));
670 outmat.write((
char*)&nelem,
sizeof (nelem));
675 return StatusCode::FAILURE;
676 outmat.setf(std::ios::fixed);
677 outmat.setf(std::ios::showpoint);
679 outmat <<
"msizz: " << std::setw(6) << msizz << std::endl;
680 outmat <<
"nelem: " << std::setw(6) << nelem << std::endl;
681 outmat <<
"AlSpaMat version: " << std::setw(6) <<
version << std::endl;
688 for (
const datamap::value_type&
p :
m_ptr_map) {
690 elem(
p.first, ii, jj);
i=ii; j=jj;
692 outmat.write((
char*)&(
i),
sizeof (
i));
693 outmat.write((
char*)&(j),
sizeof (j));
694 outmat.write((
char*)&(melem),
sizeof (melem));
697 outmat << std::setw(6) <<
i << std::setw(6) << j << std::setw(18) << melem << std::endl;
701 return StatusCode::SUCCESS;
707 std::ifstream inmat((
filename).c_str(), std::ios::binary);
709 return StatusCode::FAILURE;
714 inmat.read((
char*)&msiz,
sizeof (msiz));
723 return StatusCode::SUCCESS;
729 bool stdUnits =
true;
731 return StatusCode::FAILURE;
735 return StatusCode::FAILURE;
741 inmat.read((
char*)&msiz,
sizeof (msiz));
756 inmat.read((
char*)&nelem,
sizeof (nelem));
758 for(
int k=0;
k<nelem;
k++) {
759 inmat.read((
char*)&
i,
sizeof (
i));
760 inmat.read((
char*)&j,
sizeof (j));
761 inmat.read((
char*)&melem,
sizeof (melem));
762 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
767 for(int32_t
i=0;
i<msiz;
i++) {
768 for(int32_t j=0; j<msiz; j++) {
769 inmat.read((
char*)&melem,
sizeof (melem));
770 if(
i>=j && melem!=0. )
771 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
778 for( int32_t
i=0;
i<msiz;
i++) {
779 for( int32_t j=0; j<=
i; j++) {
780 inmat.read((
char*)&melem,
sizeof (melem));
782 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
789 return StatusCode::SUCCESS;
798 return StatusCode::FAILURE;
804 inmat.read((
char*)&msiz,
sizeof (msiz));
818 inmat.read((
char*)&nelem,
sizeof (nelem));
820 for(
int k=0;
k<nelem;
k++) {
821 inmat.read((
char*)&
i,
sizeof (
i));
822 inmat.read((
char*)&j,
sizeof (j));
823 inmat.read((
char*)&melem,
sizeof (melem));
824 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
829 for(int32_t
i=0;
i<msiz;
i++) {
830 for(int32_t j=0; j<msiz; j++) {
831 inmat.read((
char*)&melem,
sizeof (melem));
832 if(
i>=j && melem!=0. )
833 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
840 for( int32_t
i=0;
i<msiz;
i++) {
841 for( int32_t j=0; j<=
i; j++) {
842 inmat.read((
char*)&melem,
sizeof (melem));
844 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
851 return StatusCode::SUCCESS;
857 long int nonZeroElements =
m_ptr_map.size();
858 int *irow =
new int[nonZeroElements];
859 int *
icol =
new int[nonZeroElements];
860 double *
val =
new double[nonZeroElements];
864 for (
const datamap::value_type&
p :
m_ptr_map) {
874 myTMatrix->SetMatrixArray(nonZeroElements,irow,
icol,
val);