#include <dgeequ.h>

Inheritance diagram for IterativeMatrixScale< T >:

Collaboration diagram for IterativeMatrixScale< T >:

Public Member Functions
	IterativeMatrixScale (const SolutionManager::ScaleOpt &scale)

virtual	~IterativeMatrixScale ()

Public Member Functions inherited from MatrixScale< T >
	MatrixScale (const SolutionManager::ScaleOpt &scale)

virtual	~MatrixScale ()

T &	ScaleMatrix (T &mh) const

bool	ComputeScaleFactors (const T &mh)

Public Member Functions inherited from MatrixScaleBase
	MatrixScaleBase (const SolutionManager::ScaleOpt &scale)

virtual	~MatrixScaleBase ()

VectorHandler &	ScaleRightHandSide (VectorHandler &bVH) const

VectorHandler &	ScaleSolution (VectorHandler &xVH) const

std::ostream &	Report (std::ostream &os) const

const std::vector< doublereal > &	GetRowScale () const

const std::vector< doublereal > &	GetColScale () const

bool	bGetInitialized () const

Protected Member Functions
virtual bool	ComputeScaleFactors (const T &mh, std::vector< doublereal > &rowScale, std::vector< doublereal > &colScale)

virtual std::ostream &	vReport (std::ostream &os) const

Protected Member Functions inherited from MatrixScaleBase
MatrixHandler::Norm_t	GetCondNumNorm () const

void	Prepare (const MatrixHandler &mh, integer &nrows, integer &ncols)

void	PrepareRows (const MatrixHandler &mh, integer &nrows)

void	PrepareCols (const MatrixHandler &mh, integer &ncols)

bool	bReport () const

Private Attributes
std::vector< doublereal >	DR

std::vector< doublereal >	DC

std::vector< doublereal >	normR

std::vector< doublereal >	normC

const integer	iMaxIter

const doublereal	dTol

integer	iIterTaken

doublereal	maxNormR

doublereal	maxNormC

Additional Inherited Members
Static Public Member Functions inherited from MatrixScale< T >
static MatrixScale< T > *	Allocate (const SolutionManager::ScaleOpt &scale)

Protected Attributes inherited from MatrixScaleBase
std::vector< doublereal >	rowScale

std::vector< doublereal >	colScale

doublereal	dCondBefore

doublereal	dCondAfter

const unsigned	uFlags

bool	bOK

Detailed Description

template<typename T>
class IterativeMatrixScale< T >

Definition at line 176 of file dgeequ.h.

Constructor & Destructor Documentation

template<typename T >

IterativeMatrixScale< T >::IterativeMatrixScale ( const SolutionManager::ScaleOpt & scale )

inline

Definition at line 665 of file dgeequ.h.

 :MatrixScale<T>(scale),
  iMaxIter(scale.iMaxIter),
  dTol(scale.dTol),
  iIterTaken(-1),
  maxNormR(-1.),
  maxNormC(-1.)
 {
 
 }

template<typename T >

IterativeMatrixScale< T >::~IterativeMatrixScale ( )

virtual

Definition at line 677 of file dgeequ.h.

 {
 
 }

Member Function Documentation

template<typename T >

bool IterativeMatrixScale< T >::ComputeScaleFactors	(	const T &	mh,
		std::vector< doublereal > &	rowScale,
		std::vector< doublereal > &	colScale
	)

protectedvirtual

Implements MatrixScale< T >.

Definition at line 683 of file dgeequ.h.

References ASSERT, MBDYN_EXCEPT_ARGS, and grad::sqrt().

 {
         const integer nrows = mh.iGetNumRows();
         const integer ncols = mh.iGetNumCols();
 
         if (nrows <= 0) {
                 // error
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS,
                         "invalid null or negative row number");
         }
 
         if (ncols <= 0) {
                 // error
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS,
                         "invalid null or negative column number");
         }
 
         if (rowScale.empty()) {
                 rowScale.resize(nrows);
                 normR.resize(nrows);
                 DR.resize(nrows);
                 std::fill(rowScale.begin(), rowScale.end(), 1.);
         } else if (rowScale.size() != static_cast<size_t>(nrows)) {
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS, "row number mismatch");
         } else {
                 // Use the scale values from the last Newton iteration
         }
 
         if (colScale.empty()) {
                 colScale.resize(ncols);
                 normC.resize(ncols);
                 DC.resize(ncols);
                 std::fill(colScale.begin(), colScale.end(), 1.);
         } else if (colScale.size() != static_cast<size_t>(ncols)) {
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS, "column number mismatch");
         } else {
                 // Use the scale values from the last Newton iteration
         }
 
         int i;
         bool bConverged = false, bFirstTime = true;
 
         while (true) {
                 for (i = 1; i <= iMaxIter; ++i) {
                         std::fill(DR.begin(), DR.end(), 0.);
                         std::fill(DC.begin(), DC.end(), 0.);
 
                         for (typename T::const_iterator im = mh.begin(); im != mh.end(); ++im) {
                                 doublereal d = im->dCoef;
 
                                 if (d == 0.) {
                                         continue;
                                 }
 
                                 d = std::abs(d * rowScale[im->iRow] * colScale[im->iCol]);
 
                                 if (d > DR[im->iRow]) {
                                         DR[im->iRow] = d;
                                 }
 
                                 if (d > DC[im->iCol]) {
                                         DC[im->iCol] = d;
                                 }
                         }
 
                         for (int j = 0; j < nrows; ++j) {
                                 rowScale[j] /= sqrt(DR[j]);
                         }
 
                         for (int j = 0; j < ncols; ++j) {
                                 colScale[j] /= sqrt(DC[j]);
                         }
 
                         std::fill(normR.begin(), normR.end(), 0.);
                         std::fill(normC.begin(), normC.end(), 0.);
 
                         for (typename T::const_iterator im = mh.begin(); im != mh.end(); ++im) {
                                 doublereal d = im->dCoef;
 
                                 if (d == 0.) {
                                         continue;
                                 }
 
                                 d = std::abs(d * rowScale[im->iRow] * colScale[im->iCol]);
 
                                 if (d > normR[im->iRow]) {
                                         normR[im->iRow] = d;
                                 }
 
                                 if (d > normC[im->iCol]) {
                                         normC[im->iCol] = d;
                                 }
                         }
 
                         ASSERT(!normR.empty());
                         ASSERT(!normC.empty());
 
                         maxNormR = 0.;
 
                         for (std::vector<doublereal>::const_iterator ir = normR.begin();
                                  ir != normR.end(); ++ir) {
                                 maxNormR = std::max(maxNormR, std::abs(1. - *ir));
                         }
 
                         maxNormC = 0.;
 
                         for (std::vector<doublereal>::const_iterator ic = normC.begin();
                                  ic != normC.end(); ++ic) {
                                 maxNormC = std::max(maxNormC, std::abs(1. - *ic));
                         }
 
                         if (maxNormR < dTol && maxNormC < dTol) {
                                 bConverged = true;
                                 break;
                         }
                 }
 
                 if (bConverged) {
                         break;  // Scale factors have been computed successfully!
                 } else if (bFirstTime) {
                         // No convergence: Reset the scale factors and try again!
                         std::fill(rowScale.begin(), rowScale.end(), 1.);
                         std::fill(colScale.begin(), colScale.end(), 1.);
                         bFirstTime = false;
                 } else {
                         // Still no convergence: Bail out!
                         break;
                 }
         }
 
         iIterTaken = i;
 
         return bConverged;
 }

Here is the call graph for this function:

template<typename T >

std::ostream & IterativeMatrixScale< T >::vReport ( std::ostream & os ) const

protectedvirtual

Implements MatrixScaleBase.

Definition at line 819 of file dgeequ.h.

References MatrixScaleBase::bOK.

 {
         if (!MatrixScaleBase::bOK) {
                 os << "Warning: matrix scale did not converge\n";
         }
 
         os << "row scale: " << maxNormR << std::endl
            << "col scale: " <<  maxNormC << std::endl
            << "iter scale: " << iIterTaken << std::endl;
 
         return os;
 }