25 #include <Base/Common/W32Compat.hh>
27 #include "SparseMatrix.hh"
28 #include <MathAlgo/SVD.hh>
33 #define ABS(x) (((x) < 0) ? -(x) : (x))
34 #define MAX(x,y) ((x) > (y)) ? (x) : (y)
47 for (
unsigned int i = 0; i < numRow_; i++) {
48 std::list<std::pair<unsigned int, double> > row;
49 for (
unsigned int j = 0; j < numCol_; j++) {
51 if (ABS(val) > SPARSE_MATRIX_EPSILON) {
52 row.push_back(std::pair<unsigned int, double>(j, val));
53 maxRow_ = MAX(maxRow_, i);
54 maxCol_ = MAX(maxCol_, j);
57 if (!row.empty()) rows_[i] = row;
101 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
102 const_iterator rowsIter = rows_.find(row);
103 if (rowsIter != rows_.end()) {
104 std::list<std::pair<unsigned int, double> >::const_iterator colsIter;
105 std::list<std::pair<unsigned int, double> >::const_iterator colsEnd
106 = rowsIter->second.end();
107 for (colsIter = rowsIter->second.begin();
108 colsIter != colsEnd && colsIter->first < col; colsIter++);
109 if (colsIter != colsEnd && colsIter->first == col) {
110 return colsIter->second;
119 maxRow_ = MAX(maxRow_, row);
120 maxCol_ = MAX(maxCol_, col);
121 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
122 iterator rowsIter = rows_.find(row);
123 if (rowsIter == rows_.end() && ABS(val) <= SPARSE_MATRIX_EPSILON) {
126 std::list<std::pair<unsigned int, double> >::iterator it;
127 std::list<std::pair<unsigned int, double> >::iterator end = rows_[row].end();
128 for (it = rows_[row].begin(); it != end && it->first < col; it++);
129 if (it != end && it->first == col) {
130 if (ABS(val) > SPARSE_MATRIX_EPSILON) {
133 rows_[row].erase(it);
134 if (rows_[row].empty()) rows_.erase(row);
136 }
else if (ABS(val) > SPARSE_MATRIX_EPSILON) {
137 rows_[row].insert(it, std::pair<unsigned int,double>(col, val));
140 if (numRow_ < (row + 1) || numCol_ < (col + 1)){
141 BIASWARNONCE(
"Element was set out of current matrix dimensions. "
142 <<
"Dimensions have been increased!")
143 numRow_ = MAX(numRow_, (row + 1));
144 numCol_ = MAX(numCol_, (col + 1));
152 unsigned int startRowInA,
unsigned int startColInA,
153 unsigned int row_count,
unsigned int col_count)
155 BIASASSERT(startRowInA+row_count <= (
unsigned int)(A.
GetRows()));
156 BIASASSERT(startColInA+col_count <= (
unsigned int)(A.
GetCols()));
162 maxRow_ = MAX(maxRow_, row+row_count-1);
163 maxCol_ = MAX(maxCol_, col+col_count-1);
164 for (
unsigned int i = 0; i < row_count; i++) {
165 std::list<std::pair<unsigned int, double> >::iterator it;
166 std::list<std::pair<unsigned int, double> >::iterator end
167 = rows_[row+i].end();
168 for (it = rows_[row+i].begin(); it != end && it->first < col; it++);
171 BIASASSERT(it == end || it->first >= col+col_count);
173 for (
unsigned int j = 0; j < col_count; j++) {
174 if (ABS(A[startRowInA+i][startColInA+j]) > SPARSE_MATRIX_EPSILON) {
175 rows_[row+i].insert(it, std::pair<unsigned int, double>
176 (col+j, A[startRowInA+i][startColInA+j]));
181 if (numRow_ < (row+row_count) || numCol_ < (col+col_count)){
182 BIASWARNONCE(
"Sub matrix exceeds current matrix dimensions. "
183 <<
"Dimensions have been increased!");
184 numRow_ = MAX(numRow_, (row+row_count));
185 numCol_ = MAX(numCol_, (col+col_count));
196 result.
rows_.clear();
197 std::map< unsigned int, std::list<std::pair<unsigned int, double> > >::
199 std::map< unsigned int, std::list<std::pair<unsigned int, double> > >::
200 iterator rowsEnd = rows_.end();
201 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
202 std::list<std::pair<unsigned int, double> >::iterator colsIter;
203 std::list<std::pair<unsigned int, double> >::iterator colsEnd
204 = rowsIter->second.end();
205 for (colsIter = rowsIter->second.begin();
206 colsIter != colsEnd; colsIter++) {
209 result.
rows_[colsIter->first].push_back(
210 std::pair<unsigned int, double>(rowsIter->first, colsIter->second));
223 std::list<std::pair<unsigned int, double> >::iterator it;
224 std::list<std::pair<unsigned int, double> >::iterator end = rows_[row].end();
225 for (it = rows_[row].begin(); it != end && it->first < col; it++);
228 BIASASSERT(it == end || it->first > col);
230 rows_[row].insert(it, std::pair<unsigned int, double>(col, value));
236 std::map<
unsigned int,
237 std::list<std::pair<unsigned int, double> > >::iterator row;
238 std::list< std::pair<unsigned int, double> >::iterator elem;
239 for (row = rows_.begin(); row != rows_.end(); row++) {
240 for (elem = row->second.begin(); elem != row->second.end(); elem++) {
241 elem->second *= multiplier;
253 MultiplyWithTransposeOf(AT, result);
261 BIASERR(
"does not work for some reason, returing with doing anything...");
268 result.
rows_.clear();
269 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
271 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
272 iterator rowsEnd1 = rows_.end();
273 for (rowsIter1 = rows_.begin(); rowsIter1 != rowsEnd1; rowsIter1++) {
274 std::list<std::pair<unsigned int, double> > row;
275 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
277 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
278 iterator rowsEnd2 = A.
rows_.end();
279 for (rowsIter2 = A.
rows_.begin(); rowsIter2 != rowsEnd2; rowsIter2++) {
281 std::list<std::pair<unsigned int, double> >::iterator colsIter1 =
282 rowsIter1->second.begin();
283 std::list<std::pair<unsigned int, double> >::iterator colsIter2 =
284 rowsIter2->second.begin();
285 std::list<std::pair<unsigned int, double> >::iterator colsEnd1 =
286 rowsIter1->second.end();
287 std::list<std::pair<unsigned int, double> >::iterator colsEnd2 =
288 rowsIter2->second.end();
289 while (colsIter1 != colsEnd1 && colsIter2 != colsEnd2) {
290 if (colsIter1->first < colsIter2->first) {
292 }
else if (colsIter1->first > colsIter2->first) {
295 sum += colsIter1->second * colsIter2->second;
300 if (ABS(sum) > SPARSE_MATRIX_EPSILON) {
303 row.push_back(std::pair<unsigned int, double>(rowsIter2->first, sum));
306 if (!row.empty()) result.
rows_[rowsIter1->first] = row;
318 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
319 const_iterator rowsIter;
320 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
321 const_iterator rowsEnd = rows_.end();
322 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
323 BIASASSERT(rowsIter->first < M.
GetRows());
324 std::list<std::pair<unsigned int, double> >::const_iterator colsIter;
325 std::list<std::pair<unsigned int, double> >::const_iterator colsEnd
326 = rowsIter->second.end();
327 for (colsIter = rowsIter->second.begin();
328 colsIter != colsEnd; colsIter++) {
329 BIASASSERT(colsIter->first < M.
GetCols());
330 M[rowsIter->first][colsIter->first] = colsIter->second;
339 double maxValue = 0.0;
340 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
341 const_iterator rowsIter;
342 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
343 const_iterator rowsEnd = rows_.end();
344 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
345 std::list<std::pair<unsigned int, double> >::const_iterator colsIter;
346 std::list<std::pair<unsigned int, double> >::const_iterator colsEnd
347 = rowsIter->second.end();
348 for (colsIter = rowsIter->second.begin();
349 colsIter != colsEnd; colsIter++) {
351 if (colsIter->first > rowsIter->first)
break;
352 if (colsIter->first == rowsIter->first) {
353 if (colsIter->second > maxValue) {
354 maxValue = colsIter->second;
355 index = colsIter->first;
367 double maxValue = 0.0;
368 map<unsigned int, std::list<std::pair<unsigned int, double> > >::
369 const_iterator rowsIter = rows_.find(row);
370 if (rowsIter != rows_.end()) {
371 std::list<std::pair<unsigned int, double> >::const_iterator colsIter;
372 std::list<std::pair<unsigned int, double> >::const_iterator colsEnd
373 = rowsIter->second.end();
374 for (colsIter = rowsIter->second.begin();
375 colsIter != colsEnd; colsIter++) {
376 if (colsIter->second > maxValue) {
377 maxValue = colsIter->second;
378 col = colsIter->first;
390 double maxValue = 0.0;
391 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
392 const_iterator rowsIter;
393 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
394 const_iterator rowsEnd = rows_.end();
395 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
396 std::list<std::pair<unsigned int, double> >::const_iterator colsIter;
397 std::list<std::pair<unsigned int, double> >::const_iterator colsEnd
398 = rowsIter->second.end();
399 for (colsIter = rowsIter->second.begin();
400 colsIter != colsEnd; colsIter++) {
402 if (colsIter->first > col)
break;
403 if (colsIter->first == col) {
404 if (colsIter->second > maxValue) {
405 maxValue = colsIter->second;
406 row = rowsIter->first;
416 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
418 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
419 iterator rowsEnd = rows_.end();
421 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
422 std::list<std::pair<unsigned int, double> >::iterator colsIter;
423 std::list<std::pair<unsigned int, double> >::iterator colsEnd
424 = rowsIter->second.end();
426 for (colsIter = rowsIter->second.begin();
427 colsIter != colsEnd; colsIter++) {
429 if (colsIter->first > rowsIter->first)
break;
431 if (colsIter->first == rowsIter->first) {
432 colsIter->second *= value;
444 BIASASSERT(maxRow_ == maxCol_);
445 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
446 iterator rowsIter = rows_.begin();
449 for (
unsigned int rowId = 0; rowId <= maxRow_; rowId++) {
451 if (rowsIter == rows_.end() || rowsIter->first != rowId) {
452 std::list<pair<unsigned int, double> > newRow;
453 newRow.push_back(std::pair<unsigned int, double>(rowId, value));
454 rowsIter = rows_.insert(rowsIter,
455 std::pair<
unsigned int,
456 list<pair<unsigned int, double> > >(rowId, newRow));
460 std::list<std::pair<unsigned int, double> >::iterator colsIter;
461 std::list<std::pair<unsigned int, double> >::iterator colsEnd
462 = rowsIter->second.end();
463 bool missingCol =
true;
464 for (colsIter = rowsIter->second.begin();
465 colsIter != colsEnd && missingCol; colsIter++) {
467 if (colsIter->first > rowId)
break;
469 if (colsIter->first == rowId) {
470 colsIter->second += value;
477 rowsIter->second.insert(colsIter,
478 std::pair<unsigned int, double>(rowId, value));
485 if (numRow_ < (maxRow_+1) || numCol_ < (maxCol_+1)){
486 BIASWARNONCE(
"Sub matrix exceeds current matrix dimensions. "
487 <<
"Dimensions have been increased!");
488 numRow_ = MAX(numRow_, (maxRow_+1));
489 numCol_ = MAX(numCol_, (maxCol_+1));
495 unsigned int numRows,
unsigned int numCols,
500 map<unsigned int, std::list<std::pair<unsigned int, double> > >::
501 const_iterator rowsIter = rows_.find(row);
502 map<unsigned int, std::list<std::pair<unsigned int, double> > >::
503 const_iterator rowsEnd = rows_.end();
506 unsigned int rowOffset = 0;
507 while (rowsIter == rowsEnd && rowOffset < numRows) {
508 rowsIter = rows_.find(row + (++rowOffset));
516 for (; rowsIter != rowsEnd; rowsIter++) {
518 if (rowsIter->first >= numRows+row)
break;
519 std::list<std::pair<unsigned int, double> >::const_iterator colsIter;
520 std::list<std::pair<unsigned int, double> >::const_iterator colsEnd
521 = rowsIter->second.end();
523 for (colsIter = rowsIter->second.begin();
524 colsIter != colsEnd; colsIter++) {
526 if (colsIter == colsEnd) {
527 BIASERR(
"Sparse matrix is not consistent (contains empty rows)!");
531 if (colsIter == colsEnd)
break;
534 if (colsIter->first >= col)
break;
536 for (; colsIter != colsEnd; colsIter++) {
537 if (colsIter->first >= numCols+col)
break;
538 BIASASSERT(rowsIter->first >= row && rowsIter->first < M.
GetRows()+row);
539 BIASASSERT(colsIter->first >= col && colsIter->first < M.
GetCols()+col);
540 M[rowsIter->first-row][colsIter->first-col] = colsIter->second;
546 Print(std::string
const &name, std::ostream &stream)
const
548 stream << name <<
" m " << GetRowsNum() <<
"n " << GetColsNum() << endl;
549 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
550 const_iterator rowsIter;
551 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
552 const_iterator rowsEnd = rows_.end();
553 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
554 std::list<std::pair<unsigned int, double> >::const_iterator colsIter;
555 std::list<std::pair<unsigned int, double> >::const_iterator colsEnd
556 = rowsIter->second.end();
557 for (colsIter = rowsIter->second.begin();
558 colsIter != colsEnd; colsIter++) {
559 stream << name <<
"[" << rowsIter->first <<
"][" << colsIter->first
560 <<
"] = " << colsIter->second << endl;
572 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
574 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
575 iterator rowsEnd = rows_.end();
576 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
578 std::list<std::pair<unsigned int, double> >::iterator colsIter;
579 std::list<std::pair<unsigned int, double> >::iterator colsEnd
580 = rowsIter->second.end();
581 for (colsIter = rowsIter->second.begin();
582 colsIter != colsEnd; colsIter++) {
583 sum += colsIter->second * x[colsIter->first];
585 result[rowsIter->first] = sum;
592 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
594 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
595 iterator rowsEnd = H.
rows_.end();
598 for (rowsIter = H.
rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
599 rows_[numRow_+rowsIter->first] = rowsIter->second;
610 maxCol_ = MAX(maxCol_, H.
maxCol_);
631 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
633 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
634 iterator rowsEnd = H.
rows_.end();
635 for (rowsIter = H.
rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
636 std::list<std::pair<unsigned int, double> > row = rows_[rowsIter->first];
637 std::list<std::pair<unsigned int, double> >::iterator colsIter;
638 std::list<std::pair<unsigned int, double> >::iterator colsEnd
639 = rowsIter->second.end();
640 for (colsIter = rowsIter->second.begin();
641 colsIter != colsEnd; colsIter++) {
642 row.push_back(std::pair<unsigned int,double>(maxCol_+colsIter->first +1,
645 rows_[rowsIter->first] = row;
654 maxRow_ = MAX(maxRow_, H.
maxRow_);
675 BIASASSERT ((numRow_ == numCol_) && (b.
Size() == numRow_));
678 for (
unsigned int i = 0; i <= maxRow_; i++) {
680 rows_[i].push_back(std::pair<unsigned int, double>(maxCol_+i+1, 1.0));
681 if (ABS(b[i]) > SPARSE_MATRIX_EPSILON) {
683 rows_[i].push_back(std::pair<unsigned int, double>(2*maxCol_+2, b[i]));
688 int res = GaussJordan_();
689 if (res != 0)
return res;
692 inverse.
newsize(maxRow_+1, maxCol_+1);
696 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
698 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
699 iterator rowsEnd = rows_.end();
700 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
701 std::list<std::pair<unsigned int, double> >::iterator colsIter;
702 std::list<std::pair<unsigned int, double> >::iterator colsEnd
703 = rowsIter->second.end();
704 for (colsIter = rowsIter->second.begin();
705 colsIter != colsEnd && colsIter->first <= maxCol_; colsIter++);
707 if (colsIter == colsEnd) {
708 BIASERR(
"Sparse matrix is not consistent (contains empty rows)!");
711 for (; colsIter != colsEnd; colsIter++) {
712 if (colsIter->first < 2*maxCol_+2) {
713 inverse[rowsIter->first][colsIter->first-maxCol_-1] = colsIter->second;
715 result[rowsIter->first] = colsIter->second;
721 maxCol_ = 2*maxCol_+2;
723 if ((numCol_ < maxCol_ + 1)){
724 BIASWARNONCE(
"Call to 'InvertAndSolve' enlarged current matrix which now exceeds "
725 <<
"its dimensions. Dimensions have been increased!");
726 numCol_ = maxCol_ + 1;
738 BIASASSERT ((numRow_ == numCol_) && (b.
Size() == numRow_));
741 for (
unsigned int i = 0; i <= maxRow_; i++) {
743 rows_[i].push_back(std::pair<unsigned int, double>(maxCol_+i+1, 1));
744 if (ABS(b[i]) > SPARSE_MATRIX_EPSILON) {
746 rows_[i].push_back(std::pair<unsigned int,double>(2*maxCol_+2, b[i]));
751 int res = GaussJordan_();
752 if (res != 0)
return res;
757 inverse.
rows_.clear();
760 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
762 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
763 iterator rowsEnd = rows_.end();
764 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
765 std::list<std::pair<unsigned int, double> >::iterator colsIter;
766 std::list<std::pair<unsigned int, double> >::iterator colsEnd
767 = rowsIter->second.end();
768 for (colsIter = rowsIter->second.begin();
769 colsIter != colsEnd && colsIter->first <= maxCol_; colsIter++);
771 if (colsIter == colsEnd) {
772 BIASERR(
"Sparse matrix is not consistent (contains empty rows)!");
775 std::list<std::pair<unsigned int, double> > row;
776 for (; colsIter != colsEnd; colsIter++) {
777 if (colsIter->first < 2*maxCol_+2) {
778 row.push_back(std::pair<unsigned int, double>(colsIter->first-maxCol_-1,
781 result[rowsIter->first] = colsIter->second;
784 if (!row.empty()) inverse.
rows_[rowsIter->first] = row;
788 maxCol_ = maxCol_*2+2;
791 if ((numCol_ < maxCol_ + 1) || (numRow_ < maxRow_ + 1)){
792 BIASWARNONCE(
"Call to 'InvertAndSolve' enlarged current matrix which now exceeds "
793 <<
"its dimensions. Dimensions have been increased!");
794 numCol_ = maxCol_ + 1;
795 numRow_ = maxRow_ + 1;
806 BIASASSERT ((numRow_ == numCol_) && (b.
Size() == numRow_));
810 for (
unsigned int i = 0; i <= maxRow_; i++) {
811 if (ABS(b[i]) > SPARSE_MATRIX_EPSILON) {
813 rows_[i].push_back(std::pair<unsigned int, double>(maxCol_+1, b[i]));
818 int res = GaussJordan_();
819 if (res != 0)
return res;
824 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
826 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
827 iterator rowsEnd = rows_.end();
828 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
829 std::pair<unsigned int, double> back = rowsIter->second.back();
830 if (back.first == maxCol_+1) {
831 result[rowsIter->first] = back.second;
839 if (((numCol_ < maxCol_ + 1) || (numRow_ < maxRow_ + 1))){
840 BIASWARNONCE(
"Call to 'Invert' enlarged current matrix which now exceeds "
841 <<
"its dimensions. Dimensions have been increased!");
842 numCol_ = maxCol_ + 1;
843 numRow_ = maxRow_ + 1;
873 if(numRow_ > numCol_){
883 MultiplyWithTransposeOf(*
this, AAT);
899 BIASASSERT(numRow_ == numCol_);
902 for (
unsigned int i = 0; i <= maxRow_; i++) {
904 rows_[i].push_back(std::pair<unsigned int, double>(maxCol_+i+1, 1));
908 int res = GaussJordan_();
909 if (res != 0)
return res;
912 inverse.
newsize(numRow_, numCol_);
914 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
916 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
917 iterator rowsEnd = rows_.end();
918 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
919 std::list<std::pair<unsigned int, double> >::iterator colsIter;
920 std::list<std::pair<unsigned int, double> >::iterator colsEnd
921 = rowsIter->second.end();
922 for (colsIter = rowsIter->second.begin();
923 colsIter != colsEnd && colsIter->first <= maxCol_; colsIter++);
925 if (colsIter == colsEnd) {
926 BIASERR(
"Sparse matrix is not consistent (contains empty rows)!");
929 for (; colsIter != colsEnd; colsIter++) {
930 inverse[rowsIter->first][colsIter->first-maxCol_-1] = colsIter->second;
935 maxCol_ = 2*maxCol_+1;
937 if ((numCol_ < maxCol_ + 1) || (numRow_ < maxRow_)){
938 BIASWARNONCE(
"Call to 'Invert' enlarged current matrix which now exceeds "
939 <<
"its dimensions. Dimensions have been increased!");
940 numCol_ = maxCol_ + 1;
941 numRow_ = maxRow_ + 1;
952 cout <<
"numRow " << numRow_ <<
" numCol " << numCol_ << endl;
953 BIASASSERT (numRow_ == numCol_);
956 for (
unsigned int i = 0; i <= maxRow_; i++) {
958 rows_[i].push_back(std::pair<unsigned int, double>(maxCol_+i+1, 1));
962 int res = GaussJordan_();
963 if (res != 0)
return res;
966 inverse.
rows_.clear();
971 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
973 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
974 iterator rowsEnd = rows_.end();
975 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
976 std::list<std::pair<unsigned int, double> >::iterator colsIter;
977 std::list<std::pair<unsigned int, double> >::iterator colsEnd
978 = rowsIter->second.end();
979 for (colsIter = rowsIter->second.begin();
980 colsIter != colsEnd && colsIter->first <= maxCol_; colsIter++);
982 if (colsIter == colsEnd) {
983 BIASERR(
"Sparse matrix is not consistent (contains empty rows)!");
986 std::list<std::pair<unsigned int, double> > row;
987 for (; colsIter != colsEnd; colsIter++) {
988 row.push_back(std::pair<unsigned int, double>(colsIter->first-maxCol_-1,
991 if (!row.empty()) inverse.
rows_[rowsIter->first] = row;
995 maxCol_ = 2*maxCol_+1;
997 if ((numCol_ < maxCol_ + 1) || (numRow_ < maxRow_ + 1)){
998 BIASWARNONCE(
"Call to 'Invert' enlarged current matrix which now exceeds "
999 <<
"its dimensions. Dimensions have been increased! maxCol " << maxCol_+1 <<
" numCol " << numCol_ <<
" maxRows " << maxRow_+1 <<
" numRows " << numRow_ << endl);
1000 numCol_ = maxCol_ + 1;
1001 numRow_ = maxRow_ + 1;
1009 #ifdef GAUSS_JORDAN_DEBUG
1010 cout <<
"-- SparseMatrix::GaussJordan_() : Debugging --" << endl << endl;
1013 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
1015 for (rowsIter1 = rows_.begin(); rowsIter1 != rows_.end(); rowsIter1++) {
1016 #ifdef GAUSS_JORDAN_DEBUG
1017 cout <<
"- Processing line " << rowsIter1->first << endl << flush;
1020 double pivot_value = 0.0;
1021 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
1022 iterator pivot_index;
1023 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
1026 for (rowsIter2 = rowsIter1; rowsIter2 != rows_.end(); rowsIter2++) {
1028 if (rowsIter2->second.empty())
continue;
1029 #ifdef GAUSS_JORDAN_DEBUG
1030 cout <<
" - Checking for row echelon form in line "
1031 << rowsIter2->first << endl << flush;
1033 std::pair<unsigned int, double> front = rowsIter2->second.front();
1036 BIASASSERT(front.first >= rowsIter1->first);
1038 if (front.first == rowsIter1->first) {
1039 if (ABS(front.second) > ABS(pivot_value)) {
1040 #ifdef GAUSS_JORDAN_DEBUG
1041 cout <<
" - Changed pivot index to " << rowsIter2->first
1044 pivot_value = front.second;
1045 pivot_index = rowsIter2;
1049 if (pivot_value == 0.0) {
1050 BIASERR(
"Matrix is singular!");
1054 if (pivot_index != rowsIter1) {
1058 #ifdef GAUSS_JORDAN_DEBUG
1059 cout <<
"- Swapping lines " << rowsIter1->first <<
" and "
1060 << pivot_index->first <<
" with sizes " << rowsIter1->second.size()
1061 <<
" and " << pivot_index->second.size() <<
"... " << flush;
1063 swap(rowsIter1->second, pivot_index->second);
1064 #ifdef GAUSS_JORDAN_DEBUG
1065 cout <<
"done!" << flush << endl;
1068 #ifdef GAUSS_JORDAN_DEBUG
1069 cout <<
"- Eliminating column downwards..." << flush << endl;
1072 rowsIter2 = rowsIter1;
1073 for (rowsIter2++; rowsIter2 != rows_.end(); rowsIter2++) {
1074 if (rowsIter2->second.empty())
continue;
1075 std::pair<unsigned int, double> front = rowsIter2->second.front();
1078 if (front.first == rowsIter1->first) {
1079 rowsIter2->second.pop_front();
1080 double factor = front.second/pivot_value;
1084 if (ABS(factor) < 1.0/(SPARSE_MATRIX_EPSILON*SPARSE_MATRIX_EPSILON)) {
1089 std::list<std::pair<unsigned int, double> >::iterator colsIter1 =
1090 rowsIter1->second.begin();
1091 std::list<std::pair<unsigned int, double> >::iterator colsIter2 =
1092 rowsIter2->second.begin();
1095 while (colsIter1 != rowsIter1->second.end()) {
1096 if (colsIter2 == rowsIter2->second.end()) {
1097 double value = -factor*colsIter1->second;
1098 if (ABS(value) > SPARSE_MATRIX_EPSILON) {
1099 rowsIter2->second.insert(colsIter2,
1100 std::pair<unsigned int, double>(colsIter1->first, value));
1103 }
else if (colsIter1->first>colsIter2->first) {
1105 }
else if (colsIter1->first<colsIter2->first) {
1106 double value = -factor*colsIter1->second;
1107 if (ABS(value) > SPARSE_MATRIX_EPSILON) {
1108 rowsIter2->second.insert(colsIter2,
1109 std::pair<unsigned int, double>(colsIter1->first, value));
1113 colsIter2->second -= factor*colsIter1->second;
1114 if (ABS(colsIter2->second) <= SPARSE_MATRIX_EPSILON) {
1115 colsIter2 = rowsIter2->second.erase(colsIter2);
1124 #ifdef GAUSS_JORDAN_DEBUG
1125 cout <<
"- Eliminating column upwards..." << flush << endl;
1128 for (rowsIter2 = rows_.begin(); rowsIter2 != rowsIter1; rowsIter2++) {
1129 std::list<std::pair<unsigned int, double> >::iterator colsIter2 =
1130 rowsIter2->second.begin();
1132 if (colsIter2 != rowsIter2->second.end()) {
1133 std::list<std::pair<unsigned int, double> >::iterator colsIter1 =
1134 rowsIter1->second.begin();
1135 if (colsIter2->first == colsIter1->first) {
1136 double factor = colsIter2->second/pivot_value;
1140 if (ABS(factor) < 1.0/(SPARSE_MATRIX_EPSILON*SPARSE_MATRIX_EPSILON)) {
1145 colsIter2 = rowsIter2->second.erase(colsIter2);
1147 while (colsIter1 != rowsIter1->second.end()) {
1148 if (colsIter2 == rowsIter2->second.end()) {
1149 double value = -factor*colsIter1->second;
1150 if (ABS(value) > SPARSE_MATRIX_EPSILON) {
1151 rowsIter2->second.insert(colsIter2,
1152 std::pair<unsigned int, double>(colsIter1->first, value));
1155 }
else if (colsIter1->first>colsIter2->first) {
1157 }
else if (colsIter1->first<colsIter2->first) {
1158 double value = -factor*colsIter1->second;
1159 if (ABS(value) > SPARSE_MATRIX_EPSILON) {
1160 rowsIter2->second.insert(colsIter2,
1161 std::pair<unsigned int, double>(colsIter1->first, value));
1165 colsIter2->second -= factor*colsIter1->second;
1166 if (ABS(colsIter2->second) <= SPARSE_MATRIX_EPSILON) {
1167 colsIter2 = rowsIter2->second.erase(colsIter2);
1177 #ifdef GAUSS_JORDAN_DEBUG
1178 cout <<
"- Normalizing row..." << flush << endl;
1181 std::list<std::pair<unsigned int, double> >::iterator colsIter1
1182 = rowsIter1->second.begin();
1183 while (colsIter1 != rowsIter1->second.end()) {
1184 colsIter1->second /= pivot_value;
1185 if (ABS(colsIter1->second) <= SPARSE_MATRIX_EPSILON) {
1186 colsIter1 = rowsIter1->second.erase(colsIter1);
1194 #ifdef GAUSS_JORDAN_DEBUG
1195 cout <<
"- Removing empty rows from resulting matrix..." << flush << endl;
1197 rowsIter1 = rows_.begin();
1198 while (rowsIter1 != rows_.end()) {
1199 if (rowsIter1->second.empty()) {
1200 rows_.erase(rowsIter1);
1211 unsigned int maxRow = 0, maxCol = 0;
1212 unsigned int lastRow = 0;
1213 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
1214 const_iterator rowsIter;
1215 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
1216 const_iterator rowsEnd = rows_.end();
1217 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
1218 if (rowsIter != rows_.begin()) {
1219 if (rowsIter->first <= lastRow) {
1220 BIASERR(
"Sparse matrix has rows in wrong order: Row "
1221 << rowsIter->first <<
" follows row " << lastRow <<
"!");
1226 lastRow = rowsIter->first;
1227 maxRow = MAX(maxRow, rowsIter->first);
1229 if (rowsIter->second.empty()) {
1230 BIASERR(
"Sparse matrix contains empty row " << rowsIter->first<<
"!");
1235 unsigned int lastCol = 0;
1236 std::list<std::pair<unsigned int, double> >::const_iterator colsIter;
1237 std::list<std::pair<unsigned int, double> >::const_iterator colsEnd
1238 = rowsIter->second.end();
1239 for (colsIter = rowsIter->second.begin();
1240 colsIter != colsEnd; colsIter++) {
1241 if (colsIter != rowsIter->second.begin()) {
1242 if (colsIter->first <= lastCol) {
1243 BIASERR(
"Sparse matrix has columns in wrong order: Column "
1244 << colsIter->first <<
" follows column " << lastCol <<
"!");
1249 lastCol = colsIter->first;
1250 maxCol = MAX(maxCol, colsIter->first);
1252 if (ABS(colsIter->second) <= SPARSE_MATRIX_EPSILON) {
1253 BIASERR(
"Sparse matrix contains to small value " << colsIter->second
1254 <<
" at [" << rowsIter->first <<
"][" << colsIter->first
1263 if (maxRow > maxRow_ || maxCol > maxCol_) {
1264 BIASERR(
"Sparse matrix contains values out of bounds (max. row ="
1265 << maxRow_ <<
", max. col = " << maxCol_ <<
", value at ["
1266 << maxRow <<
"][" << maxCol <<
"] found!)");
1275 std::vector<unsigned int> &cols,
1276 std::vector<double> &values)
const
1281 unsigned int num = 0;
1282 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
1283 const_iterator rowsIter;
1284 std::map<unsigned int, std::list<std::pair<unsigned int, double> > >::
1285 const_iterator rowsEnd = rows_.end();
1286 for (rowsIter = rows_.begin(); rowsIter != rowsEnd; rowsIter++) {
1287 std::list<std::pair<unsigned int, double> >::const_iterator colsIter;
1288 std::list<std::pair<unsigned int, double> >::const_iterator colsEnd
1289 = rowsIter->second.end();
1290 for (colsIter = rowsIter->second.begin(); colsIter != colsEnd; colsIter++) {
1291 rows.push_back(rowsIter->first);
1292 cols.push_back(colsIter->first);
1293 values.push_back(colsIter->second);
1297 BIASASSERT(rows.size() == num && cols.size() == num && values.size() == num);
double GetElement(unsigned int row, unsigned int col) const
Returns value of an element.
void GetAllElements(std::vector< unsigned int > &rows, std::vector< unsigned int > &cols, std::vector< double > &values) const
Return rows, columns and values for all elements in matrix.
bool CheckConsistency()
Checks if this sparse matrix is still consistent, i.e.
void Reinit(unsigned rows, unsigned cols)
clears old matrix and reinits with given size
int Invert(BIAS::Matrix< double > &inverse)
Returns inverse matrix as dense matrix.
void Print(std::string const &name, std::ostream &stream=std::cout) const
Prints the contents of this matrix.
unsigned int maxRow_
Current max.
Matrix< T > & newsize(Subscript M, Subscript N)
std::map< unsigned int, std::list< std::pair< unsigned int, double > > > rows_
Maps row indices to mappings from column indices to values.
void Multiply(SparseMatrix &A, SparseMatrix &result)
Multiplies with A and returns result M*A.
unsigned int Size() const
length of the vector
double GetMaxDiagonalElement() const
void SetZero()
equivalent to matrix call
unsigned int GetRows() const
Returns the maximum number of rows (max non-zero entry).
Implementation of sparse matrix operations.
Vector< T > & newsize(Subscript N)
SparseMatrix & operator=(const SparseMatrix &S)
void SetZero()
Sets all values to zero.
unsigned int numRow_
Indicates the current number of rows and clumns.
unsigned int GetRows() const
unsigned int GetCols() const
Returns the maximum number of columns (max non-zero entry).
void MultiplyIP(const double multiplier)
Multiplies this elementwise with 'multiplier'.
int GaussJordan_()
Performs Gauss-Jordan algorithm on sparse matrix.
void InsertElement_(unsigned int row, unsigned int col, double value)
Insert an element.
double GetMaxRowElement(unsigned int row, unsigned int &col) const
Returns maximum element at given row (at least zero).
unsigned int GetCols() const
void MultiplyWithTransposeOf(SparseMatrix &A, SparseMatrix &result)
Multiplies with transpose of A and returns result M*A'.
void GetAsDense(BIAS::Matrix< double > &M) const
Returns this matrix as dense BIAS matrix in M.
void MultiplyDiagonalBy(const double &value)
Runs across diagonal and multiplies diagonal element with value.
void SetSubMatrix(unsigned int row, unsigned int col, BIAS::Matrix< double > const &A, unsigned int startRowInA, unsigned int startColInA, unsigned int row_count, unsigned int col_count)
Sets a submatrix beginning at (row,col) from a dense matrix.
int Solve(BIAS::Vector< double > &b, BIAS::Vector< double > &result)
Returns solution of linear equation system A*x = b.
void AddToDiagonal(const double &value)
Runs across diagonal and adds value to diagonal elements.
SparseMatrix()
Default constructor for zero sparse matrix.
void AppendMatrixRight(SparseMatrix &H)
double GetMaxColumnElement(unsigned int col, unsigned int &row) const
Returns maximum element at given column (at least zero).
unsigned int GetRowsNum() const
Returns the total number of rows (matrix size).
void AppendMatrixBottom(SparseMatrix &H)
int PseudoInverse(BIAS::Matrix< double > &inverse)
computes pseudo inverse and returns result in inverse
void Transpose(SparseMatrix &result)
Returns transposed matrix in result (as sparse matrix).
int InvertAndSolve(BIAS::Vector< double > &b, BIAS::Matrix< double > &inverse, BIAS::Vector< double > &result)
Inverts matrix and returns inverse (as dense matrix) and solution of linear equation system A^(-1)*x ...
void SetElement(unsigned int row, unsigned int col, double val)
Sets a new element or changes the value of an existing element.
unsigned int GetColsNum() const
Returns the total number of columns (matrix size).