Belos Package Browser (Single Doxygen Collection) Development
BelosNewtonBasis.hpp
Go to the documentation of this file.
00001 #ifndef __Belos_NewtonBasis_hpp
00002 #define __Belos_NewtonBasis_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 namespace Belos {
00046 
00050   template<class Scalar>
00051   static std::ostream&
00052   printShift (std::ostream& out, 
00053         const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& realPart,
00054         const typename Teuchos::ScalarTraits<Scalar>::magnitudeType& imagPart)
00055   {
00056     out << "(" << realPart << ", " << imagPart << ")";
00057     return out;
00058   }
00059 
00060   template<class Scalar, class ExceptionType>
00061   void
00062   checkModifiedNewtonShifts (const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& realParts,
00063          const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& imagParts,
00064          const Teuchos::ArrayView<const int>& mults)
00065   {
00066     using Teuchos::ArrayView;
00067     typedef Teuchos::ScalarTraits<Scalar> STS;
00068     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00069     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00070     typedef Teuchos::ArrayView<const int>::size_type size_type;
00071 
00072     TEUCHOS_TEST_FOR_EXCEPTION(realParts.size() != imagParts.size(), 
00073            ExceptionType, 
00074            "Newton basis shifts: There are " << realParts.size() 
00075            << " real values, but " << imagParts.size() 
00076            << " corresponding imaginary values.  realParts[k] + "
00077            "i*imagParts[k] comprises a single complex-valued shift,"
00078            " so the realParts and imagParts arrays should always "
00079            "have the same length.";
00080     TEUCHOS_TEST_FOR_EXCEPTION(mults.size() != realParts.size() || 
00081            mults.size() != imagParts.size(),
00082            ExceptionType, 
00083            "Newton basis shifts: There are " << realParts.size()
00084            << " unique shift values, but " << mults.size()
00085            << " multiplicities.  The realParts, imagParts, and "
00086            "mults arrays should all have the same length.";
00087       int k = 0;
00088       while (k < mults.size())
00089   {
00090     if (mults[k] <= 0)
00091       {
00092         std::ostringstream os;
00093         os << "Error in shifts for the modified Newton basis: "
00094      << "Shift number " << (k+1) << " of " << mults.size()
00095      << " has nonpositive multiplicity " << mults[k] 
00096      << ".  All shifts for the Newton basis must have "
00097     "positive multiplicity.";
00098         throw ExceptionType (os.str());
00099       }
00100     else if (imagParts[k] > STM::zero())
00101       {
00102         const bool bad = k == multiplicities_.size() - 1 ||
00103     realParts[k] != realParts[k+1] ||
00104     imagParts[k] != -imagParts[k+1] ||
00105     multiplicities[k] != multiplicities_[k+1];
00106         if (bad)
00107     {
00108       std::ostringstream os;
00109       os << "Error in shifts for the modified Newton basis: ";
00110       if (multiplicities_.size() - 1)
00111         os << "The last (of " << mults.size() << ") shift " 
00112            << printShift(realParts[k], imagParts[k])
00113            << " has nonzero imaginary part.  "
00114           "This is forbidden, since complex shifts must occur in "
00115           "consecutive complex conjugate pairs in the modified "
00116           "Newton basis, but this complex shift occurs on its own.";
00117       else if (realParts[k] != realParts[k+1] ||
00118          imagParts[k] != -imagParts[k+1])
00119         os << "The current (" << (k+1) << " of " << mults.size() 
00120            << ") shift "
00121            << printShift(realParts[k], imagParts[k])
00122            << " is complex, but the shift that follows "
00123            << printShift(realParts[k+1], imagParts[k+1])
00124            << " is not its complex conjugate.  This is forbidden, "
00125           "since shifts must occur in consecutive complex "
00126           "conjugate pairs.";
00127       else if (multiplicities[k] != multiplicities_[k+1])
00128         os << "Shifts number " << (k+1) << " and " << (k+2) 
00129            << " are consecutive shifts that form a complex "
00130           "conjugate pair, but they have different "
00131           "multiplicities.  This is forbidden in the modified "
00132           "Newton basis, since complex shifts must be the "
00133           "eigenvalues of a real matrix, and thus may only occur "
00134           "in complex conjugate pairs with the same multiplicity.";
00135       throw ExceptionType (os.str());
00136     }
00137         k = k + 2;
00138       }
00139     else if (imagParts[k] < STM::zero())
00140       {
00141         std::ostringstream os;
00142         os << "Error in shifts for the modified Newton basis: "
00143      << "Shift number " << (k+1) << " of " << mults.size()
00144      << " is complex with negative imaginary part.  Complex "
00145     "shifts are allowed, but must occur in complex conjugate "
00146     "pairs, with the first member of each pair having positive "
00147     "imaginary part.  In this case, the other member of the pair "
00148     "should have occurred first.";
00149         throw ExceptionType (os.str());
00150       }
00151     else // imagParts[k] == 0
00152       ++k;
00153   }
00154     }
00155 
00175   template<class Scalar, bool isComplex>
00176   Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,Scalar> >
00177   makeNewtonChangeOfBasisMatrix (const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& realParts,
00178          const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& imagParts,
00179          const Teuchos::ArrayView<const int>& multiplicities,
00180          const bool modified);
00181 
00190   //
00195   template<class Scalar>
00196   Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,Scalar> >
00197   makeModifiedNewtonChangeOfBasisMatrix (const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& realParts,
00198            const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& imagParts,
00199            const Teuchos::ArrayView<const int>& multiplicities)
00200   {
00201     using Teuchos::RCP;
00202     typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00203     typedef Teuchos::ScalarTraits<Scalar> STS;
00204     typedef typename STS::magnitudeType magnitude_type;
00205 
00206     RCP<mat_type> B (new mat_type (s+1, s));
00207 
00208     int j = 0; // Current column of the B matrix
00209     int kk = 0;
00210     while (kk < multiplicities.size())
00211       {
00212   const magnitude_type curRealPart = realParts[kk];
00213   const magnitude_type curImagPart = imagParts[kk];
00214   const int curMult = multiplicities[kk];
00215   
00216   // If the current shift (kk) is the first element of a complex
00217   // conjugate pair, then the following must hold, else there is
00218   // an error:
00219   //
00220   // 1. This is at least one more shift left
00221   // 2. lambda(kk) == conj(lambda(kk+1))
00222   // 3. multiplicities[kk] == multiplicities[kk+1]
00223   //
00224   // We assume that the given shifts have been run through
00225   // checkModifiedNewtonShifts(), which will detect any of the possible error
00226   // conditions.
00227   if (imagParts[kk] > 0)
00228     {
00229       // Handle both elements of the complex conjugate pair
00230       // at the same time.
00231       for (int kkk = 0; kkk < mults[kk]; ++kkk)
00232         {
00233     (*B)(j, j) = curRealPart;
00234     (*B)(j+1, j) = STS::one();
00235     ++j;
00236     (*B)(j-1, j) = -(curImagPart * curImagPart);
00237     (*B)(j, j) = curRealPart;
00238     (*B)(j+1, j) = STS::one();
00239     ++j;
00240         }
00241       // We've already used lambda(kk+1) implicitly, so skip ahead.
00242       kk = kk + 2;
00243     }
00244   else if (imagParts[kk] == 0)
00245     {
00246       for (int kkk = 0; kkk < curMult; ++kkk)
00247         { // This is a real shift.  Hopefully Scalar has an operator=
00248     // that does the right thing for magnitude_type inputs.
00249     (*B)(j, j) = curRealPart;
00250     (*B)(j+1, j) = STS::one();
00251     ++j;
00252         }
00253       kk = kk + 1;
00254     }
00255   else
00256     throw std::logic_error("Should never get here!");
00257       }
00258     return B;
00259   }
00260 
00261 
00262   // Specialization for real Scalar
00263   template<class Scalar>
00264   Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,Scalar> >
00265   makeNewtonChangeOfBasisMatrix<Scalar, false> (const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& realParts,
00266             const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& imagParts,
00267             const Teuchos::ArrayView<const int>& multiplicities,
00268             const bool modified)
00269   {
00270     using Teuchos::RCP;
00271     typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00272     typedef Teuchos::ScalarTraits<Scalar> STS;
00273     typedef typename STS::magnitudeType magnitude_type;
00274 
00275     RCP<mat_type> B (new mat_type (s+1, s));
00276     if (modified)
00277       makeModifiedNewtonChangeOfBasisMatrix<Scalar> (realParts, imagParts, multiplicities);
00278     else
00279       {
00280   typedef typename Teuchos::ArrayView<const magnitude_type>::size_type size_type;
00281   int j = 0; // Current column index of B
00282   // For each (multiple) shift value
00283   for (size_type k = 0; k < multiplicities.size(); ++k)
00284     {
00285       TEUCHOS_TEST_FOR_EXCEPTION(imagParts[k] != 0, std::invalid_argument,
00286              "The (non-modified) Newton basis does not work when "
00287              "Scalar is a real type, but the shifts are complex. "
00288              " Shift number " << k << " = " 
00289              << printShift(realParts[k], imagParts[k]) << ".";
00290       // For each multiplicity of that shift value
00291       for (size_type kk = 0; kk < multiplicities[k]; ++kk)
00292         {
00293     (*B)(j, j) = realParts[k];
00294     (*B)(j+1, j) = STS::one();
00295     ++j;
00296         }
00297     }
00298       }
00299     return B;
00300   }
00301 
00302 
00303   // Specialization for complex Scalar
00304   template<class Scalar>
00305   Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,Scalar> >
00306   makeNewtonChangeOfBasisMatrix<Scalar, true> (const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& realParts,
00307                  const Teuchos::ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& imagParts,
00308                  const Teuchos::ArrayView<const int>& multiplicities,
00309                  const bool modified)
00310   {
00311     using Teuchos::RCP;
00312     typedef Teuchos::SerialDenseMatrix<int,Scalar> mat_type;
00313     typedef Teuchos::ScalarTraits<Scalar> STS;
00314     typedef typename STS::magnitudeType magnitude_type;
00315 
00316     RCP<mat_type> B (new mat_type (s+1, s));
00317 
00318     if (modified)
00319       return makeModifiedNewtonChangeOfBasisMatrix<Scalar> (realParts, imagParts, multiplicities);
00320     else
00321       {
00322   typedef typename Teuchos::ArrayView<const magnitude_type>::size_type size_type;
00323   int j = 0; // Current column index of B
00324   for (size_type k = 0; k < multiplicities.size(); ++k)
00325     for (size_type kk = 0; kk < multiplicities[k]; ++kk)
00326       {
00327         // We have to template this function on whether or not
00328         // Scalar is a complex-valued type, because the
00329         // two-argument constructor for Scalar only works
00330         // (syntactically) if Scalar is complex.
00331         (*B)(j, j) = Scalar (realParts[k], imagParts[k]);
00332         (*B)(j+1, j) = STS::one();
00333         ++j;
00334       }
00335       }
00336     return B;
00337   }
00338 
00339 } // namespace Belos
00340 
00341 #endif // __Belos_NewtonBasis_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines