diff --git a/Source/lsmrBase.cxx b/Source/lsmrBase.cxx index e8bcac7..bb924f9 100755 --- a/Source/lsmrBase.cxx +++ b/Source/lsmrBase.cxx @@ -322,7 +322,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x ) std::fill( v, v+n, zero); std::fill( w, v+n, zero); std::fill( x, x+n, zero); - + this->Scale( m, (-1.0), u ); this->Aprod1( m, n, x, u ); this->Scale( m, (-1.0), u ); diff --git a/Source/lsmrBase.h b/Source/lsmrBase.h index 1f935ad..3a786eb 100755 --- a/Source/lsmrBase.h +++ b/Source/lsmrBase.h @@ -445,7 +445,7 @@ class lsmrBase */ unsigned int GetStoppingReason() const; - /** + /** * Returns an string giving the reason for termination. * Expands on GetStoppingReason * diff --git a/Source/lsmrDense.cxx b/Source/lsmrDense.cxx index 547d96e..5cf270a 100644 --- a/Source/lsmrDense.cxx +++ b/Source/lsmrDense.cxx @@ -75,5 +75,3 @@ Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const x[col] += sum; } } - - diff --git a/Source/lsmrDense.h b/Source/lsmrDense.h index 92ec6ed..c8f9795 100644 --- a/Source/lsmrDense.h +++ b/Source/lsmrDense.h @@ -24,9 +24,9 @@ /** \class lsmrDense * * Specific implementation of the solver for a type of dense Matrix. - * + * */ -class lsmrDense : public lsmrBase +class lsmrDense : public lsmrBase { public: @@ -57,4 +57,4 @@ class lsmrDense : public lsmrBase double ** A; }; -#endif +#endif diff --git a/Source/lsqrBase.cxx b/Source/lsqrBase.cxx index 4f7b9c0..b19d325 100644 --- a/Source/lsqrBase.cxx +++ b/Source/lsqrBase.cxx @@ -118,14 +118,14 @@ lsqrBase::GetFinalEstimateOfNormRbar() const } -double +double lsqrBase::GetFinalEstimateOfNormOfResiduals() const { return this->Arnorm; } -double +double lsqrBase::GetFinalEstimateOfNormOfX() const { return this->xnorm; @@ -209,7 +209,7 @@ lsqrBase::D2Norm( double a, double b ) const { return zero; } - + const double sa = a / scale; const double sb = b / scale; @@ -239,12 +239,12 @@ lsqrBase::Dnrm2( unsigned int n, const double *x ) const for ( unsigned int i = 0; i < n; i++ ) { - if ( x[i] != 0.0 ) + if ( x[i] != 0.0 ) { double dx = x[i]; const double absxi = std::abs(dx); - if ( magnitudeOfLargestElement < absxi ) + if ( magnitudeOfLargestElement < absxi ) { // rescale the sum to the range of the new element dx = magnitudeOfLargestElement / absxi; @@ -266,7 +266,7 @@ lsqrBase::Dnrm2( unsigned int n, const double *x ) const } -/** +/** * * The array b must have size m * @@ -281,16 +281,16 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x ) { (*this->nout) << "Enter LSQR " << std::endl; (*this->nout) << m << ", " << n << std::endl; - (*this->nout) << this->damp << ", " << this->wantse << std::endl; - (*this->nout) << this->atol << ", " << this->conlim << std::endl; - (*this->nout) << this->btol << ", " << this->itnlim << std::endl; + (*this->nout) << this->damp << ", " << this->wantse << std::endl; + (*this->nout) << this->atol << ", " << this->conlim << std::endl; + (*this->nout) << this->btol << ", " << this->itnlim << std::endl; } this->damped = ( this->damp > zero ); this->itn = 0; this->istop = 0; - + unsigned int nstop = 0; this->maxdx = 0; @@ -363,7 +363,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x ) return; } - + double rhobar = alpha; double phibar = beta; @@ -372,7 +372,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x ) double test1 = 0.0; double test2 = 0.0; - + if ( this->nout ) { @@ -416,7 +416,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x ) // do { - this->itn++; + this->itn++; //---------------------------------------------------------------- // Perform the next step of the bidiagonalization to obtain the @@ -488,7 +488,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x ) double t3 = one / rho; double dknorm = zero; - if ( this->wantse ) + if ( this->wantse ) { for ( unsigned int i = 0; i < n; i++ ) { @@ -649,7 +649,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x ) } } - } while ( istop == 0); + } while ( istop == 0); //=================================================================== // End of iteration loop. diff --git a/Source/lsqrBase.h b/Source/lsqrBase.h index 2da1b47..ce17331 100644 --- a/Source/lsqrBase.h +++ b/Source/lsqrBase.h @@ -26,44 +26,44 @@ * \brief implement a solver for a set of linear equations. * * LSQR finds a solution x to the following problems: - * + * * 1. Unsymmetric equations: Solve A*x = b - * + * * 2. Linear least squares: Solve A*x = b * in the least-squares sense - * + * * 3. Damped least squares: Solve ( A )*x = ( b ) * ( damp*I ) ( 0 ) * in the least-squares sense - * + * * where A is a matrix with m rows and n columns, b is an m-vector, * and damp is a scalar. (All quantities are real.) * The matrix A is treated as a linear operator. It is accessed * by means of subroutine calls with the following purpose: - * + * * call Aprod1(m,n,x,y) must compute y = y + A*x without altering x. * call Aprod2(m,n,x,y) must compute x = x + A'*y without altering y. - * + * * LSQR uses an iterative method to approximate the solution. * The number of iterations required to reach a certain accuracy * depends strongly on the scaling of the problem. Poor scaling of * the rows or columns of A should therefore be avoided where * possible. - * + * * For example, in problem 1 the solution is unaltered by * row-scaling. If a row of A is very small or large compared to * the other rows of A, the corresponding row of ( A b ) should be * scaled up or down. - * + * * In problems 1 and 2, the solution x is easily recovered * following column-scaling. Unless better information is known, * the nonzero columns of A should be scaled so that they all have * the same Euclidean norm (e.g., 1.0). - * + * * In problem 3, there is no freedom to re-scale if damp is * nonzero. However, the value of damp should be assigned only * after attention has been paid to the scaling of A. - * + * * The parameter damp is intended to help regularize * ill-conditioned systems, by preventing the true solution from * being very large. Another aid to regularization is provided by @@ -74,13 +74,13 @@ * of the solver that is available at * http://www.stanford.edu/group/SOL/software.html * distributed under a BSD license. - * + * * This class is a replacement for the lsqr code taken from netlib. * That code had to be removed because it is copyrighted by ACM and * its license was incompatible with a BSD license. - * + * */ -class lsqrBase +class lsqrBase { public: @@ -102,7 +102,7 @@ class lsqrBase * The size of the vector y is m. */ virtual void Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const = 0; - + /** * returns sqrt( a**2 + b**2 ) * with precautions to avoid overflow. @@ -114,12 +114,12 @@ class lsqrBase * with precautions to avoid overflow. */ double Dnrm2( unsigned int n, const double *x ) const; - + /** * Scale a vector by multiplying with a constant */ void Scale( unsigned int n, double factor, double *x ) const; - + /** A logical variable to say if the array se(*) of standard error estimates * should be computed. If m > n or damp > 0, the system is overdetermined * and the standard errors may be useful. (See the first LSQR reference.) @@ -127,7 +127,7 @@ class lsqrBase * storage can be saved by setting wantse = .false. and using any convenient * array for se(*), which won't be touched. If you call this method with the * flag ON, then you MUST provide a working memory array to store the standard - * error estimates, via the method SetStandardErrorEstimates() + * error estimates, via the method SetStandardErrorEstimates() */ void SetStandardErrorEstimatesFlag( bool ); @@ -137,8 +137,8 @@ class lsqrBase */ void SetToleranceA( double ); - /** An estimate of the relative error in the data - * defining the rhs b. For example, if b is + /** An estimate of the relative error in the data + * defining the rhs b. For example, if b is * accurate to about 6 digits, set btol = 1.0e-6. */ void SetToleranceB( double ); @@ -181,9 +181,9 @@ class lsqrBase * of damp in the range 0 to sqrt(eps)*norm(A) * will probably have a negligible effect. * Larger values of damp will tend to decrease - * the norm of x and reduce the number of + * the norm of x and reduce the number of * iterations required by LSQR. - * + * * The work per iteration and the storage needed * by LSQR are the same for all values of damp. * @@ -198,47 +198,47 @@ class lsqrBase */ void SetMaximumNumberOfIterations( unsigned int ); - /** + /** * If provided, a summary will be printed out to this stream during * the execution of the Solve function. */ void SetOutputStream( std::ostream & os ); - /** Provide the array where the standard error estimates will be stored. + /** Provide the array where the standard error estimates will be stored. * You MUST provide this working memory array if you turn on the computation * of standard error estimates with teh method SetStandardErrorEstimatesFlag(). */ void SetStandardErrorEstimates( double * array ); - /** + /** * Returns an integer giving the reason for termination: - * + * * 0 x = 0 is the exact solution. * No iterations were performed. - * + * * 1 The equations A*x = b are probably compatible. * Norm(A*x - b) is sufficiently small, given the * values of atol and btol. - * + * * 2 damp is zero. The system A*x = b is probably * not compatible. A least-squares solution has * been obtained that is sufficiently accurate, * given the value of atol. - * + * * 3 damp is nonzero. A damped least-squares * solution has been obtained that is sufficiently * accurate, given the value of atol. - * + * * 4 An estimate of cond(Abar) has exceeded conlim. * The system A*x = b appears to be ill-conditioned, * or there could be an error in Aprod1 or Aprod2. - * + * * 5 The iteration limit itnlim was reached. * */ unsigned int GetStoppingReason() const; - /** + /** * Returns an string giving the reason for termination. * Expands on GetStoppingReason * @@ -250,7 +250,7 @@ class lsqrBase unsigned int GetNumberOfIterationsPerformed() const; - /** + /** * An estimate of the Frobenius norm of Abar. * This is the square-root of the sum of squares * of the elements of Abar. @@ -263,7 +263,7 @@ class lsqrBase double GetFrobeniusNormEstimateOfAbar() const; - /** + /** * An estimate of cond(Abar), the condition * number of Abar. A very high value of Acond * may again indicate an error in Aprod1 or Aprod2. @@ -297,7 +297,7 @@ class lsqrBase /** * Execute the solver - * + * * solves Ax = b or min ||Ax - b|| with or without damping, * * m is the size of the input vector b @@ -338,4 +338,4 @@ class lsqrBase double * se; }; -#endif +#endif diff --git a/Source/lsqrDense.cxx b/Source/lsqrDense.cxx index 5175ae2..c62c395 100644 --- a/Source/lsqrDense.cxx +++ b/Source/lsqrDense.cxx @@ -77,7 +77,7 @@ Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const } -/* +/* returns x = (I - 2*z*z')*x. @@ -111,5 +111,3 @@ HouseholderTransformation(unsigned int n, const double * z, double * x ) const *xp++ -= scalarProduct * (*zp++); } } - - diff --git a/Source/lsqrDense.h b/Source/lsqrDense.h index c01c599..a4429f1 100644 --- a/Source/lsqrDense.h +++ b/Source/lsqrDense.h @@ -24,9 +24,9 @@ /** \class lsqrDense * * Specific implementation of the solver for a type of dense Matrix. - * + * */ -class lsqrDense : public lsqrBase +class lsqrDense : public lsqrBase { public: @@ -48,7 +48,7 @@ class lsqrDense : public lsqrBase * The size of the vector y is m. */ void Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const; - + /** Householder Transformation: reflects the vector "x" across the * hyperplane whose normal is defined by vector "z". The dimension of * the hyperspace is given by "n". */ @@ -62,4 +62,4 @@ class lsqrDense : public lsqrBase double ** A; }; -#endif +#endif diff --git a/Testing/lsqrTest1.cxx b/Testing/lsqrTest1.cxx index ffb8dcd..b8f839f 100644 --- a/Testing/lsqrTest1.cxx +++ b/Testing/lsqrTest1.cxx @@ -57,7 +57,7 @@ int lsqrTest1( int , char * [] ) const double norm =solver.Dnrm2( n1, x1 ); const double expectedNorm = sqrt(5.0); - const double ratioOfDifference = + const double ratioOfDifference = fabs( norm - expectedNorm ) / expectedNorm; if( ratioOfDifference > tolerance ) @@ -67,7 +67,7 @@ int lsqrTest1( int , char * [] ) std::cerr << "Received = " << norm << std::endl; std::cerr << "ratioOfDifference = " << ratioOfDifference << std::endl; return EXIT_FAILURE; - } + } std::cout << "Dnrm2 test 1 passed " << std::endl; } @@ -86,7 +86,7 @@ int lsqrTest1( int , char * [] ) const double norm =solver.Dnrm2( n2, x2 ); const double expectedNorm = dominantValue; - const double ratioOfDifference = + const double ratioOfDifference = fabs( norm - expectedNorm ) / expectedNorm; if( ratioOfDifference > tolerance ) @@ -96,7 +96,7 @@ int lsqrTest1( int , char * [] ) std::cerr << "Received = " << norm << std::endl; std::cerr << "ratioOfDifference = " << ratioOfDifference << std::endl; return EXIT_FAILURE; - } + } std::cout << "Dnrm2 test 2 passed " << std::endl; } @@ -108,7 +108,7 @@ int lsqrTest1( int , char * [] ) const double b = minorValue; const double d2norm = solver.D2Norm( a, b ); - const double ratioOfDifference = + const double ratioOfDifference = fabs( d2norm - dominantValue ) / dominantValue; if( ratioOfDifference > tolerance ) @@ -117,7 +117,7 @@ int lsqrTest1( int , char * [] ) std::cerr << "Expected = " << dominantValue << std::endl; std::cerr << "Received = " << d2norm << std::endl; return EXIT_FAILURE; - } + } std::cout << "D2Norm test 1 passed " << std::endl; } @@ -129,7 +129,7 @@ int lsqrTest1( int , char * [] ) const double d2norm = solver.D2Norm( a, b ); const double expectedD2norm = sqrt( a*a + b*b ); - const double ratioOfDifference = + const double ratioOfDifference = fabs( d2norm - expectedD2norm ) / expectedD2norm; if( ratioOfDifference > tolerance ) @@ -138,7 +138,7 @@ int lsqrTest1( int , char * [] ) std::cerr << "Expected = " << expectedD2norm << std::endl; std::cerr << "Received = " << d2norm << std::endl; return EXIT_FAILURE; - } + } std::cout << "D2Norm test 2 passed " << std::endl; } @@ -157,7 +157,7 @@ int lsqrTest1( int , char * [] ) std::cerr << "Expected = " << zero << std::endl; std::cerr << "Received = " << d2norm << std::endl; return EXIT_FAILURE; - } + } std::cout << "D2Norm test 3 passed " << std::endl; } @@ -182,7 +182,7 @@ int lsqrTest1( int , char * [] ) solver.SetStandardErrorEstimatesFlag( true ); - { // basic test for Scale() + { // basic test for Scale() const double factor = 8.0; const unsigned int mm = 5; double xx[mm]; @@ -196,7 +196,7 @@ int lsqrTest1( int , char * [] ) for ( unsigned int i = 0; i < mm; i++ ) { const double expectedValue = ( i*5.0*factor ); - const double ratioOfDifference = + const double ratioOfDifference = fabs( xx[i] - expectedValue ) / expectedValue; if ( ratioOfDifference > tolerance ) @@ -208,7 +208,7 @@ int lsqrTest1( int , char * [] ) } - { // basic test for Aprod1() + { // basic test for Aprod1() const unsigned int mm = 2; const unsigned int nn = 2; double xx[nn]; @@ -231,7 +231,7 @@ int lsqrTest1( int , char * [] ) std::cout << "yy = " << yy[0] << " " << yy[1] << std::endl; } - { // basic test for Aprod1() + { // basic test for Aprod1() const unsigned int mm = 2; const unsigned int nn = 3; double xx[nn]; @@ -257,7 +257,7 @@ int lsqrTest1( int , char * [] ) std::cout << "yy = " << yy[0] << " " << yy[1] << std::endl; } - { // basic test for Aprod1() + { // basic test for Aprod1() const unsigned int mm = 3; const unsigned int nn = 2; double xx[nn]; @@ -284,7 +284,7 @@ int lsqrTest1( int , char * [] ) std::cout << "yy = " << yy[0] << " " << yy[1] << " " << yy[2] << std::endl; } - { // basic test for Aprod2() + { // basic test for Aprod2() const unsigned int mm = 2; const unsigned int nn = 2; double xx[nn]; @@ -307,7 +307,7 @@ int lsqrTest1( int , char * [] ) std::cout << "xx = " << xx[0] << " " << xx[1] << std::endl; } - { // basic test for Aprod2() + { // basic test for Aprod2() const unsigned int mm = 2; const unsigned int nn = 3; double xx[nn]; @@ -333,7 +333,7 @@ int lsqrTest1( int , char * [] ) std::cout << "xx = " << xx[0] << " " << xx[1] << " " << xx[2] << std::endl; } - { // basic test for Aprod2() + { // basic test for Aprod2() const unsigned int mm = 3; const unsigned int nn = 2; double xx[nn]; @@ -360,7 +360,7 @@ int lsqrTest1( int , char * [] ) std::cout << "xx = " << xx[0] << " " << xx[1] << std::endl; } - { // basic test for Solve() + { // basic test for Solve() const unsigned int mm = 2; const unsigned int nn = 2; double bb[nn]; @@ -383,7 +383,7 @@ int lsqrTest1( int , char * [] ) std::cout << "xx = " << xx[0] << " " << xx[1] << std::endl; } - { // basic test for Solve() + { // basic test for Solve() const unsigned int mm = 2; const unsigned int nn = 2; double bb[nn]; @@ -410,7 +410,7 @@ int lsqrTest1( int , char * [] ) } - { // basic test for Solve() + { // basic test for Solve() const unsigned int mm = 2; const unsigned int nn = 2; double bb[nn];