Belos Version of the Day
BelosNewtonOpAkx.hpp
Go to the documentation of this file.
00001 #ifndef __Belos_NewtonOpAkx_hpp
00002 #define __Belos_NewtonOpAkx_hpp
00003 
00004 //@HEADER
00005 // ************************************************************************
00006 //
00007 //                 Belos: Block Linear Solvers Package
00008 //                  Copyright 2004 Sandia Corporation
00009 //
00010 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00011 // the U.S. Government retains certain rights in this software.
00012 //
00013 // Redistribution and use in source and binary forms, with or without
00014 // modification, are permitted provided that the following conditions are
00015 // met:
00016 //
00017 // 1. Redistributions of source code must retain the above copyright
00018 // notice, this list of conditions and the following disclaimer.
00019 //
00020 // 2. Redistributions in binary form must reproduce the above copyright
00021 // notice, this list of conditions and the following disclaimer in the
00022 // documentation and/or other materials provided with the distribution.
00023 //
00024 // 3. Neither the name of the Corporation nor the names of the
00025 // contributors may be used to endorse or promote products derived from
00026 // this software without specific prior written permission.
00027 //
00028 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00029 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00030 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00031 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00032 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00033 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00034 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00035 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00036 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00037 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00038 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00039 //
00040 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00041 //
00042 // ************************************************************************
00043 //@HEADER
00044 
00045 #include <BelosAkx.hpp>
00046 #include <BelosNewtonBasis.hpp>
00047 #include <Teuchos_BLAS.hpp>
00048 
00049 // #include <numeric>
00050 
00051 namespace Belos {
00052 
00060   template<class Scalar, class MV, class OP>
00061   class NewtonOpAkx : public OpAkx<Scalar, MV, OP> {
00062   public:
00063     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00064 
00065     NewtonOpAkx (const Teuchos::RCP<const OP>& A, 
00066      const Teuchos::RCP<const OP>& M_left,
00067      const Teuchos::RCP<const OP>& M_right,
00068      const bool modified) :
00069       OpAkx (A, M_left, M_right), modified_ (modified) {}
00070 
00071     void 
00072     computeBasis (const MV& q_last, MV& V_cur, const int s) 
00073     {
00074       using Teuchos::Range1D;
00075       using Teuchos::RCP;
00076       using Teuchos::rcp_const_cast;
00077       typedef MultiVecTraits<Scalar, MV> MVT;
00078       typedef Teuchos::ScalarTraits<Scalar> STS;
00079 
00080       TEUCHOS_TEST_FOR_EXCEPTION(s < 0, std::invalid_argument, 
00081        "Number of basis vectors to compute (s = " 
00082        << s << ") is invalid; it must be positive.");
00083       if (s == 0)
00084   return; // Nothing to do
00085       TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(V_cur) < s, std::invalid_argument,
00086        "You're asking to compute s = " << s << " basis "
00087        "vectors, but only have room for " 
00088        << MVT::GetNumberVecs(V_cur) << ".");
00089       RCP<const MV> v_prv = Teuchos::rcpFromRef (q_last);
00090       int k = 0;
00091       while (k < s)
00092   {
00093     RCP<MV> v_cur = MVT::CloneViewNonConst (V_cur, Range1D(k,k));
00094     if (imagParts_[k] == STS::zero())
00095       { // The current shift is real-valued.
00096         applyOp (v_prv, v_cur);
00097         MVT::MvAddMv (-realParts_[k], v_prv, STS::one(), v_cur, v_cur);
00098         // Shift the v_prv view over to the next column of
00099         // V_cur.  v_cur gets created on demand, so we don't
00100         // have to reassign it.
00101         v_prv = rcp_const_cast<const MV> (v_cur);
00102         ++k;
00103       }
00104     else // The current shift has nonzero imaginary part.
00105       {  
00106         // "Modified" Newton basis trick to ensure real arithmetic,
00107         // even if the shifts are complex (as long as any complex
00108         // shifts occur in consecutive complex conjugate pairs with
00109         // the positive-real-part member of the pair first).
00110         //
00111         // v_cur := (A - \Re \theta I) v_prv
00112         applyOp (v_prv, v_cur);
00113         MVT::MvAddMv (-realParts_[k], v_prv, STS::one(), v_cur, v_cur);
00114         // v_nxt := (A - \Re \theta I) v_cur + (\Im \theta)^2 v_prv
00115         RCP<MV> v_nxt = MVT::CloneViewNonConst (V_cur, Range1D(k+1, k+1));
00116         applyOp (v_nxt, v_cur);
00117         MVT::MvAddMv (-realParts_[k], v_cur, STS::one(), v_nxt, v_nxt);
00118         // FIXME (mfh 14 Feb 2011)
00119         //
00120         // You might think that we should be concerned that
00121         // squaring the imaginary part of the current shift
00122         // might overflow.  This is indeed possible, but the
00123         // nonmodified Newton basis doesn't really promise any
00124         // better (two steps more or less multiply by the square
00125         // of the magnitude).  This probably needs some more
00126         // analysis, but I don't think it's really broken.
00127         {
00128     const magnitude_type imSqu = imagParts_[k] * imagParts_[k];
00129     MVT::MvAddMv (imSqu, v_prv, STS::one(), v_nxt, v_nxt);
00130         }
00131         // Shift the v_prv view over to the next column of
00132         // V_cur.  v_cur and v_nxt get created on demand, so we
00133         // don't have to reassign them.
00134         v_prv = rcp_const_cast<const MV> (v_cur);
00135         k += 2;
00136       }
00137   }
00138     }
00139 
00140     void 
00141     computeFlexibleBasis (const MV& q_last, 
00142         MV& Z_cur, 
00143         MV& V_cur, 
00144         const int s)
00145     {
00146       using Teuchos::Range1D;
00147       using Teuchos::RCP;
00148       using Teuchos::rcp_const_cast;
00149       typedef MultiVecTraits<Scalar, MV> MVT;
00150 
00151       TEUCHOS_TEST_FOR_EXCEPTION(s < 0, std::invalid_argument, 
00152        "Number of basis vectors to compute (s = " 
00153        << s << ") is invalid; it must be positive.");
00154       if (s == 0)
00155   return; // Nothing to do
00156       TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(V_cur) < s, std::invalid_argument,
00157        "You're asking to compute s = " << s << " basis "
00158        "vectors, but only have room for " 
00159        << MVT::GetNumberVecs(V_cur) << " in V_cur.");
00160       TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(Z_cur) < s, std::invalid_argument,
00161        "You're asking to compute s = " << s << " basis "
00162        "vectors, but only have room for " 
00163        << MVT::GetNumberVecs(Z_cur) << " in Z_cur.");
00164       RCP<const MV> v_prv = Teuchos::rcpFromRef (q_last);
00165       int k = 0;
00166       while (k < s)
00167   {
00168     RCP<MV> z_cur = MVT::CloneViewNonConst (Z_cur, Range1D(k,k));
00169     RCP<MV> v_cur = MVT::CloneViewNonConst (V_cur, Range1D(k,k));
00170     applyFlexibleOp (v_prv, z_cur, v_cur);
00171     if (imagParts_[k] == STS::zero())
00172       {
00173         // We've already applied the full operator to v_prv,
00174         // producing v_cur.  Now we just need to apply the shift
00175         // to v_cur.
00176         MVT::MvAddMv (-realParts_[k], v_prv, STS::one(), v_cur, v_cur);
00177         // Shift the v_prv view over to the next column of
00178         // V_cur.  z_cur and v_cur get created on demand, so we
00179         // don't have to reassign them.
00180         v_prv = rcp_const_cast<const MV> (v_cur);
00181         ++k;
00182       }
00183     else
00184       {
00185         // Apply the first shift (realParts_[k]) to v_cur.
00186         MVT::MvAddMv (-realParts_[k], v_prv, STS::one(), v_cur, v_cur);
00187         RCP<MV> z_nxt = MVT::CloneViewNonConst (Z_cur, Range1D(k+1, k+1));
00188         RCP<MV> v_nxt = MVT::CloneViewNonConst (V_cur, Range1D(k+1, k+1));
00189         // After the next line, z_nxt will be just the result of
00190         // applying the right preconditioner to v_cur (to which
00191         // the first shift has already been applied).  Thus,
00192         // it's only necessary to apply the second shift to
00193         // v_nxt.
00194         applyFlexibleOp (v_cur, z_nxt, v_nxt);
00195         MVT::MvAddMv (-realParts_[k], v_cur, STS::one(), v_nxt, v_nxt);
00196         {
00197     const magnitude_type imSqu = imagParts_[k] * imagParts_[k];
00198     MVT::MvAddMv (imSqu, v_prv, STS::one(), v_nxt, v_nxt);
00199         }
00200         // Shift the v_prv view over to the next column of
00201         // V_cur.  z_cur, v_cur, z_nxt, and v_nxt get created on
00202         // demand, so we don't have to reassign them.
00203         v_prv = rcp_const_cast<const MV> (v_cur);
00204         k += 2;
00205       }
00206   }
00207     }
00208 
00212     std::pair<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,Scalar>, int>
00213     changeOfBasisMatrix (int s)
00214     {
00215       using Teuchos::Array;
00216       using Teuchos::ArrayView;
00217       typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00218       typedef Teuchos::ScalarTraits<Scalar> STS;
00219       typedef Teuchos::ScalarTraits<magnitude_type> STM;
00220 
00221       // Only keep the first s shifts, counting multiplicities.
00222       int numShifts = 0, numUniqueShifts = 0;
00223       Array<int> mults;
00224       while (numShifts < s)
00225   {
00226     const int k = numUniqueShifts;
00227 
00228     // If the current shift is the first member of a complex
00229     // conjugate pair (i.e., has positive imaginary part), see
00230     // if there is room to add all multiplicities for both
00231     // members of the complex conjugate pair.
00232     if (imagParts_[k] > STM::zero())
00233       { 
00234         int curMult = 1;
00235         // Can we add at least two more shifts (counting multiplicities)?
00236         if (numShifts + 1 >= s)
00237     {
00238       // We either have to increase the room so we can add
00239       // both, or not add either member of the pair.  We
00240       // decide what to do based on how many shifts we've
00241       // added thus far: if <= 1, we add the pair and
00242       // increase s accordingly, else we don't add the
00243       // pair and stop.
00244       if (numShifts > 1)
00245         break; // Don't add any more shifts
00246     }
00247         else 
00248     { // We can add both members of the pair, with
00249       // multiplicity at least 1.  Let's see how many
00250       // pairs we can add.
00251       const int leftover = numShifts + 2*multiplicities_[k] - s;
00252       int curMult;
00253       if (leftover > 0)
00254         // If even, accept, else round up one.
00255         curMult = (leftover % 2 == 0) ? (leftover/2) : (leftover/2 + 1);
00256       else // We can add both shifts with all their multiplicities.
00257         curMult = multiplicities_[k];
00258     }
00259         numUniqueShifts += 2;
00260         numShifts += 2*curMult;
00261         mults.push_back (curMult);
00262         mults.push_back (curMult);
00263       }
00264     else
00265       {
00266         const int curMult = 
00267     (numShifts + multiplicities_[k] > s) ? 
00268     (numShifts + multiplicities_[k] - s) : 
00269     multiplicities_[k];
00270         numUniqueShifts++;
00271         numShifts += curMult;
00272         mults.push_back (curMult);
00273       }
00274   }
00275       // Views of the shifts and their multiplicities of interest.
00276       ArrayView<const magnitude_type> realParts = realParts_.view (0, numUniqueShifts);
00277       ArrayView<const magnitude_type> imagParts = imagParts_.view (0, numUniqueShifts);
00278 
00279       // FIXME (mfh 09 Feb 2011) We must also recompute the matrix
00280       // when the shifts have changed.
00281       if (B_.is_null() || 
00282     (B_->numRows() != s+1 || B_->numCols() != s) ||
00283     shiftsChanged_)
00284   B_ = makeChangeOfBasisMatrix (realParts, imagParts, mults, modified_);
00285       return B_;
00286     }
00287 
00288 
00289     void
00290     updateBasis (const Teuchos::ArrayView<const magnitude_type>& realParts,
00291      const Teuchos::ArrayView<const magnitude_type>& imagParts,
00292      const Teuchos::ArrayView<const int>& multiplicities,
00293      const int numValues)
00294     {
00295       using Teuchos::Array;
00296       using Teuchos::ArrayView;
00297       typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00298       typedef Teuchos::ScalarTraits<Scalar> STS;
00299       typedef Teuchos::ScalarTraits<magnitude_type> STM;
00300 
00301       checkModifiedNewtonShifts (realParts, imagParts, multiplicities);
00302 
00303       // Only take in new shifts if there are more new shifts than old
00304       // ones.
00305       if (numValues < realParts_.size())
00306   return;
00307       else
00308   {
00309     realParts_.resize (numValues);
00310     imagParts_.resize (numValues);
00311     mults_.resize (numValues);
00312     
00313     realParts_.assign (realParts.begin(), realParts.end());
00314     imagParts_.assign (imagParts.begin(), imagParts.end());
00315     mults_.assign (multiplicities.begin(), multiplicities.end());
00316 
00317     // Recompute the change-of-basis matrix, since the shifts
00318     // have changed.
00319     typedef Teuchos::ScalarTraits<Scalar> STS;
00320     if (B_.is_null() || B_->numRows() != numValues+1 || 
00321         B_->numCols() != numValues)
00322       B_ = makeNewtonChangeOfBasisMatrix<Scalar, STS::isComplex> (realParts_, imagParts_, mults_, modified_);
00323   }
00324     }
00325 
00326   private:
00327 
00329     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,Scalar> > B_;
00330 
00332     Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> realParts_;
00336     Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> imagParts_;
00337 
00339     bool modified_;
00340   };
00341 
00342 } // namespace Belos
00343 
00344 #endif // __Belos_NewtonOpAkx_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines