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!" );
269 if(
size() !=
m.size()) {
270 throw std::range_error(
"AlSpaMat::operator+: size do not match!" );
276 b.m_ptr_map[
pos->first] +=
pos->second;
277 b.m_nele =
b.m_ptr_map.size();
285 if(
size() !=
m.size()) {
286 throw std::range_error(
"AlSpaMat::operator+=: size do not match!" );
290 for (
pos =
m.m_ptr_map.begin();
pos!=
m.m_ptr_map.end(); ++
pos)
291 (*this).m_ptr_map[
pos->first] +=
pos->second;
300 if(
size() !=
m.size()) {
301 throw std::range_error(
"AlSpaMat::operator-: size do not match!" );
307 b.m_ptr_map[
pos->first] -=
pos->second;
308 b.m_nele =
b.m_ptr_map.size();
316 if(
size() !=
m.size()) {
317 throw std::range_error(
"AlSpaMat::operator-=: size do not match!" );
321 for (
pos =
m.m_ptr_map.begin();
pos!=
m.m_ptr_map.end(); ++
pos)
322 (*this).m_ptr_map[
pos->first] -=
pos->second;
331 if(
size() !=
m.size() ) {
332 throw std::range_error(
"AlSpaMat::operator*: size do not match!" );
335 long int isiz(
size());
338 for(
long int i=0;
i<isiz;
i++ )
339 for(
long int j=0; j<isiz; j++ )
340 for(
int k=0;
k<isiz;
k++)
348 if(
size() !=
m.nrow() ) {
349 throw std::range_error(
"AlSpaMat::operator*: size do not match!" );
352 long int isiz(
size());
353 long int icol(
m.ncol());
356 for(
long int i=0;
i<isiz;
i++ )
357 for(
long int j=0; j<
icol; j++ )
358 for(
long int k=0;
k<isiz;
k++)
367 if(
size() !=
v.size() ) {
368 throw std::range_error(
"AlSpaMat::operator*: size do not match! " );
371 long int isiz(
size());
374 for(
long int i=0;
i<isiz;
i++ )
375 for(
long int j=0; j<isiz; j++ )
397 a.m_ptr_map.insert(std::make_pair(
pos->first, (
pos->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);
428 tripletList.emplace_back(
i,j,
pos->second);
429 if(
i!=j) tripletList.emplace_back(j,
i,
pos->second);
431 eigenBigMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
439 Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
441 solver.compute(eigenBigMatrix);
442 if(solver.info()!=Eigen::Success) {
444 throw std::domain_error(
"AlSpaMat::SolveWithEigen: failed to compute: -- is the input matrix singular?" );
448 Eigen::VectorXd
x = solver.solve( eigenBigVector );
449 if(solver.info()!=Eigen::Success) {
451 throw std::domain_error(
"AlSpaMat::SolveWithEigen: Failed to solve: -- is your matrix singular? ");
456 for(
int i(0);
i< RHS.
size(); ++
i ){
461 Eigen::VectorXd
residual = eigenBigMatrix *
x - eigenBigVector;
462 double sumresidual = 0;
468 if( sumresidual > 1
e-3 ){
469 throw std::overflow_error(
"AlSpaMat::SolveWithEigen: your solution is no good! ");
484 throw std::invalid_argument(
"AlSpaMat::diagonalize: not implemented!" );
490 throw std::invalid_argument(
"AlSpaMat::invert: not implemented!" );
539 throw std::invalid_argument(
"AlSpaMat::RemoveCollsRows: Vector of indices larger than matrix size." );
545 for (
int i=0;
i<
n-1;
i++)
546 for (
int j=
i+1; j<
n; j++)
554 for (
int i=0;
i<
n;
i++) {
556 throw std::out_of_range(
"AlSpaMat::RemoveCollsRows: Index goes beyond matrix." );
601 if (counterCol==5-control) {
608 if (counterRow==5-control) {
638 for (
pos =
m.m_ptr_map.begin();
pos!=
m.m_ptr_map.end(); ++
pos) {
639 m.elem(
pos->first,
i, j);
666 std::ofstream outmat;
668 int32_t msizz = 1000000+
size();
674 return StatusCode::FAILURE;
677 outmat.write((
char*)&msizz,
sizeof (msizz));
679 outmat.write((
char*)&nelem,
sizeof (nelem));
684 return StatusCode::FAILURE;
685 outmat.setf(std::ios::fixed);
686 outmat.setf(std::ios::showpoint);
688 outmat <<
"msizz: " << std::setw(6) << msizz << std::endl;
689 outmat <<
"nelem: " << std::setw(6) << nelem << std::endl;
690 outmat <<
"AlSpaMat version: " << std::setw(6) <<
version << std::endl;
700 elem(
pos->first, ii, jj);
i=ii; j=jj;
702 outmat.write((
char*)&(
i),
sizeof (
i));
703 outmat.write((
char*)&(j),
sizeof (j));
704 outmat.write((
char*)&(melem),
sizeof (melem));
707 outmat << std::setw(6) <<
i << std::setw(6) << j << std::setw(18) << melem << std::endl;
711 return StatusCode::SUCCESS;
717 std::ifstream inmat((
filename).c_str(), std::ios::binary);
719 return StatusCode::FAILURE;
724 inmat.read((
char*)&msiz,
sizeof (msiz));
733 return StatusCode::SUCCESS;
739 bool stdUnits =
true;
741 return StatusCode::FAILURE;
745 return StatusCode::FAILURE;
751 inmat.read((
char*)&msiz,
sizeof (msiz));
766 inmat.read((
char*)&nelem,
sizeof (nelem));
768 for(
int k=0;
k<nelem;
k++) {
769 inmat.read((
char*)&
i,
sizeof (
i));
770 inmat.read((
char*)&j,
sizeof (j));
771 inmat.read((
char*)&melem,
sizeof (melem));
772 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
777 for(int32_t
i=0;
i<msiz;
i++) {
778 for(int32_t j=0; j<msiz; j++) {
779 inmat.read((
char*)&melem,
sizeof (melem));
780 if(
i>=j && melem!=0. )
781 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
788 for( int32_t
i=0;
i<msiz;
i++) {
789 for( int32_t j=0; j<=
i; j++) {
790 inmat.read((
char*)&melem,
sizeof (melem));
792 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
799 return StatusCode::SUCCESS;
808 return StatusCode::FAILURE;
814 inmat.read((
char*)&msiz,
sizeof (msiz));
828 inmat.read((
char*)&nelem,
sizeof (nelem));
830 for(
int k=0;
k<nelem;
k++) {
831 inmat.read((
char*)&
i,
sizeof (
i));
832 inmat.read((
char*)&j,
sizeof (j));
833 inmat.read((
char*)&melem,
sizeof (melem));
834 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
839 for(int32_t
i=0;
i<msiz;
i++) {
840 for(int32_t j=0; j<msiz; j++) {
841 inmat.read((
char*)&melem,
sizeof (melem));
842 if(
i>=j && melem!=0. )
843 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
850 for( int32_t
i=0;
i<msiz;
i++) {
851 for( int32_t j=0; j<=
i; j++) {
852 inmat.read((
char*)&melem,
sizeof (melem));
854 m_ptr_map.insert(std::make_pair(std::make_pair(
i,j), melem));
861 return StatusCode::SUCCESS;
867 long int nonZeroElements =
m_ptr_map.size();
868 int *irow =
new int[nonZeroElements];
869 int *
icol =
new int[nonZeroElements];
870 double *
val =
new double[nonZeroElements];
876 i =
pos->first.first;
877 j =
pos->first.second;
885 myTMatrix->SetMatrixArray(nonZeroElements,irow,
icol,
val);