MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
dgeequ.h
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/dgeequ.h,v 1.21 2017/01/12 14:43:53 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 #ifndef DGEEQU_H
33 #define DGEEQU_H
34 
35 #include <algorithm>
36 #include <cassert>
37 #include <limits>
38 #include <cmath>
39 #include <vector>
40 #include <limits>
41 #include "mh.h"
42 #include "fullmh.h"
43 #include "solman.h"
44 
45 #ifdef USE_LAPACK
46 #include "ac/lapack.h"
47 #endif // USE_LAPACK
48 
50 {
51 public:
52  explicit MatrixScaleBase(const SolutionManager::ScaleOpt& scale);
53  virtual ~MatrixScaleBase();
55  inline VectorHandler& ScaleSolution(VectorHandler& xVH) const;
56  std::ostream& Report(std::ostream& os) const;
57  const std::vector<doublereal>& GetRowScale()const{ return rowScale; }
58  const std::vector<doublereal>& GetColScale()const{ return colScale; }
59  bool bGetInitialized()const{ return !(rowScale.empty() && colScale.empty()); } // Allow row only or column only scaling
60 
61 protected:
63  virtual std::ostream& vReport(std::ostream& os) const=0;
64  inline void Prepare(const MatrixHandler& mh, integer& nrows, integer& ncols);
65  void PrepareRows(const MatrixHandler& mh, integer& nrows);
66  void PrepareCols(const MatrixHandler& mh, integer& ncols);
67  inline bool bReport() const;
68  std::vector<doublereal> rowScale, colScale;
70  const unsigned uFlags;
71  bool bOK;
72 
73 private:
74  static VectorHandler&
75  ScaleVector(VectorHandler& v, const std::vector<doublereal>& s);
76 };
77 
78 template <typename T>
80 {
81 public:
82  inline explicit MatrixScale(const SolutionManager::ScaleOpt& scale);
83  virtual ~MatrixScale();
84  inline T& ScaleMatrix(T& mh) const;
85  inline bool ComputeScaleFactors(const T& mh);
87 
88 protected:
89  virtual bool ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale)=0;
90 };
91 
92 template <typename T>
94 {
95 public:
96  inline RowSumMatrixScale(const SolutionManager::ScaleOpt& scale);
97  virtual ~RowSumMatrixScale();
98 
99 protected:
100  virtual bool ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale);
101  virtual std::ostream& vReport(std::ostream& os) const;
102 
103 private:
104  std::vector<doublereal> normRow;
105 };
106 
107 template <typename T>
109 {
110 public:
111  inline RowMaxMatrixScale(const SolutionManager::ScaleOpt& scale);
112  virtual ~RowMaxMatrixScale();
113 
114 protected:
115  virtual bool ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale);
116  virtual std::ostream& vReport(std::ostream& os) const;
117 
118 private:
119  std::vector<doublereal> normRow;
120 };
121 
122 template <typename T>
124 {
125 public:
126  inline ColSumMatrixScale(const SolutionManager::ScaleOpt& scale);
127  virtual ~ColSumMatrixScale();
128 
129 protected:
130  virtual bool ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale);
131  virtual std::ostream& vReport(std::ostream& os) const;
132 
133 private:
134  std::vector<doublereal> normCol;
135 };
136 
137 template <typename T>
139 {
140 public:
141  inline ColMaxMatrixScale(const SolutionManager::ScaleOpt& scale);
142  virtual ~ColMaxMatrixScale();
143 
144 protected:
145  virtual bool ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale);
146  virtual std::ostream& vReport(std::ostream& os) const;
147 
148 private:
149  std::vector<doublereal> normCol;
150 };
151 
152 // computes scaling factors for a matrix handler that has an iterator
153 // based on lapack's dgeequ
154 
155 template <typename T>
157 {
158 public:
159  inline LapackMatrixScale(const SolutionManager::ScaleOpt& scale);
160  virtual ~LapackMatrixScale();
161 
162 protected:
163  virtual bool ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale);
164  virtual std::ostream& vReport(std::ostream& os) const;
165 
166 private:
169 };
170 
171 // computes scaling factors for a matrix handler that has an iterator
172 // based on `A parallel Matrix Scaling Algorithm'
173 // from Patrick R. Amestoy, Iain S. Duff, Daniel Ruiz and Bora Ucar
174 
175 template <typename T>
177 {
178 public:
179  inline IterativeMatrixScale(const SolutionManager::ScaleOpt& scale);
180  virtual ~IterativeMatrixScale();
181 
182 protected:
183  virtual bool ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale);
184  virtual std::ostream& vReport(std::ostream& os) const;
185 
186 private:
187  std::vector<doublereal> DR, DC, normR, normC;
192 };
193 
194 // computes scaling factors for a matrix handler that has an iterator
195 // based on
196 // R. Sinkhorn and P. Knopp. Concerning nonnegative matrices and doubly stochastic
197 // matrices. Pacific Journal of Mathematics, 21(2): 1967.
198 
199 template <typename T>
201 {
202 public:
204  virtual ~RowMaxColMaxMatrixScale();
205 
206 protected:
207  virtual bool ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale);
208  virtual std::ostream& vReport(std::ostream& os) const;
209 
210 private:
213  std::vector<doublereal> normRow, normCol;
214 };
215 
217 {
218  if (!rowScale.empty()) {
219  ScaleVector(bVH, rowScale);
220  }
221 
222  return bVH;
223 }
224 
226 {
227  if (!colScale.empty()) {
228  ScaleVector(xVH, colScale);
229  }
230 
231  return xVH;
232 }
233 
235 {
238  return MatrixHandler::NORM_1;
239 
242 
243  default:
244  ASSERT(0);
246  }
247 }
248 
250 {
252 }
253 
254 void MatrixScaleBase::Prepare(const MatrixHandler& mh, integer& nrows, integer& ncols)
255 {
256  PrepareRows(mh, nrows);
257  PrepareCols(mh, ncols);
258 }
259 
260 template <typename T>
262  :MatrixScaleBase(scale)
263 {
264 
265 }
266 
267 template <typename T>
269 {
270 
271 }
272 
273 template <typename T>
275 {
276  if (uFlags & SolutionManager::SCALEF_COND_NUM) {
277  dCondBefore = mh.ConditionNumber(GetCondNumNorm());
278  }
279 
280  const bool bScaleRows = !rowScale.empty();
281  const bool bScaleCols = !colScale.empty();
282 
283  for (typename T::const_iterator i = mh.begin(); i != mh.end(); ++i) {
284  doublereal dCoef = i->dCoef;
285 
286  if (bScaleRows) {
287  dCoef *= rowScale[i->iRow];
288  }
289 
290  if (bScaleCols) {
291  dCoef *= colScale[i->iCol];
292  }
293 
294  mh(i->iRow + 1, i->iCol + 1) = dCoef;
295  }
296 
297  if (uFlags & SolutionManager::SCALEF_COND_NUM) {
298  dCondAfter = mh.ConditionNumber(GetCondNumNorm());
299  }
300 
301  return mh;
302 }
303 
304 template <typename T>
306 {
307  bOK = ComputeScaleFactors(mh, rowScale, colScale);
308 
309  return bOK;
310 }
311 
312 template <typename T>
314 {
315  MatrixScale<T>* pMatScale = 0;
316 
317  switch (scale.algorithm) {
319  SAFENEWWITHCONSTRUCTOR(pMatScale,
321  RowMaxMatrixScale<T>(scale));
322  break;
323 
325  SAFENEWWITHCONSTRUCTOR(pMatScale,
327  RowSumMatrixScale<T>(scale));
328  break;
329 
331  SAFENEWWITHCONSTRUCTOR(pMatScale,
333  ColMaxMatrixScale<T>(scale));
334  break;
335 
337  SAFENEWWITHCONSTRUCTOR(pMatScale,
339  ColSumMatrixScale<T>(scale));
340  break;
341 
343  SAFENEWWITHCONSTRUCTOR(pMatScale,
345  LapackMatrixScale<T>(scale));
346  break;
347 
349  SAFENEWWITHCONSTRUCTOR(pMatScale,
351  IterativeMatrixScale<T>(scale));
352  break;
353 
355  SAFENEWWITHCONSTRUCTOR(pMatScale,
358  break;
359  default:
360  ASSERT(0);
362  }
363 
364  return pMatScale;
365 }
366 
367 template <typename T>
369  :MatrixScale<T>(scale)
370 {
371 
372 }
373 
374 template <typename T>
376 {
377 
378 }
379 
380 template <typename T>
381 bool RowSumMatrixScale<T>::ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale)
382 {
383  integer nrows;
384 
385  MatrixScaleBase::PrepareRows(mh, nrows);
386 
387  if (normRow.empty()) {
388  normRow.resize(nrows, 0.);
389  } else {
390  std::fill(normRow.begin(), normRow.end(), 0.);
391  }
392 
393  for (typename T::const_iterator i = mh.begin(); i != mh.end(); ++i) {
394  normRow[i->iRow] += std::abs(i->dCoef);
395  }
396 
397  for (int i = 0; i < nrows; ++i) {
398  rowScale[i] = 1. / normRow[i];
399  }
400 
401  return true;
402 }
403 
404 template <typename T>
405 std::ostream& RowSumMatrixScale<T>::vReport(std::ostream& os) const
406 {
407  return os;
408 }
409 
410 template <typename T>
412  :MatrixScale<T>(scale)
413 {
414 
415 }
416 
417 template <typename T>
419 {
420 
421 }
422 
423 template <typename T>
424 bool RowMaxMatrixScale<T>::ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale)
425 {
426  integer nrows;
427 
428  MatrixScaleBase::PrepareRows(mh, nrows);
429 
430  if (normRow.empty()) {
431  normRow.resize(nrows, 0.);
432  } else {
433  std::fill(normRow.begin(), normRow.end(), 0.);
434  }
435 
436  for (typename T::const_iterator i = mh.begin(); i != mh.end(); ++i) {
437  const doublereal d = std::abs(i->dCoef);
438 
439  if (d > normRow[i->iRow]) {
440  normRow[i->iRow] = d;
441  }
442  }
443 
444  for (int i = 0; i < nrows; ++i) {
445  rowScale[i] = 1. / normRow[i];
446  }
447 
448  return true;
449 }
450 
451 template <typename T>
452 std::ostream& RowMaxMatrixScale<T>::vReport(std::ostream& os) const
453 {
454  return os;
455 }
456 
457 template <typename T>
459  :MatrixScale<T>(scale)
460 {
461 
462 }
463 
464 template <typename T>
466 {
467 
468 }
469 
470 template <typename T>
471 bool ColSumMatrixScale<T>::ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale)
472 {
473  integer ncols;
474 
475  MatrixScaleBase::PrepareCols(mh, ncols);
476 
477  if (normCol.empty()) {
478  normCol.resize(ncols, 0.);
479  } else {
480  std::fill(normCol.begin(), normCol.end(), 0.);
481  }
482 
483  for (typename T::const_iterator i = mh.begin(); i != mh.end(); ++i) {
484  normCol[i->iCol] += std::abs(i->dCoef);
485  }
486 
487  for (int i = 0; i < ncols; ++i) {
488  colScale[i] = 1. / normCol[i];
489  }
490 
491  return true;
492 }
493 
494 template <typename T>
495 std::ostream& ColSumMatrixScale<T>::vReport(std::ostream& os) const
496 {
497  return os;
498 }
499 
500 template <typename T>
502  :MatrixScale<T>(scale)
503 {
504 
505 }
506 
507 template <typename T>
509 {
510 
511 }
512 
513 template <typename T>
514 bool ColMaxMatrixScale<T>::ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale)
515 {
516  integer ncols;
517 
518  MatrixScaleBase::PrepareCols(mh, ncols);
519 
520  if (normCol.empty()) {
521  normCol.resize(ncols, 0.);
522  } else {
523  std::fill(normCol.begin(), normCol.end(), 0.);
524  }
525 
526  for (typename T::const_iterator i = mh.begin(); i != mh.end(); ++i) {
527  const doublereal d = std::abs(i->dCoef);
528 
529  if (d > normCol[i->iCol]) {
530  normCol[i->iCol] = d;
531  }
532  }
533 
534  for (int i = 0; i < ncols; ++i) {
535  colScale[i] = 1. / normCol[i];
536  }
537 
538  return true;
539 }
540 
541 template <typename T>
542 std::ostream& ColMaxMatrixScale<T>::vReport(std::ostream& os) const
543 {
544  return os;
545 }
546 
547 template <typename T>
549  :MatrixScale<T>(scale)
550 {
551 #if defined(HAVE_DLAMCH) || defined(HAVE_DLAMCH_)
552  // Use dlamch according to Netlib's dgeequ.f
553  SMLNUM = __FC_DECL__(dlamch)("S");
554 #else
555  // According to Netlib's dlamch for x86-64 machines
556  SMLNUM = std::numeric_limits<doublereal>::min();
557 #endif
558  BIGNUM = 1./SMLNUM;
559 
560  rowcnd = colcnd = amax = -1;
561 }
562 
563 template <typename T>
565 {
566 
567 }
568 
569 template <typename T>
570 bool LapackMatrixScale<T>::ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale)
571 {
572  integer nrows, ncols;
573  MatrixScale<T>::Prepare(mh, nrows, ncols);
574 
575  doublereal rcmin;
576  doublereal rcmax;
577 
578  for (typename T::const_iterator i = mh.begin(); i != mh.end(); ++i) {
579  doublereal d = std::abs(i->dCoef);
580  if (d > rowScale[i->iRow]) {
581  rowScale[i->iRow] = d;
582  }
583  }
584 
585  rcmin = BIGNUM;
586  rcmax = 0.;
587  for (std::vector<doublereal>::iterator i = rowScale.begin(); i != rowScale.end(); ++i) {
588  if (*i > rcmax) {
589  rcmax = *i;
590  }
591  if (*i < rcmin) {
592  rcmin = *i;
593  }
594  }
595 
596  amax = rcmax;
597 
598  if (rcmin == 0.) {
600  "null min row value in dgeequ");
601  }
602 
603  for (std::vector<doublereal>::iterator i = rowScale.begin(); i != rowScale.end(); ++i) {
604  *i = 1./(std::min(std::max(*i, SMLNUM), BIGNUM));
605  }
606 
607  rowcnd = std::max(rcmin, SMLNUM)/std::min(rcmax, BIGNUM);
608 
609  for (typename T::const_iterator i = mh.begin(); i != mh.end(); ++i) {
610  doublereal d = std::abs(i->dCoef)*rowScale[i->iRow];
611  if (d > colScale[i->iCol]) {
612  colScale[i->iCol] = d;
613  }
614  }
615 
616  rcmin = BIGNUM;
617  rcmax = 0.;
618  for (std::vector<doublereal>::iterator i = colScale.begin(); i != colScale.end(); ++i) {
619  if (*i > rcmax) {
620  rcmax = *i;
621  }
622  if (*i < rcmin) {
623  rcmin = *i;
624  }
625  }
626 
627  if (rcmin == 0.) {
629  "null min column value in dgeequ");
630  }
631 
632  for (std::vector<doublereal>::iterator i = colScale.begin(); i != colScale.end(); ++i) {
633  *i = 1./(std::min(std::max(*i, SMLNUM), BIGNUM));
634  }
635 
636  colcnd = std::max(rcmin, SMLNUM)/std::min(rcmax, BIGNUM);
637 
638  return true;
639 }
640 
641 template <typename T>
642 std::ostream& LapackMatrixScale<T>::vReport(std::ostream& os) const
643 {
644  if (amax < std::numeric_limits<doublereal>::epsilon()
645  || amax > 1./std::numeric_limits<doublereal>::epsilon())
646  {
647  os << "Warning: The matrix should be scaled\n";
648  }
649 
650  if (colcnd >= 0.1) {
651  os << "Warning: it is not worth scaling the columns\n";
652  }
653 
654  if (rowcnd >= 0.1
655  && amax >= std::numeric_limits<doublereal>::epsilon()
656  && amax <= 1./std::numeric_limits<doublereal>::epsilon())
657  {
658  os << "Warning: it is not worth scaling the rows\n";
659  }
660 
661  return os;
662 }
663 
664 template <typename T>
666 :MatrixScale<T>(scale),
667  iMaxIter(scale.iMaxIter),
668  dTol(scale.dTol),
669  iIterTaken(-1),
670  maxNormR(-1.),
671  maxNormC(-1.)
672 {
673 
674 }
675 
676 template <typename T>
678 {
679 
680 }
681 
682 template <typename T>
683 bool IterativeMatrixScale<T>::ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale)
684 {
685  const integer nrows = mh.iGetNumRows();
686  const integer ncols = mh.iGetNumCols();
687 
688  if (nrows <= 0) {
689  // error
691  "invalid null or negative row number");
692  }
693 
694  if (ncols <= 0) {
695  // error
697  "invalid null or negative column number");
698  }
699 
700  if (rowScale.empty()) {
701  rowScale.resize(nrows);
702  normR.resize(nrows);
703  DR.resize(nrows);
704  std::fill(rowScale.begin(), rowScale.end(), 1.);
705  } else if (rowScale.size() != static_cast<size_t>(nrows)) {
706  throw ErrGeneric(MBDYN_EXCEPT_ARGS, "row number mismatch");
707  } else {
708  // Use the scale values from the last Newton iteration
709  }
710 
711  if (colScale.empty()) {
712  colScale.resize(ncols);
713  normC.resize(ncols);
714  DC.resize(ncols);
715  std::fill(colScale.begin(), colScale.end(), 1.);
716  } else if (colScale.size() != static_cast<size_t>(ncols)) {
717  throw ErrGeneric(MBDYN_EXCEPT_ARGS, "column number mismatch");
718  } else {
719  // Use the scale values from the last Newton iteration
720  }
721 
722  int i;
723  bool bConverged = false, bFirstTime = true;
724 
725  while (true) {
726  for (i = 1; i <= iMaxIter; ++i) {
727  std::fill(DR.begin(), DR.end(), 0.);
728  std::fill(DC.begin(), DC.end(), 0.);
729 
730  for (typename T::const_iterator im = mh.begin(); im != mh.end(); ++im) {
731  doublereal d = im->dCoef;
732 
733  if (d == 0.) {
734  continue;
735  }
736 
737  d = std::abs(d * rowScale[im->iRow] * colScale[im->iCol]);
738 
739  if (d > DR[im->iRow]) {
740  DR[im->iRow] = d;
741  }
742 
743  if (d > DC[im->iCol]) {
744  DC[im->iCol] = d;
745  }
746  }
747 
748  for (int j = 0; j < nrows; ++j) {
749  rowScale[j] /= sqrt(DR[j]);
750  }
751 
752  for (int j = 0; j < ncols; ++j) {
753  colScale[j] /= sqrt(DC[j]);
754  }
755 
756  std::fill(normR.begin(), normR.end(), 0.);
757  std::fill(normC.begin(), normC.end(), 0.);
758 
759  for (typename T::const_iterator im = mh.begin(); im != mh.end(); ++im) {
760  doublereal d = im->dCoef;
761 
762  if (d == 0.) {
763  continue;
764  }
765 
766  d = std::abs(d * rowScale[im->iRow] * colScale[im->iCol]);
767 
768  if (d > normR[im->iRow]) {
769  normR[im->iRow] = d;
770  }
771 
772  if (d > normC[im->iCol]) {
773  normC[im->iCol] = d;
774  }
775  }
776 
777  ASSERT(!normR.empty());
778  ASSERT(!normC.empty());
779 
780  maxNormR = 0.;
781 
782  for (std::vector<doublereal>::const_iterator ir = normR.begin();
783  ir != normR.end(); ++ir) {
784  maxNormR = std::max(maxNormR, std::abs(1. - *ir));
785  }
786 
787  maxNormC = 0.;
788 
789  for (std::vector<doublereal>::const_iterator ic = normC.begin();
790  ic != normC.end(); ++ic) {
791  maxNormC = std::max(maxNormC, std::abs(1. - *ic));
792  }
793 
794  if (maxNormR < dTol && maxNormC < dTol) {
795  bConverged = true;
796  break;
797  }
798  }
799 
800  if (bConverged) {
801  break; // Scale factors have been computed successfully!
802  } else if (bFirstTime) {
803  // No convergence: Reset the scale factors and try again!
804  std::fill(rowScale.begin(), rowScale.end(), 1.);
805  std::fill(colScale.begin(), colScale.end(), 1.);
806  bFirstTime = false;
807  } else {
808  // Still no convergence: Bail out!
809  break;
810  }
811  }
812 
813  iIterTaken = i;
814 
815  return bConverged;
816 }
817 
818 template <typename T>
819 std::ostream& IterativeMatrixScale<T>::vReport(std::ostream& os) const
820 {
821  if (!MatrixScaleBase::bOK) {
822  os << "Warning: matrix scale did not converge\n";
823  }
824 
825  os << "row scale: " << maxNormR << std::endl
826  << "col scale: " << maxNormC << std::endl
827  << "iter scale: " << iIterTaken << std::endl;
828 
829  return os;
830 }
831 
832 template <typename T>
834 :MatrixScale<T>(scale),
835  dTol(scale.dTol),
836  maxNormRow(-1.),
837  maxNormCol(-1.)
838 {
839 
840 }
841 
842 template <typename T>
844 {
845 
846 }
847 
848 template <typename T>
849 bool RowMaxColMaxMatrixScale<T>::ComputeScaleFactors(const T& mh, std::vector<doublereal>& rowScale, std::vector<doublereal>& colScale)
850 {
851  const integer nrows = mh.iGetNumRows();
852  const integer ncols = mh.iGetNumCols();
853 
854  if (nrows <= 0) {
856  "invalid null or negative row number");
857  }
858 
859  if (ncols <= 0) {
861  "invalid null or negative column number");
862  }
863 
864  if (rowScale.empty()) {
865  rowScale.resize(nrows, 1.);
866  normRow.resize(nrows);
867  } else if (rowScale.size() != static_cast<size_t>(nrows)) {
868  throw ErrGeneric(MBDYN_EXCEPT_ARGS, "row number mismatch");
869  }
870 
871  if (colScale.empty()) {
872  colScale.resize(ncols, 1.);
873  normCol.resize(ncols);
874  } else if (colScale.size() != static_cast<size_t>(ncols)) {
875  throw ErrGeneric(MBDYN_EXCEPT_ARGS, "column number mismatch");
876  }
877 
878  std::fill(normRow.begin(), normRow.end(), 0.);
879 
880  for (typename T::const_iterator im = mh.begin(); im != mh.end(); ++im) {
881  doublereal d = im->dCoef;
882 
883  if (d == 0.) {
884  continue;
885  }
886 
887  d = std::abs(d);
888 
889  if (d > normRow[im->iRow]) {
890  normRow[im->iRow] = d;
891  }
892  }
893 
894  for (integer i = 0; i < nrows; ++i) {
895  rowScale[i] = 1. / normRow[i];
896  }
897 
898  std::fill(normCol.begin(), normCol.end(), 0.);
899 
900  for (typename T::const_iterator im = mh.begin(); im != mh.end(); ++im) {
901  doublereal d = im->dCoef;
902 
903  if (d == 0.) {
904  continue;
905  }
906 
907  d = std::abs(rowScale[im->iRow] * d);
908 
909  if (d > normCol[im->iCol]) {
910  normCol[im->iCol] = d;
911  }
912  }
913 
914  for (integer i = 0; i < ncols; ++i) {
915  colScale[i] = 1. / normCol[i];
916  }
917 
918  if (!this->bReport()) {
919  // Test for convergence will not be checked
920  // It would be a waste of time to compute the row and column norms
921  return true;
922  }
923 
924  std::fill(normRow.begin(), normRow.end(), 0.);
925  std::fill(normCol.begin(), normCol.end(), 0.);
926 
927  for (typename T::const_iterator im = mh.begin(); im != mh.end(); ++im) {
928  doublereal d = im->dCoef;
929 
930  if (d == 0.) {
931  continue;
932  }
933 
934  d = std::abs(colScale[im->iCol] * rowScale[im->iRow] * d);
935 
936  if (d > normRow[im->iRow]) {
937  normRow[im->iRow] = d;
938  }
939 
940  if (d > normCol[im->iCol]) {
941  normCol[im->iCol] = d;
942  }
943  }
944 
945  maxNormRow = 0.;
946 
947  for (std::vector<doublereal>::const_iterator ir = normRow.begin();
948  ir != normRow.end(); ++ir) {
949  maxNormRow = std::max(maxNormRow, std::abs(1. - *ir));
950  }
951 
952  maxNormCol = 0.;
953 
954  for (std::vector<doublereal>::const_iterator ic = normCol.begin();
955  ic != normCol.end(); ++ic) {
956  maxNormCol = std::max(maxNormCol, std::abs(1. - *ic));
957  }
958 
959  return (maxNormRow < dTol && maxNormCol < dTol);
960 }
961 
962 template <typename T>
963 std::ostream& RowMaxColMaxMatrixScale<T>::vReport(std::ostream& os) const
964 {
965  if (!MatrixScaleBase::bOK) {
966  os << "Warning: matrix scale did not converge\n";
967  }
968 
969  os << "row scale: " << maxNormRow << std::endl
970  << "col scale: " << maxNormCol << std::endl;
971 
972  return os;
973 }
974 
975 #endif // DGEEQU_H
virtual ~ColSumMatrixScale()
Definition: dgeequ.h:465
const integer iMaxIter
Definition: dgeequ.h:188
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
bool bReport() const
Definition: dgeequ.h:249
doublereal amax
Definition: dgeequ.h:168
doublereal dCondBefore
Definition: dgeequ.h:69
const unsigned uFlags
Definition: dgeequ.h:70
MatrixHandler::Norm_t GetCondNumNorm() const
Definition: dgeequ.h:234
std::vector< doublereal > colScale
Definition: dgeequ.h:68
integer iIterTaken
Definition: dgeequ.h:190
bool bGetInitialized() const
Definition: dgeequ.h:59
virtual bool ComputeScaleFactors(const T &mh, std::vector< doublereal > &rowScale, std::vector< doublereal > &colScale)
Definition: dgeequ.h:514
std::vector< doublereal > DC
Definition: dgeequ.h:187
bool ComputeScaleFactors(const T &mh)
Definition: dgeequ.h:305
MatrixScale(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:261
std::vector< doublereal > normCol
Definition: dgeequ.h:213
const std::vector< doublereal > & GetRowScale() const
Definition: dgeequ.h:57
MatrixScaleBase(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.cc:38
std::ostream & Report(std::ostream &os) const
Definition: dgeequ.cc:52
static VectorHandler & ScaleVector(VectorHandler &v, const std::vector< doublereal > &s)
Definition: dgeequ.cc:68
virtual ~RowMaxMatrixScale()
Definition: dgeequ.h:418
RowMaxColMaxMatrixScale(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:833
std::vector< doublereal > normRow
Definition: dgeequ.h:213
void PrepareCols(const MatrixHandler &mh, integer &ncols)
Definition: dgeequ.cc:104
virtual ~MatrixScaleBase()
Definition: dgeequ.cc:47
virtual bool ComputeScaleFactors(const T &mh, std::vector< doublereal > &rowScale, std::vector< doublereal > &colScale)
Definition: dgeequ.h:424
const doublereal dTol
Definition: dgeequ.h:189
virtual std::ostream & vReport(std::ostream &os) const
Definition: dgeequ.h:495
RowSumMatrixScale(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:368
virtual bool ComputeScaleFactors(const T &mh, std::vector< doublereal > &rowScale, std::vector< doublereal > &colScale)
Definition: dgeequ.h:381
const std::vector< doublereal > & GetColScale() const
Definition: dgeequ.h:58
std::vector< doublereal > normCol
Definition: dgeequ.h:134
std::vector< doublereal > normRow
Definition: dgeequ.h:119
doublereal rowcnd
Definition: dgeequ.h:168
virtual bool ComputeScaleFactors(const T &mh, std::vector< doublereal > &rowScale, std::vector< doublereal > &colScale)
Definition: dgeequ.h:570
void Prepare(const MatrixHandler &mh, integer &nrows, integer &ncols)
Definition: dgeequ.h:254
doublereal maxNormCol
Definition: dgeequ.h:212
virtual ~LapackMatrixScale()
Definition: dgeequ.h:564
doublereal maxNormRow
Definition: dgeequ.h:212
ColSumMatrixScale(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:458
doublereal dCondAfter
Definition: dgeequ.h:69
virtual std::ostream & vReport(std::ostream &os) const
Definition: dgeequ.h:542
std::vector< doublereal > normR
Definition: dgeequ.h:187
VectorHandler & ScaleRightHandSide(VectorHandler &bVH) const
Definition: dgeequ.h:216
virtual std::ostream & vReport(std::ostream &os) const
Definition: dgeequ.h:405
virtual bool ComputeScaleFactors(const T &mh, std::vector< doublereal > &rowScale, std::vector< doublereal > &colScale)
Definition: dgeequ.h:849
std::vector< doublereal > DR
Definition: dgeequ.h:187
virtual std::ostream & vReport(std::ostream &os) const
Definition: dgeequ.h:819
virtual ~RowMaxColMaxMatrixScale()
Definition: dgeequ.h:843
#define ASSERT(expression)
Definition: colamd.c:977
virtual ~IterativeMatrixScale()
Definition: dgeequ.h:677
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
RowMaxMatrixScale(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:411
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
std::vector< doublereal > rowScale
Definition: dgeequ.h:68
IterativeMatrixScale(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:665
const doublereal dTol
Definition: dgeequ.h:211
LapackMatrixScale(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:548
virtual std::ostream & vReport(std::ostream &os) const
Definition: dgeequ.h:642
virtual ~MatrixScale()
Definition: dgeequ.h:268
virtual bool ComputeScaleFactors(const T &mh, std::vector< doublereal > &rowScale, std::vector< doublereal > &colScale)
Definition: dgeequ.h:683
doublereal maxNormR
Definition: dgeequ.h:191
virtual ~RowSumMatrixScale()
Definition: dgeequ.h:375
static MatrixScale< T > * Allocate(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:313
std::vector< doublereal > normRow
Definition: dgeequ.h:104
virtual std::ostream & vReport(std::ostream &os) const =0
void PrepareRows(const MatrixHandler &mh, integer &nrows)
Definition: dgeequ.cc:82
std::vector< doublereal > normC
Definition: dgeequ.h:187
virtual ~ColMaxMatrixScale()
Definition: dgeequ.h:508
ScaleAlgorithm algorithm
Definition: solman.h:145
doublereal SMLNUM
Definition: dgeequ.h:167
ColMaxMatrixScale(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:501
doublereal maxNormC
Definition: dgeequ.h:191
doublereal colcnd
Definition: dgeequ.h:168
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
doublereal BIGNUM
Definition: dgeequ.h:167
virtual std::ostream & vReport(std::ostream &os) const
Definition: dgeequ.h:963
virtual bool ComputeScaleFactors(const T &mh, std::vector< doublereal > &rowScale, std::vector< doublereal > &colScale)
Definition: dgeequ.h:471
VectorHandler & ScaleSolution(VectorHandler &xVH) const
Definition: dgeequ.h:225
std::vector< doublereal > normCol
Definition: dgeequ.h:149
virtual std::ostream & vReport(std::ostream &os) const
Definition: dgeequ.h:452
T & ScaleMatrix(T &mh) const
Definition: dgeequ.h:274