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 
00057   template<class Scalar, class MV, class OP>
00058   class NewtonOpAkx : public OpAkx<Scalar, MV, OP> {
00059   public:
00060     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00061 
00062     NewtonOpAkx (const Teuchos::RCP<const OP>& A, 
00063      const Teuchos::RCP<const OP>& M_left,
00064      const Teuchos::RCP<const OP>& M_right,
00065      const bool modified) :
00066       OpAkx (A, M_left, M_right), modified_ (modified) {}
00067 
00068     void 
00069     computeBasis (const MV& q_last, MV& V_cur, const int s) 
00070     {
00071       using Teuchos::Range1D;
00072       using Teuchos::RCP;
00073       using Teuchos::rcp_const_cast;
00074       typedef MultiVecTraits<Scalar, MV> MVT;
00075       typedef Teuchos::ScalarTraits<Scalar> STS;
00076 
00077       TEST_FOR_EXCEPTION(s < 0, std::invalid_argument, 
00078        "Number of basis vectors to compute (s = " 
00079        << s << ") is invalid; it must be positive.");
00080       if (s == 0)
00081   return; // Nothing to do
00082       TEST_FOR_EXCEPTION(MVT::GetNumberVecs(V_cur) < s, std::invalid_argument,
00083        "You're asking to compute s = " << s << " basis "
00084        "vectors, but only have room for " 
00085        << MVT::GetNumberVecs(V_cur) << ".");
00086       RCP<const MV> v_prv = Teuchos::rcpFromRef (q_last);
00087       int k = 0;
00088       while (k < s)
00089   {
00090     RCP<MV> v_cur = MVT::CloneViewNonConst (V_cur, Range1D(k,k));
00091     if (imagParts_[k] == STS::zero())
00092       { // The current shift is real-valued.
00093         applyOp (v_prv, v_cur);
00094         MVT::MvAddMv (-realParts_[k], v_prv, STS::one(), v_cur, v_cur);
00095         // Shift the v_prv view over to the next column of
00096         // V_cur.  v_cur gets created on demand, so we don't
00097         // have to reassign it.
00098         v_prv = rcp_const_cast<const MV> (v_cur);
00099         ++k;
00100       }
00101     else // The current shift has nonzero imaginary part.
00102       {  
00103         // "Modified" Newton basis trick to ensure real arithmetic,
00104         // even if the shifts are complex (as long as any complex
00105         // shifts occur in consecutive complex conjugate pairs with
00106         // the positive-real-part member of the pair first).
00107         //
00108         // v_cur := (A - \Re \theta I) v_prv
00109         applyOp (v_prv, v_cur);
00110         MVT::MvAddMv (-realParts_[k], v_prv, STS::one(), v_cur, v_cur);
00111         // v_nxt := (A - \Re \theta I) v_cur + (\Im \theta)^2 v_prv
00112         RCP<MV> v_nxt = MVT::CloneViewNonConst (V_cur, Range1D(k+1, k+1));
00113         applyOp (v_nxt, v_cur);
00114         MVT::MvAddMv (-realParts_[k], v_cur, STS::one(), v_nxt, v_nxt);
00115         // FIXME (mfh 14 Feb 2011)
00116         //
00117         // You might think that we should be concerned that
00118         // squaring the imaginary part of the current shift
00119         // might overflow.  This is indeed possible, but the
00120         // nonmodified Newton basis doesn't really promise any
00121         // better (two steps more or less multiply by the square
00122         // of the magnitude).  This probably needs some more
00123         // analysis, but I don't think it's really broken.
00124         {
00125     const magnitude_type imSqu = imagParts_[k] * imagParts_[k];
00126     MVT::MvAddMv (imSqu, v_prv, STS::one(), v_nxt, v_nxt);
00127         }
00128         // Shift the v_prv view over to the next column of
00129         // V_cur.  v_cur and v_nxt get created on demand, so we
00130         // don't have to reassign them.
00131         v_prv = rcp_const_cast<const MV> (v_cur);
00132         k += 2;
00133       }
00134   }
00135     }
00136 
00137     void 
00138     computeFlexibleBasis (const MV& q_last, 
00139         MV& Z_cur, 
00140         MV& V_cur, 
00141         const int s)
00142     {
00143       using Teuchos::Range1D;
00144       using Teuchos::RCP;
00145       using Teuchos::rcp_const_cast;
00146       typedef MultiVecTraits<Scalar, MV> MVT;
00147 
00148       TEST_FOR_EXCEPTION(s < 0, std::invalid_argument, 
00149        "Number of basis vectors to compute (s = " 
00150        << s << ") is invalid; it must be positive.");
00151       if (s == 0)
00152   return; // Nothing to do
00153       TEST_FOR_EXCEPTION(MVT::GetNumberVecs(V_cur) < s, std::invalid_argument,
00154        "You're asking to compute s = " << s << " basis "
00155        "vectors, but only have room for " 
00156        << MVT::GetNumberVecs(V_cur) << " in V_cur.");
00157       TEST_FOR_EXCEPTION(MVT::GetNumberVecs(Z_cur) < s, std::invalid_argument,
00158        "You're asking to compute s = " << s << " basis "
00159        "vectors, but only have room for " 
00160        << MVT::GetNumberVecs(Z_cur) << " in Z_cur.");
00161       RCP<const MV> v_prv = Teuchos::rcpFromRef (q_last);
00162       int k = 0;
00163       while (k < s)
00164   {
00165     RCP<MV> z_cur = MVT::CloneViewNonConst (Z_cur, Range1D(k,k));
00166     RCP<MV> v_cur = MVT::CloneViewNonConst (V_cur, Range1D(k,k));
00167     applyFlexibleOp (v_prv, z_cur, v_cur);
00168     if (imagParts_[k] == STS::zero())
00169       {
00170         // We've already applied the full operator to v_prv,
00171         // producing v_cur.  Now we just need to apply the shift
00172         // to v_cur.
00173         MVT::MvAddMv (-realParts_[k], v_prv, STS::one(), v_cur, v_cur);
00174         // Shift the v_prv view over to the next column of
00175         // V_cur.  z_cur and v_cur get created on demand, so we
00176         // don't have to reassign them.
00177         v_prv = rcp_const_cast<const MV> (v_cur);
00178         ++k;
00179       }
00180     else
00181       {
00182         // Apply the first shift (realParts_[k]) to v_cur.
00183         MVT::MvAddMv (-realParts_[k], v_prv, STS::one(), v_cur, v_cur);
00184         RCP<MV> z_nxt = MVT::CloneViewNonConst (Z_cur, Range1D(k+1, k+1));
00185         RCP<MV> v_nxt = MVT::CloneViewNonConst (V_cur, Range1D(k+1, k+1));
00186         // After the next line, z_nxt will be just the result of
00187         // applying the right preconditioner to v_cur (to which
00188         // the first shift has already been applied).  Thus,
00189         // it's only necessary to apply the second shift to
00190         // v_nxt.
00191         applyFlexibleOp (v_cur, z_nxt, v_nxt);
00192         MVT::MvAddMv (-realParts_[k], v_cur, STS::one(), v_nxt, v_nxt);
00193         {
00194     const magnitude_type imSqu = imagParts_[k] * imagParts_[k];
00195     MVT::MvAddMv (imSqu, v_prv, STS::one(), v_nxt, v_nxt);
00196         }
00197         // Shift the v_prv view over to the next column of
00198         // V_cur.  z_cur, v_cur, z_nxt, and v_nxt get created on
00199         // demand, so we don't have to reassign them.
00200         v_prv = rcp_const_cast<const MV> (v_cur);
00201         k += 2;
00202       }
00203   }
00204     }
00205 
00209     std::pair<Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,Scalar>, int>
00210     changeOfBasisMatrix (int s)
00211     {
00212       using Teuchos::Array;
00213       using Teuchos::ArrayView;
00214       typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00215       typedef Teuchos::ScalarTraits<Scalar> STS;
00216       typedef Teuchos::ScalarTraits<magnitude_type> STM;
00217 
00218       // Only keep the first s shifts, counting multiplicities.
00219       int numShifts = 0, numUniqueShifts = 0;
00220       Array<int> mults;
00221       while (numShifts < s)
00222   {
00223     const int k = numUniqueShifts;
00224 
00225     // If the current shift is the first member of a complex
00226     // conjugate pair (i.e., has positive imaginary part), see
00227     // if there is room to add all multiplicities for both
00228     // members of the complex conjugate pair.
00229     if (imagParts_[k] > STM::zero())
00230       { 
00231         int curMult = 1;
00232         // Can we add at least two more shifts (counting multiplicities)?
00233         if (numShifts + 1 >= s)
00234     {
00235       // We either have to increase the room so we can add
00236       // both, or not add either member of the pair.  We
00237       // decide what to do based on how many shifts we've
00238       // added thus far: if <= 1, we add the pair and
00239       // increase s accordingly, else we don't add the
00240       // pair and stop.
00241       if (numShifts > 1)
00242         break; // Don't add any more shifts
00243     }
00244         else 
00245     { // We can add both members of the pair, with
00246       // multiplicity at least 1.  Let's see how many
00247       // pairs we can add.
00248       const int leftover = numShifts + 2*multiplicities_[k] - s;
00249       int curMult;
00250       if (leftover > 0)
00251         // If even, accept, else round up one.
00252         curMult = (leftover % 2 == 0) ? (leftover/2) : (leftover/2 + 1);
00253       else // We can add both shifts with all their multiplicities.
00254         curMult = multiplicities_[k];
00255     }
00256         numUniqueShifts += 2;
00257         numShifts += 2*curMult;
00258         mults.push_back (curMult);
00259         mults.push_back (curMult);
00260       }
00261     else
00262       {
00263         const int curMult = 
00264     (numShifts + multiplicities_[k] > s) ? 
00265     (numShifts + multiplicities_[k] - s) : 
00266     multiplicities_[k];
00267         numUniqueShifts++;
00268         numShifts += curMult;
00269         mults.push_back (curMult);
00270       }
00271   }
00272       // Views of the shifts and their multiplicities of interest.
00273       ArrayView<const magnitude_type> realParts = realParts_.view (0, numUniqueShifts);
00274       ArrayView<const magnitude_type> imagParts = imagParts_.view (0, numUniqueShifts);
00275 
00276       // FIXME (mfh 09 Feb 2011) We must also recompute the matrix
00277       // when the shifts have changed.
00278       if (B_.is_null() || 
00279     (B_->numRows() != s+1 || B_->numCols() != s) ||
00280     shiftsChanged_)
00281   B_ = makeChangeOfBasisMatrix (realParts, imagParts, mults, modified_);
00282       return B_;
00283     }
00284 
00285 
00286     void
00287     updateBasis (const Teuchos::ArrayView<const magnitude_type>& realParts,
00288      const Teuchos::ArrayView<const magnitude_type>& imagParts,
00289      const Teuchos::ArrayView<const int>& multiplicities,
00290      const int numValues)
00291     {
00292       using Teuchos::Array;
00293       using Teuchos::ArrayView;
00294       typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00295       typedef Teuchos::ScalarTraits<Scalar> STS;
00296       typedef Teuchos::ScalarTraits<magnitude_type> STM;
00297 
00298       checkModifiedNewtonShifts (realParts, imagParts, multiplicities);
00299 
00300       // Only take in new shifts if there are more new shifts than old
00301       // ones.
00302       if (numValues < realParts_.size())
00303   return;
00304       else
00305   {
00306     realParts_.resize (numValues);
00307     imagParts_.resize (numValues);
00308     mults_.resize (numValues);
00309     
00310     realParts_.assign (realParts.begin(), realParts.end());
00311     imagParts_.assign (imagParts.begin(), imagParts.end());
00312     mults_.assign (multiplicities.begin(), multiplicities.end());
00313 
00314     // Recompute the change-of-basis matrix, since the shifts
00315     // have changed.
00316     typedef Teuchos::ScalarTraits<Scalar> STS;
00317     if (B_.is_null() || B_->numRows() != numValues+1 || 
00318         B_->numCols() != numValues)
00319       B_ = makeNewtonChangeOfBasisMatrix<Scalar, STS::isComplex> (realParts_, imagParts_, mults_, modified_);
00320   }
00321     }
00322 
00323   private:
00324 
00326     Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,Scalar> > B_;
00327 
00329     Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> realParts_;
00333     Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> imagParts_;
00334 
00336     bool modified_;
00337   };
00338 
00339 } // namespace Belos
00340 
00341 #endif // __Belos_NewtonOpAkx_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines