Compute the length of the current solution by the recursion formula from...

Compute the length of the current solution by the recursion formula from Conn/Gould/Toint.  That way the trust region adapts to the preconditioner.  But now I cannot specify a custom norm for the trust-region

[[Imported from SVN: r2451]]

dune-solvers/solvers/tcgsolver.cc

@@ -79,6 +79,10 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
 
     rholast = p*b;         // orthogonalization
 
+    double solutionNormSquared = 0;
+    double solutionTimesCorrection = 0;
+    double correctionNormSquared = rholast;
+    
     // Loop until desired tolerance or maximum number of iterations is reached
     for (i=0; i<this->maxIterations_ && (error>this->tolerance_ || std::isnan(error)); i++) {
         
@@ -91,14 +95,25 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
 
         field_type energyNormSquared = EnergyNorm<MatrixType,VectorType>::normSquared(p, *matrix_);
 
+#if 0  // I used this for debugging
+        std::cout << "Computed values: " << std::endl;
+        std::cout << "   " << trustRegionScalarProduct(p,p)
+                  << "   " << trustRegionScalarProduct(*x_,p)
+                  << "   " << trustRegionScalarProduct(*x_,*x_) << std::endl;
+
+        std::cout << "Recursive values: " << std::endl;
+        std::cout << "   " << correctionNormSquared
+                  << "   " << solutionTimesCorrection
+                  << "   " << solutionNormSquared
+                  << std::endl;
+#endif
+
         // if energy is concave in the search direction go to the boundary and stop there
         if (energyNormSquared <= 0) {
 
-            // Compute r such that \eta = \eta^j + r p minimizes m_k and satisfies
-            // || \eta || = \Delta
-            field_type tau = positiveRoot(trustRegionScalarProduct(p,p),
-                                          2*trustRegionScalarProduct(*x_,p), 
-                                          trustRegionScalarProduct(*x_,*x_) - trustRegionRadius_*trustRegionRadius_);
+            field_type tau = positiveRoot(correctionNormSquared,
+                                          2*solutionTimesCorrection, 
+                                          solutionNormSquared - trustRegionRadius_*trustRegionRadius_);
 
             // add scaled correction to current iterate
             x_->axpy(tau, p);
@@ -124,15 +139,20 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
         VectorType tentativeNewIterate = *x_;
         tentativeNewIterate.axpy(alpha,p);           // update solution
 
-        std::cout << "tentative radius: " << trustRegionScalarProduct(tentativeNewIterate,tentativeNewIterate)
-                  << std::endl;
+        // Compute length of this tentative new iterate
+        double tentativeLengthSquared = solutionNormSquared 
+            + 2*alpha*solutionTimesCorrection 
+            + alpha*alpha*correctionNormSquared;
 
-        if (trustRegionScalarProduct(tentativeNewIterate,tentativeNewIterate) >= trustRegionRadius_*trustRegionRadius_) {
+//         std::cout << "tentative radius^2: " << trustRegionScalarProduct(tentativeNewIterate,tentativeNewIterate)
+//                   << std::endl;
+        std::cout << "tentative length^2: " << tentativeLengthSquared << std::endl;
 
-            // Compute r >= 0 such that \eta = \eta_j + \tau \delta_j satisfies || \eta || = \Delta
-            field_type tau = positiveRoot(trustRegionScalarProduct(p,p), 
-                                          2*trustRegionScalarProduct(*x_,p), 
-                                          trustRegionScalarProduct(*x_,*x_) - trustRegionRadius_*trustRegionRadius_);
+        if (tentativeLengthSquared >= trustRegionRadius_*trustRegionRadius_) {
+
+            field_type tau = positiveRoot(correctionNormSquared, 
+                                          2*solutionTimesCorrection, 
+                                          solutionNormSquared - trustRegionRadius_*trustRegionRadius_);
             
             
             
@@ -149,6 +169,8 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
         *x_ = tentativeNewIterate;           // update solution
         rhs_->axpy(-alpha,q);          // update defect
 
+        solutionNormSquared += 2*alpha*solutionTimesCorrection + alpha*alpha*correctionNormSquared;
+
         // determine new search direction
         q = 0;                      // clear correction
         
@@ -166,11 +188,20 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
         p += q;                     // orthogonalization with correction
         rholast = rho;              // remember rho for recurrence
 
+        // recursion for the preconditioner norm of the correction
+        solutionTimesCorrection = beta * ( solutionTimesCorrection + alpha*correctionNormSquared);
+
+        correctionNormSquared = rholast  + beta*beta * correctionNormSquared;
+
+        // //////////////////////////////////////////////////
         //   write iteration to file, if requested
+        // //////////////////////////////////////////////////
         if (this->historyBuffer_!="") 
             this->writeIterate(*x_, i);
 
+        // //////////////////////////////////////////////////
         //   Compute error        
+        // //////////////////////////////////////////////////
         double oldNorm = errorNorm_->operator()(oldSolution);
         oldSolution -= *x_;