diff --git a/Source/lsmrBase.cxx b/Source/lsmrBase.cxx
index e8bcac7..95b9f37 100755
--- a/Source/lsmrBase.cxx
+++ b/Source/lsmrBase.cxx
@@ -15,7 +15,6 @@
  *  limitations under the License.
  *
  *=========================================================================*/
-
 #include "lsmrBase.h"
 
 #include <algorithm>
@@ -41,7 +40,7 @@ lsmrBase::lsmrBase()
   this->btol = 1e-6;
   this->conlim = 1.0 / ( 10 * sqrt( this->eps ) );
   this->itnlim = 10;
-  this->nout = NULL;
+  this->nout = nullptr;
   this->istop = 0;
   this->itn = 0;
   this->normA = 0.0;
@@ -58,9 +57,7 @@ lsmrBase::lsmrBase()
 }
 
 
-lsmrBase::~lsmrBase()
-{
-}
+lsmrBase::~lsmrBase() = default;
 
 
 unsigned int
@@ -249,24 +246,24 @@ lsmrBase::Dnrm2( unsigned int n, const double *x ) const
   for ( unsigned int i = 0; i < n; i++ )
     {
       if ( x[i] != 0.0 )
-	{
-	  double dx = x[i];
-	  const double absxi = std::abs(dx);
-
-	  if ( magnitudeOfLargestElement < absxi )
-	    {
-	      // rescale the sum to the range of the new element
-	      dx = magnitudeOfLargestElement / absxi;
-	      sumOfSquaresScaled = sumOfSquaresScaled * (dx * dx) + 1.0;
-	      magnitudeOfLargestElement = absxi;
-	    }
-	  else
-	    {
-	      // rescale the new element to the range of the sum
-	      dx = absxi / magnitudeOfLargestElement;
-	      sumOfSquaresScaled += dx * dx;
-	    }
-	}
+        {
+          double dx = x[i];
+          const double absxi = std::abs(dx);
+
+          if ( magnitudeOfLargestElement < absxi )
+            {
+              // rescale the sum to the range of the new element
+              dx = magnitudeOfLargestElement / absxi;
+              sumOfSquaresScaled = sumOfSquaresScaled * (dx * dx) + 1.0;
+              magnitudeOfLargestElement = absxi;
+            }
+          else
+            {
+              // rescale the new element to the range of the sum
+              dx = absxi / magnitudeOfLargestElement;
+              sumOfSquaresScaled += dx * dx;
+            }
+        }
     }
 
   const double norm = magnitudeOfLargestElement * sqrt( sumOfSquaresScaled );
@@ -290,7 +287,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
 
   // Initialize.
 
-  unsigned int localVecs = std::min( localSize, std::min( m,n ) );
+  unsigned int const localVecs = std::min( localSize, std::min( m,n ) );
 
   if( this->nout )
     {
@@ -302,7 +299,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
     (*this->nout) << " localSize (no. of vectors for local reorthogonalization) = " << this->localSize << std::endl;
     }
 
-  int pfreq = 20;
+  int const pfreq = 20;
   int pcount = 0;
   this->damped = ( this->damp > zero );
 
@@ -322,7 +319,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 );
@@ -400,15 +397,15 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
   if ( this->nout )
     {
       if ( damped )
-	{
-	  (*this->nout) << "   Itn       x(1)           norm rbar    Abar'rbar"
-	    " Compatible    LS    norm Abar cond Abar\n";
-	}
+        {
+          (*this->nout) << "   Itn       x(1)           norm rbar    Abar'rbar"
+            " Compatible    LS    norm Abar cond Abar\n";
+        }
       else
-	{
-	  (*this->nout) << "   Itn       x(1)            norm r         A'r   "
-	    " Compatible    LS      norm A    cond A\n";
-	}
+        {
+          (*this->nout) << "   Itn       x(1)            norm r         A'r   "
+            " Compatible    LS      norm A    cond A\n";
+        }
 
       test1 = one;
       test2 = alpha / beta;
@@ -434,34 +431,34 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
 
     if ( beta > zero )
       {
-	this->Scale( m, (one/beta), u );
-	if ( localOrtho ) {
-	  if (localPointer+1 < localVecs) {
-	    localPointer = localPointer + 1;
-	  } else {
-	    localPointer = 0;
-	    localVQueueFull = true;
-	  }
-	  std::copy( v, v+n, localV+localPointer*n );
-	}
-	this->Scale( n, (- beta), v );
-	this->Aprod2( m, n, v, u );    // v = A'*u
-	if ( localOrtho ) {
-	  unsigned int localOrthoLimit = localVQueueFull ? localVecs : localPointer+1;
-
-	  for( unsigned int localOrthoCount =0; localOrthoCount<localOrthoLimit;
-	       ++localOrthoCount) {
-	    double d = std::inner_product(v,v+n,localV+n*localOrthoCount,0.0);
-	    daxpy( n, -d, localV+localOrthoCount*n, v );
-	  }
-	}
-
-	alpha  = this->Dnrm2( n, v );
-
-	if ( alpha > zero )
-	  {
+        this->Scale( m, (one/beta), u );
+        if ( localOrtho ) {
+          if (localPointer+1 < localVecs) {
+            localPointer = localPointer + 1;
+          } else {
+            localPointer = 0;
+            localVQueueFull = true;
+          }
+          std::copy( v, v+n, localV+localPointer*n );
+        }
+        this->Scale( n, (- beta), v );
+        this->Aprod2( m, n, v, u );    // v = A'*u
+        if ( localOrtho ) {
+          unsigned int localOrthoLimit = localVQueueFull ? localVecs : localPointer+1;
+
+          for( unsigned int localOrthoCount =0; localOrthoCount<localOrthoLimit;
+               ++localOrthoCount) {
+            double d = std::inner_product(v,v+n,localV+n*localOrthoCount,0.0);
+            daxpy( n, -d, localV+localOrthoCount*n, v );
+          }
+        }
+
+        alpha  = this->Dnrm2( n, v );
+
+        if ( alpha > zero )
+          {
             this->Scale( n, (one/alpha), v );
-	  }
+          }
       }
 
     // At this point, beta = beta_{k+1}, alpha = alpha_{k+1}.
@@ -470,25 +467,25 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
     //----------------------------------------------------------------
     // Construct rotation Qhat_{k,2k+1}.
 
-    double alphahat = this->D2Norm( alphabar, damp );
-    double chat     = alphabar/alphahat;
-    double shat     = damp/alphahat;
+    double const alphahat = this->D2Norm( alphabar, damp );
+    double const chat     = alphabar/alphahat;
+    double const shat     = damp/alphahat;
 
     // Use a plane rotation (Q_i) to turn B_i to R_i.
 
-    double rhoold   = rho;
+    double const rhoold   = rho;
     rho      = D2Norm(alphahat, beta);
-    double c        = alphahat/rho;
-    double s        = beta/rho;
-    double thetanew = s*alpha;
+    double const c        = alphahat/rho;
+    double const s        = beta/rho;
+    double const thetanew = s*alpha;
     alphabar = c*alpha;
 
     // Use a plane rotation (Qbar_i) to turn R_i^T into R_i^bar.
 
-    double rhobarold = rhobar;
-    double zetaold   = zeta;
-    double thetabar  = sbar*rho;
-    double rhotemp   = cbar*rho;
+    double const rhobarold = rhobar;
+    double const zetaold   = zeta;
+    double const thetabar  = sbar*rho;
+    double const rhotemp   = cbar*rho;
     rhobar    = this->D2Norm(cbar*rho, thetanew);
     cbar      = cbar*rho/rhobar;
     sbar      = thetanew/rhobar;
@@ -506,20 +503,20 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
     // Estimate ||r||.
 
     // Apply rotation Qhat_{k,2k+1}.
-    double betaacute =   chat* betadd;
-    double betacheck = - shat* betadd;
+    double const betaacute =   chat* betadd;
+    double const betacheck = - shat* betadd;
 
     // Apply rotation Q_{k,k+1}.
-    double betahat   =   c*betaacute;
+    double const betahat   =   c*betaacute;
     betadd    = - s*betaacute;
 
     // Apply rotation Qtilde_{k-1}.
     // betad = betad_{k-1} here.
 
-    double thetatildeold = thetatilde;
-    double rhotildeold   = this->D2Norm(rhodold, thetabar);
-    double ctildeold     = rhodold/rhotildeold;
-    double stildeold     = thetabar/rhotildeold;
+    double const thetatildeold = thetatilde;
+    double const rhotildeold   = this->D2Norm(rhodold, thetabar);
+    double const ctildeold     = rhodold/rhotildeold;
+    double const stildeold     = thetabar/rhotildeold;
     thetatilde    = stildeold* rhobar;
     rhodold       =   ctildeold* rhobar;
     betad         = - stildeold*betad + ctildeold*betahat;
@@ -528,7 +525,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
     // rhodold = rhod_k  here.
 
     tautildeold   = (zetaold - thetatildeold*tautildeold)/rhotildeold;
-    double taud          = (zeta - thetatilde*tautildeold)/rhodold;
+    double const taud          = (zeta - thetatilde*tautildeold)/rhodold;
     d             = d + betacheck*betacheck;
     this->normr         = sqrt(d + Sqr(betad - taud) + Sqr(betadd));
 
@@ -557,9 +554,9 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
 
     test1   = this->normr / this->normb;
     test2   = this->normAr/(this->normA*this->normr);
-    double test3   = one/this->condA;
-    double t1      = test1/(one + this->normA*this->normx/this->normb);
-    double rtol    = this->btol + this->atol*this->normA*normx/this->normb;
+    double const test3   = one/this->condA;
+    double const t1      = test1/(one + this->normA*this->normx/this->normb);
+    double const rtol    = this->btol + this->atol*this->normA*normx/this->normb;
 
     // The following tests guard against extremely small values of
     // atol, btol or ctol.  (The user may have set any or all of
@@ -593,21 +590,21 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
       if ( this->istop!=0 ) prnt = true;
 
       if ( prnt ) { // Print a line for this iteration
-	if ( pcount >= pfreq ) { // Print a heading first
-	  pcount = 0;
-	  if ( damped )
-	    {
-	      (*this->nout) << "   Itn       x(1)           norm rbar    Abar'rbar"
-		" Compatible    LS    norm Abar cond Abar\n";
-	    } else {
-	    (*this->nout) << "   Itn       x(1)            norm r         A'r   "
-	      " Compatible    LS      norm A    cond A\n";
-	  }
-	}
-	pcount = pcount + 1;
-	(*this->nout)
-	  << this->itn << ", " << x[0] << ", " <<this->normr << ", " << this->normAr << ", " << test1 << ", " << test2
-	  << ", " << this->normA << ", " << this->condA << std::endl;
+        if ( pcount >= pfreq ) { // Print a heading first
+          pcount = 0;
+          if ( damped )
+            {
+              (*this->nout) << "   Itn       x(1)           norm rbar    Abar'rbar"
+                " Compatible    LS    norm Abar cond Abar\n";
+            } else {
+            (*this->nout) << "   Itn       x(1)            norm r         A'r   "
+              " Compatible    LS      norm A    cond A\n";
+          }
+        }
+        pcount = pcount + 1;
+        (*this->nout)
+          << this->itn << ", " << x[0] << ", " <<this->normr << ", " << this->normAr << ", " << test1 << ", " << test2
+          << ", " << this->normA << ", " << this->condA << std::endl;
       }
     }
 
@@ -622,9 +619,9 @@ TerminationPrintOut()
 {
   if ( this->nout ) {
     (*this->nout) << " Exit  LSMR.       istop  = " << this->istop << "     ,itn    = " << this->itn << std::endl
-		  << " Exit  LSMR.       normA  = " << this->normA << "     ,condA  = " << this->condA << std::endl
-		  << " Exit  LSMR.       normb  = " << this->normb << "     ,normx  = " << this->normx << std::endl
-		  << " Exit  LSMR.       normr  = " << this->normr << "     ,normAr = " << this->normAr << std::endl
-		  << " Exit  LSMR.       " << this->GetStoppingReasonMessage() << std::endl;
+                  << " Exit  LSMR.       normA  = " << this->normA << "     ,condA  = " << this->condA << std::endl
+                  << " Exit  LSMR.       normb  = " << this->normb << "     ,normx  = " << this->normx << std::endl
+                  << " Exit  LSMR.       normr  = " << this->normr << "     ,normAr = " << this->normAr << std::endl
+                  << " Exit  LSMR.       " << this->GetStoppingReasonMessage() << std::endl;
   }
 }
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..485d41e 100644
--- a/Source/lsmrDense.cxx
+++ b/Source/lsmrDense.cxx
@@ -15,18 +15,15 @@
  *  limitations under the License.
  *
  *=========================================================================*/
-
 #include "lsmrDense.h"
 
 lsmrDense::lsmrDense()
 {
-  this->A = 0;
+  this->A = nullptr;
 }
 
 
-lsmrDense::~lsmrDense()
-{
-}
+lsmrDense::~lsmrDense() = default;
 
 
 void
@@ -75,5 +72,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..97025d8 100644
--- a/Source/lsmrDense.h
+++ b/Source/lsmrDense.h
@@ -24,14 +24,14 @@
 /** \class lsmrDense
  *
  * Specific implementation of the solver for a type of dense Matrix.
- *  
+ *
  */
-class lsmrDense : public lsmrBase 
+class lsmrDense : public lsmrBase
 {
 public:
 
   lsmrDense();
-  virtual ~lsmrDense();
+  ~lsmrDense() override;
 
   /**
    * computes y = y + A*x without altering x,
@@ -39,7 +39,7 @@ class lsmrDense : public lsmrBase
    * The size of the vector x is n.
    * The size of the vector y is m.
    */
-  void Aprod1(unsigned int m, unsigned int n, const double * x, double * y ) const;
+  void Aprod1(unsigned int m, unsigned int n, const double * x, double * y ) const override;
 
   /**
    * computes x = x + A'*y without altering y,
@@ -47,7 +47,7 @@ class lsmrDense : public lsmrBase
    * The size of the vector x is n.
    * The size of the vector y is m.
    */
-  void Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const;
+  void Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const override;
 
   /** Set the matrix A of the equation to be solved A*x = b. */
   void SetMatrix( double ** A );
@@ -57,4 +57,4 @@ class lsmrDense : public lsmrBase
   double ** A;
 };
 
-#endif 
+#endif
diff --git a/Source/lsqrBase.cxx b/Source/lsqrBase.cxx
index 4f7b9c0..5d4816e 100644
--- a/Source/lsqrBase.cxx
+++ b/Source/lsqrBase.cxx
@@ -28,7 +28,7 @@ lsqrBase::lsqrBase()
   this->btol = 1e-6;
   this->conlim = 1.0 / ( 10 * sqrt( this->eps ) );
   this->itnlim = 10;
-  this->nout = NULL;
+  this->nout = nullptr;
   this->istop = 0;
   this->itn = 0;
   this->Anorm = 0.0;
@@ -40,15 +40,13 @@ lsqrBase::lsqrBase()
   this->dxmax = 0.0;
   this->maxdx = 0;
   this->wantse = false;
-  this->se = NULL;
+  this->se = nullptr;
   this->damp = 0.0;
   this->damped = false;
 }
 
 
-lsqrBase::~lsqrBase()
-{
-}
+lsqrBase::~lsqrBase() = default;
 
 
 unsigned int
@@ -118,14 +116,14 @@ lsqrBase::GetFinalEstimateOfNormRbar() const
 }
 
 
-double 
+double
 lsqrBase::GetFinalEstimateOfNormOfResiduals() const
 {
   return this->Arnorm;
 }
 
 
-double 
+double
 lsqrBase::GetFinalEstimateOfNormOfX() const
 {
   return this->xnorm;
@@ -209,7 +207,7 @@ lsqrBase::D2Norm( double a, double b ) const
     {
     return zero;
     }
-  
+
   const double sa = a / scale;
   const double sb = b / scale;
 
@@ -239,12 +237,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 +264,7 @@ lsqrBase::Dnrm2( unsigned int n, const double *x ) const
 }
 
 
-/** 
+/**
  *
  *  The array b must have size m
  *
@@ -281,16 +279,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 +361,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
 
     return;
     }
-  
+
   double rhobar = alpha;
   double phibar = beta;
 
@@ -372,7 +370,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 +414,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
@@ -470,14 +468,14 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
     // Use a plane rotation to eliminate the subdiagonal element (beta)
     // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
     //----------------------------------------------------------------
-    double rho    =   this->D2Norm( rhbar1, beta );
-    double cs     =   rhbar1/rho;
-    double sn     =   beta  /rho;
-    double theta  =   sn * alpha;
+    double const rho    =   this->D2Norm( rhbar1, beta );
+    double const cs     =   rhbar1/rho;
+    double const sn     =   beta  /rho;
+    double const theta  =   sn * alpha;
     rhobar = - cs * alpha;
-    double phi    =   cs * phibar;
+    double const phi    =   cs * phibar;
     phibar =   sn * phibar;
-    double tau    =   sn * phi;
+    double const tau    =   sn * phi;
 
 
     //----------------------------------------------------------------
@@ -488,7 +486,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++ )
         {
@@ -504,7 +502,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
       {
       for ( unsigned int i = 0; i < n; i++ )
         {
-        double t = w[i];
+        double const t = w[i];
         x[i]   = t1 * t + x[i];
         w[i]   = t2 * t + v[i];
         dknorm = ( t3 * t )*( t3 * t ) + dknorm;
@@ -520,7 +518,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
     //----------------------------------------------------------------
     dknorm = sqrt( dknorm );
     dnorm  = this->D2Norm( dnorm, dknorm );
-    double dxk  = fabs( phi* dknorm );
+    double const dxk  = fabs( phi* dknorm );
     if ( this->dxmax < dxk)
       {
       this->dxmax  = dxk;
@@ -594,7 +592,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
     // See if it is time to print something.
     //----------------------------------------------------------------
     bool prnt = false;
-    if (this->nout)
+    if ( this->nout )
       {
       if (n     <=        40) prnt = true;
       if (this->itn   <=        10) prnt = true;
@@ -649,7 +647,7 @@ Solve( unsigned int m, unsigned int n, const double * b, double * x )
         }
       }
 
-    } while ( istop == 0); 
+    } while ( istop == 0);
 
   //===================================================================
   // End of iteration loop.
@@ -698,7 +696,7 @@ TerminationPrintOut()
 
   if ( this->nout )
     {
-    std::string exitt = " Exit LSQR. ";
+    std::string const exitt = " Exit LSQR. ";
 
     (*this->nout) << exitt.c_str();
     (*this->nout) << "istop = ";
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..f699016 100644
--- a/Source/lsqrDense.cxx
+++ b/Source/lsqrDense.cxx
@@ -20,13 +20,11 @@
 
 lsqrDense::lsqrDense()
 {
-  this->A = 0;
+  this->A = nullptr;
 }
 
 
-lsqrDense::~lsqrDense()
-{
-}
+lsqrDense::~lsqrDense() = default;
 
 
 void
@@ -77,7 +75,7 @@ Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const
 }
 
 
-/* 
+/*
 
   returns x = (I - 2*z*z')*x.
 
@@ -111,5 +109,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..51db3ac 100644
--- a/Source/lsqrDense.h
+++ b/Source/lsqrDense.h
@@ -24,14 +24,14 @@
 /** \class lsqrDense
  *
  * Specific implementation of the solver for a type of dense Matrix.
- *  
+ *
  */
-class lsqrDense : public lsqrBase 
+class lsqrDense : public lsqrBase
 {
 public:
 
   lsqrDense();
-  virtual ~lsqrDense();
+  ~lsqrDense() override;
 
   /**
    * computes y = y + A*x without altering x,
@@ -39,7 +39,7 @@ class lsqrDense : public lsqrBase
    * The size of the vector x is n.
    * The size of the vector y is m.
    */
-  void Aprod1(unsigned int m, unsigned int n, const double * x, double * y ) const;
+  void Aprod1(unsigned int m, unsigned int n, const double * x, double * y ) const override;
 
   /**
    * computes x = x + A'*y without altering y,
@@ -47,8 +47,8 @@ class lsqrDense : public lsqrBase
    * The size of the vector x is n.
    * The size of the vector y is m.
    */
-  void Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const;
-  
+  void Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const override;
+
   /** 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/lsmrTest2.cxx b/Testing/lsmrTest2.cxx
index 3d4eedc..3334de9 100644
--- a/Testing/lsmrTest2.cxx
+++ b/Testing/lsmrTest2.cxx
@@ -54,7 +54,7 @@ int lsmrTest2( int , char * [] )
   bb[2] = 2.0;
 
   //  -3  5
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[6];
   A[0] = &(AA[0]);
diff --git a/Testing/lsqrTest1.cxx b/Testing/lsqrTest1.cxx
index ffb8dcd..6103669 100644
--- a/Testing/lsqrTest1.cxx
+++ b/Testing/lsqrTest1.cxx
@@ -29,7 +29,7 @@ int lsqrTest1( int , char * [] )
 
   lsqrDense solver;
 
-  const double tolerance = 1e-9;
+  constexpr double tolerance = 1e-9;
 
   const unsigned int n = 2;
   double x[n];
@@ -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,13 +208,13 @@ 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];
   xx[0] = 1.0;
   xx[1] = 7.0;
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[4];
   A[0] = &(AA[0]);
@@ -231,14 +231,14 @@ 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];
   xx[0] = 1.0;
   xx[1] = 7.0;
   xx[2] = 9.0;
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[6];
   A[0] = &(AA[0]);
@@ -257,13 +257,13 @@ 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];
   xx[0] = 1.0;
   xx[1] = 7.0;
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[6];
   A[0] = &(AA[0]);
@@ -284,13 +284,13 @@ 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];
   xx[0] = 0.0;
   xx[1] = 0.0;
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[4];
   A[0] = &(AA[0]);
@@ -307,14 +307,14 @@ 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];
   xx[0] = 0.0;
   xx[1] = 0.0;
   xx[2] = 0.0;
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[6];
   A[0] = &(AA[0]);
@@ -333,13 +333,13 @@ 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];
   xx[0] = 0.0;
   xx[1] = 0.0;
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[6];
   A[0] = &(AA[0]);
@@ -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];
@@ -368,7 +368,7 @@ int lsqrTest1( int , char * [] )
   bb[1] = 7.0;
   double xx[mm];
   solver.SetStandardErrorEstimatesFlag( false );
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[4];
   A[0] = &(AA[0]);
@@ -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];
@@ -393,7 +393,7 @@ int lsqrTest1( int , char * [] )
   solver.SetStandardErrorEstimatesFlag( true );
   double se[nn];
   solver.SetStandardErrorEstimates( se );
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[4];
   A[0] = &(AA[0]);
@@ -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];
@@ -420,7 +420,7 @@ int lsqrTest1( int , char * [] )
   solver.SetStandardErrorEstimatesFlag( true );
   double se[nn];
   solver.SetStandardErrorEstimates( se );
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[4];
   A[0] = &(AA[0]);
diff --git a/Testing/lsqrTest2.cxx b/Testing/lsqrTest2.cxx
index e759de9..657e083 100644
--- a/Testing/lsqrTest2.cxx
+++ b/Testing/lsqrTest2.cxx
@@ -60,7 +60,7 @@ int lsqrTest2( int , char * [] )
   solver.SetStandardErrorEstimates( se );
 
   //  -3  5
-  typedef double * RowType;
+  using RowType = double *;
   RowType A[mm];
   double AA[6];
   A[0] = &(AA[0]);