|
Tpetra Matrix/Vector Services
Version of the Day
|
|
Tpetra Matrix/Vector Services
Version of the Day
|
00001 /* 00002 // @HEADER 00003 // *********************************************************************** 00004 // 00005 // Tpetra: Templated Linear Algebra Services Package 00006 // Copyright (2008) Sandia Corporation 00007 // 00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, 00009 // the U.S. Government retains certain rights in this software. 00010 // 00011 // Redistribution and use in source and binary forms, with or without 00012 // modification, are permitted provided that the following conditions are 00013 // met: 00014 // 00015 // 1. Redistributions of source code must retain the above copyright 00016 // notice, this list of conditions and the following disclaimer. 00017 // 00018 // 2. Redistributions in binary form must reproduce the above copyright 00019 // notice, this list of conditions and the following disclaimer in the 00020 // documentation and/or other materials provided with the distribution. 00021 // 00022 // 3. Neither the name of the Corporation nor the names of the 00023 // contributors may be used to endorse or promote products derived from 00024 // this software without specific prior written permission. 00025 // 00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00037 // 00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 00039 // 00040 // ************************************************************************ 00041 // @HEADER 00042 */ 00043 00044 #ifndef MULTIPRECCG_HPP_ 00045 #define MULTIPRECCG_HPP_ 00046 00047 #include <Teuchos_TimeMonitor.hpp> 00048 #include <Teuchos_TypeNameTraits.hpp> 00049 #include <Teuchos_ParameterList.hpp> 00050 #include <Teuchos_XMLParameterListHelpers.hpp> 00051 #include <Teuchos_FancyOStream.hpp> 00052 00053 #include <Tpetra_CrsMatrix.hpp> 00054 #include <Tpetra_Vector.hpp> 00055 #include <Tpetra_RTI.hpp> 00056 #include <Tpetra_MatrixIO.hpp> 00057 00058 #include <iostream> 00059 #include <functional> 00060 00061 #ifdef HAVE_TPETRA_QD 00062 # include <qd/qd_real.h> 00063 #endif 00064 00065 #define XFORM1 TPETRA_UNARY_TRANSFORM 00066 #define XFORM2 TPETRA_BINARY_TRANSFORM 00067 #define XFORM3 TPETRA_TERTIARY_TRANSFORM 00068 00069 #define REDUCE1(in, gexp) \ 00070 Tpetra::RTI::reduce( *in, \ 00071 Tpetra::RTI::reductionGlob<Tpetra::RTI::ZeroOp<decltype((in)->meanValue())>>( \ 00072 [=]( decltype((in)->meanValue()) in ){ return gexp; }, \ 00073 std::plus<decltype((in)->meanValue())>() ) ) 00074 00075 #define REDUCE2(in1, in2, gexp) \ 00076 Tpetra::RTI::reduce( *in1, *in2, \ 00077 Tpetra::RTI::reductionGlob<Tpetra::RTI::ZeroOp<decltype((in1)->meanValue())>>( \ 00078 [=]( decltype((in1)->meanValue()) in1, \ 00079 decltype((in2)->meanValue()) in2 ) { return gexp; },\ 00080 std::plus<decltype((in1)->meanValue())>() ) ) 00081 00082 #define XFORM2RED(out,in,texp,gexp) \ 00083 Tpetra::RTI::binary_pre_transform_reduce( *out, *in, \ 00084 Tpetra::RTI::reductionGlob<Tpetra::RTI::ZeroOp<decltype((out)->meanValue())>>( \ 00085 [=]( decltype((out)->meanValue()) out, \ 00086 decltype((in)->meanValue()) in ) { return texp; }, \ 00087 [=]( decltype((out)->meanValue()) out, \ 00088 decltype((in)->meanValue()) in ) { return gexp; }, \ 00089 std::plus<decltype((out)->meanValue())>() ) ) 00090 00091 #define XFORM3RED(out,in2,in3, texp, gexp) \ 00092 Tpetra::RTI::tertiary_pre_transform_reduce( *out, *in2, *in3, \ 00093 Tpetra::RTI::reductionGlob<Tpetra::RTI::ZeroOp<decltype((out)->meanValue())>>( \ 00094 [=]( decltype((out)->meanValue()) out, \ 00095 decltype((in2)->meanValue()) in2, \ 00096 decltype((in3)->meanValue()) in3 ) \ 00097 { return texp; }, \ 00098 [=]( decltype((out)->meanValue()) out, \ 00099 decltype((in2)->meanValue()) in2, \ 00100 decltype((in3)->meanValue()) in3 ) \ 00101 { return gexp; }, \ 00102 std::plus<decltype((out)->meanValue())>() ) ) 00103 00104 namespace TpetraExamples { 00105 00106 using Teuchos::RCP; 00107 using Tpetra::createVector; 00108 using Teuchos::ParameterList; 00109 00111 template <class Tout, class Tin, class LO, class GO, class Node> 00112 struct convertHelp { 00113 static RCP<const Tpetra::CrsMatrix<Tout,LO,GO,Node>> doit(const RCP<const Tpetra::CrsMatrix<Tin,LO,GO,Node>> &A) 00114 { 00115 return A->template convert<Tout>(); 00116 } 00117 }; 00118 00120 template <class T, class LO, class GO, class Node> 00121 struct convertHelp<T,T,LO,GO,Node> { 00122 static RCP<const Tpetra::CrsMatrix<T,LO,GO,Node>> doit(const RCP<const Tpetra::CrsMatrix<T,LO,GO,Node>> &A) 00123 { 00124 return A; 00125 } 00126 }; 00127 00128 00129 /*************************************************** 00130 * Implicit Riemannian Trust-Region Eigensolver 00131 * Baker, Absil, Gallivan 2008 00132 ***************************************************/ 00133 00134 template <class S, class Souter, class LO, class GO, class Node> 00135 S diagPrecond(const RCP<const Tpetra::Vector<S,LO,GO,Node> > &D, 00136 const RCP<const Tpetra::Vector<S,LO,GO,Node> > &x, 00137 const RCP<const Tpetra::Vector<S,LO,GO,Node> > &iDx, 00138 const S &xiDx, 00139 const RCP<const Tpetra::Vector<S,LO,GO,Node> > &r, 00140 const RCP< Tpetra::Vector<S,LO,GO,Node> > &z) 00141 { 00142 /* 00143 Solve 00144 P D P z = r 00145 for solution 00146 r = P_{u,x} inv(D) r 00147 = (I - u inv(x'*u) x') inv(D) r 00148 = z - u inv(x'*u) x'*z 00149 for u = inv(D) x = iDx 00150 z = inv(D) r 00151 00152 for our purposes, 00153 1) z = r 00154 2) fuse z \= D 00155 and xz = x'*z 00156 3) fuse z = z - iDx ( xt / xiDx ) 00157 and return z'*r 00158 */ 00159 XFORM2(z,r, r); 00160 const S alpha = XFORM3RED(z,D,x, z/D, x*z ) / xiDx; 00161 return XFORM3RED(z,iDx,r, z - iDx*alpha, z*r ); 00162 } 00163 00164 template <class S, class LO, class GO, class Node> 00165 struct tCG_workspace { 00166 RCP<Tpetra::Vector<S,LO,GO,Node>> r, z, iDx, d, Hd; 00167 }; 00168 00169 /*************************************************** 00170 * Trust-region Model Subproblem CG Solver 00171 * Given x, compute eta that minimizes the 00172 * quadratic Taylor expansion of f around x, 00173 * subject to rho_x(eta) >= rho_prime, 00174 * or, equiv., eta'*eta <= ONE/rho_prime - ONE 00175 ***************************************************/ 00177 template <class S, class Souter, class LO, class GO, class Node> 00178 void truncateCG(const RCP<Teuchos::FancyOStream> &out, int verbose, 00179 const RCP<const Tpetra::Operator<S,LO,GO,Node> > &A, 00180 const RCP<const Tpetra::Vector<S,LO,GO,Node> > &diag, 00181 const RCP<const Tpetra::Vector<Souter,LO,GO,Node> > &x_outer, 00182 const RCP<const Tpetra::Vector<Souter,LO,GO,Node> > &grad, 00183 const Souter &lx_outer, 00184 const RCP<Tpetra::Vector<S,LO,GO,Node> > &eta, 00185 const tCG_workspace<S,LO,GO,Node> &workspace, 00186 double rho_prime, double kappa, double theta, int maxItersCG, double outertol) 00187 { 00188 using Teuchos::as; 00189 using Teuchos::NO_TRANS; 00190 typedef Teuchos::ScalarTraits<S> ST; 00191 typedef Teuchos::ScalarTraits<S> STO; 00192 Teuchos::OSTab tab(out,2); 00193 const S ZERO = ST::zero(); 00194 const S ONE = ST::one(); 00195 const S D2 = ONE/rho_prime - ONE; 00196 const S lx = as<S>(lx_outer); 00197 decltype(eta) r = workspace.r, 00198 z = workspace.z, 00199 d = workspace.d, 00200 iDx = workspace.iDx, 00201 Hd = workspace.Hd; 00202 static RCP<Teuchos::Time> timer; 00203 if (timer == Teuchos::null) { 00204 std::string methName("truncatedCG<"+Teuchos::TypeNameTraits<S>::name()+">"); 00205 timer = Teuchos::TimeMonitor::getNewTimer(methName); 00206 } 00207 timer->start(); 00208 RCP<const Tpetra::Vector<S,LO,GO,Node>> x; 00209 { 00210 auto xinit = createVector<S>(r->getMap()); 00211 XFORM2(xinit,x_outer, as<S>(x_outer)); 00212 x = xinit; 00213 } 00215 // init precond. compute: 00216 // iDx = inv(D)*x 00217 // xiDx = x'*iDx 00218 const S xiDx = XFORM3RED(iDx,diag,x, x/diag, x*iDx ); 00220 // init solver 00221 XFORM1( eta, ZERO ); // eta = 0 00222 S r_r = XFORM2RED( r,grad, as<S>(grad), r*r ); // r = grad 00223 S e_e = ZERO; 00224 S z_r = diagPrecond<S,Souter,LO,GO,Node>(diag,x,iDx,xiDx,r,z); 00225 XFORM2( d,z, -z ); // d = -z 00226 const S rnrm0 = ST::squareroot(r_r); 00227 const S kconv = rnrm0*kappa; 00228 const S tconv = ST::pow(rnrm0, (theta+1)); 00229 const S futile_tol = TEUCHOS_MAX(outertol, 100*rnrm0*ST::eps()); 00230 const S convtol = TEUCHOS_MAX(TEUCHOS_MIN(kconv,tconv), futile_tol); 00231 std::string convstr; 00232 if (kconv < tconv) convstr = "non-local convergence"; 00233 else convstr = "local convergence"; 00235 // iterate 00236 if (verbose > 2) *out << "starting tCG with |res| = " << rnrm0 << std::endl; 00237 if (verbose > 2) *out << "targeting tol = " << convtol << std::endl; 00238 int k = 0; 00239 std::string reason; 00240 while (true) 00241 { 00242 // apply Hessian Hd = A*d - lx*d 00243 A->apply(*d,*Hd); 00244 const S xHd = XFORM3RED(Hd,x,d, Hd - d*lx, x*Hd ); 00245 const S dHd = XFORM3RED(Hd,x,d, Hd - x*xHd, d*Hd ); 00246 if (verbose > 19) *out << "curvature: " << dHd << std::endl; 00247 00248 /* 00249 m_x(e + alpha * d) 00250 == (e + alpha * d)'*r + 1/2 * (e + alpha * d)' * H * (e + alpha * d) 00251 == e'*r + alpha*d'*r + 1/2*e'*H*e + 1/2*alpha^2*d'*H*d 00252 ---------- ------------------ 00253 derivative w.r.t. alpha, set to zero: 00254 0 = d'*r + alpha*d'*H*d 00255 so that 00256 alpha = -d'*r / d'*H*d 00257 substituting d = (gram-schmidt) -z 00258 alpha = z'*r / d'*H*d 00259 */ 00260 // compute step size 00261 const S alpha = z_r / dHd; 00262 if (verbose > 19) *out << "alpha: " << alpha << std::endl; 00263 if (verbose > 19) *out << "z_r: " << z_r << std::endl; 00264 00265 // anticipate e'*e for this step size 00266 const S e_d = REDUCE2(eta,d, eta*d), 00267 d_d = REDUCE1(d, d*d ); 00268 const S new_e_e = e_e + 2.0*alpha*e_d + alpha*alpha*d_d; 00269 00270 // check for truncation: positive curvature or implicit trust-region boundary 00271 if (dHd <= ZERO || new_e_e > D2) 00272 { 00273 // find truncated step to boundary 00274 // find alpha_prime 00275 // eta_new'*eta_new = (eta + alpha_prime*d)'*(eta + alpha_prime*d) 00276 // = e_e + 2*alpha_prime*e_d + alpha_prime^2*d_d 00277 // choose alpha_prime > 0 such that eta_new'*eta_new == D2 00278 // solve for alpha_prime 00279 // 00280 const S alpha_prime = (-e_d + ST::squareroot(e_d*e_d + d_d*(D2-e_e))) / d_d; 00281 if (verbose > 19) *out << "alpha_prime: " << alpha_prime << std::endl; 00282 XFORM2(eta,d, eta + alpha_prime*d ); 00283 if (dHd <= ZERO) reason="negative curvature"; 00284 else reason="exceeded trust-region"; 00285 // 00286 break; 00287 } 00288 00289 // compute new optimal point 00290 XFORM2(eta,d, eta + alpha*d ); 00291 e_e = new_e_e; 00292 00293 // update the residual 00294 const S x_r = XFORM3RED(r,Hd,x, r + alpha*Hd, x*r ); 00295 // re-tangentialize r 00296 r_r = XFORM2RED(r,x, r - x*x_r, r*r ); 00297 // compute norm 00298 S rnrm = ST::squareroot(r_r); 00299 if (verbose > 9) *out << "|r| = " << rnrm << std::endl; 00300 00301 // check stopping criteria 00302 if (rnrm <= convtol) {reason = convstr; break;} 00303 00304 // save the old z'*r 00305 const S zold_rold = z_r; 00306 // z = diagPrec(r) 00307 z_r = diagPrecond<S,Souter,LO,GO,Node>(diag,x,iDx,xiDx,r,z); 00308 // compute new search direction 00309 const S beta = z_r / zold_rold; 00310 XFORM2(d,z, -z + beta*d ); 00311 00312 k += 1; 00313 if (k >= maxItersCG) {reason="maximum number of inner iterations";break;} 00314 } 00315 timer->stop(); 00316 if (verbose > 1) { 00317 *out << "leaving tCG after " << k << " iterations: " << reason << std::endl; 00318 } 00319 } 00320 00322 template <class Sinner, class S, class LO, class GO, class Node> 00323 S IRTR(const RCP<Teuchos::FancyOStream> &out, ParameterList ¶ms, 00324 const RCP<const Tpetra::CrsMatrix<S,LO,GO,Node> > &A, 00325 const RCP<Tpetra::Vector<S,LO,GO,Node> > &x) 00326 { 00327 using Teuchos::NO_TRANS; 00328 Teuchos::OSTab tab(out,2); 00329 std::string methName("IRTR<"+Teuchos::TypeNameTraits<S>::name()+","+Teuchos::TypeNameTraits<Sinner>::name()+">"); 00330 typedef Tpetra::Vector<S,LO,GO,Node> Vec; 00331 typedef Tpetra::CrsMatrix<Sinner,LO,GO,Node> MatInner; 00332 typedef Teuchos::ScalarTraits<S> ST; 00333 const Sinner ONE = Teuchos::ScalarTraits<Sinner>::one(); 00334 // get objects from level database 00335 const int maxIters = params.get<int >("maxIters", 100); 00336 const int maxItersCG = params.get<int >("maxItersCG",A->getDomainMap()->getGlobalNumElements()); 00337 const int verbose = params.get<int >("verbose",0); 00338 const double tolerance = params.get<double>("tolerance"); 00339 const double kappa = params.get<double>("kappa",0.1); 00340 const double theta = params.get<double>("theta",2.0); 00341 const double rho_prime = params.get<double>("rho_prime",.1); 00342 const bool precond = params.get<int> ("precond",1); 00343 if (verbose > 10) *out << params << std::endl; 00344 if (verbose > 10) *out << "maxIters: " << maxIters << std::endl; 00345 TEUCHOS_TEST_FOR_EXCEPTION(rho_prime >= 1.0, std::runtime_error, "IRTR: rho_prime must be less than 1"); 00346 RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::getNewTimer(methName); 00347 // x, ax, gradfx, eta 00348 auto map = A->getDomainMap(); 00349 auto Ax = createVector<S>(map), 00350 r = createVector<S>(map), 00351 eta = createVector<Sinner>(map), 00352 diag = createVector<Sinner>(map); 00353 tCG_workspace<Sinner,LO,GO,Node> workspace; 00354 workspace.r = createVector<Sinner>(map); 00355 workspace.z = createVector<Sinner>(map); 00356 workspace.iDx = createVector<Sinner>(map); 00357 workspace.d = createVector<Sinner>(map); 00358 workspace.Hd = createVector<Sinner>(map); 00359 RCP<const MatInner> Ainner = convertHelp<Sinner,S,LO,GO,Node>::doit(A); 00360 if (precond) { 00361 Ainner->getLocalDiagCopy(*diag); 00362 } 00363 else { 00364 XFORM1(diag, ONE); 00365 } 00366 timer->start(); 00368 // init the solver 00369 S x_x = REDUCE1(x, x*x ); // x'*x 00370 S xnrm = ST::squareroot(x_x); 00371 XFORM1( x, x/xnrm); // normalize x 00372 A->apply(*x,*Ax); // Ax = A*x 00373 S lambda = REDUCE2(x,Ax, x*Ax ); // x'*A*x 00374 S rr = XFORM3RED(r,x,Ax, Ax - x*lambda, r*r ); // r = (I - x*x')Ax = Ax - x*lambda and compute its norm 00375 S rnrm = ST::squareroot(rr); 00376 if (verbose) { 00377 *out << "Beginning " << methName << " with λ = " << lambda << std::endl; 00378 } 00380 // iterate 00381 int k=0; 00382 while (k < maxIters && rnrm > tolerance) { 00384 // solve the trust-region subproblem for update eta 00385 truncateCG<Sinner,S,LO,GO,Node>(out,verbose,Ainner,diag, 00386 x,r,lambda, 00387 eta, workspace, 00388 rho_prime,kappa,theta,maxItersCG,tolerance); 00390 // update iterate, compute new eigenvalue 00391 x_x = XFORM2RED(x,eta, x + eta, x*x ); // x = x+eta and x'*x 00392 xnrm = ST::squareroot(x_x); 00393 XFORM1(x, x/xnrm); // normalize x 00394 A->apply(*x,*Ax); // Ax = A*x 00395 lambda = REDUCE2(x,Ax, x*Ax ); // lambda = x'*A*x 00396 rr = XFORM3RED(r,x,Ax, Ax - x*lambda, r*r ); // r = A*x - x*lambda and r'*r 00397 rnrm = ST::squareroot(rr) / ST::magnitude(lambda); 00398 k += 1; 00399 if (verbose) { 00400 *out << "After iteration " << k << ", λ = " << std::setprecision(2) << std::scientific << lambda << ", |Ax-xλ|/|λ| = " << rnrm << std::endl; 00401 } 00402 } 00403 timer->stop(); 00404 if (verbose) { 00405 *out << "Leaving " << methName << " after " << k << " iterations." << std::endl; 00406 } 00407 return lambda; 00408 } 00409 00410 } // end of namespace TpetraExamples 00411 00412 #endif // MULTIPRECCG_HPP_
1.7.6.1
