53 #ifndef ROL_POISSONINVERSION_HPP 54 #define ROL_POISSONINVERSION_HPP 60 #include "Teuchos_LAPACK.hpp" 74 typedef typename vector::size_type
uint;
88 Teuchos::RCP<const vector>
getVector(
const V& x ) {
89 using Teuchos::dyn_cast;
94 using Teuchos::dyn_cast;
102 : nu_(nz-1), nz_(nz), hu_(1./((Real)nu_+1.)), hz_(hu_),
103 alpha_(alpha), eps_(1.e-4), reg_type_(2) {}
112 for (uint i = 0; i <
nz_; i++) {
113 if ( reg_type_ == 2 ) {
114 val += alpha_/2.0 * hz_ * (*zp)[i]*(*zp)[i];
116 else if ( reg_type_ == 1 ) {
117 val += alpha_ * hz_ * std::sqrt((*zp)[i]*(*zp)[i] + eps_);
119 else if ( reg_type_ == 0 ) {
121 val += alpha_ * std::sqrt(std::pow((*zp)[i]-(*zp)[i+1],2.0)+eps_);
131 if ( reg_type_ == 2 ) {
135 else if ( reg_type_ == 1 ) {
139 for (uint i = 0; i <
nz_; i++) {
140 (*gp)[i] = alpha_ * hz_ * (*zp)[i]/std::sqrt(std::pow((*zp)[i],2.0)+eps_);
143 else if ( reg_type_ == 0 ) {
148 for (uint i = 0; i <
nz_; i++) {
150 diff = (*zp)[i]-(*zp)[i+1];
151 (*gp)[i] = alpha_ * diff/std::sqrt(std::pow(diff,2.0)+eps_);
153 else if ( i == nz_-1 ) {
154 diff = (*zp)[i-1]-(*zp)[i];
155 (*gp)[i] = -alpha_ * diff/std::sqrt(std::pow(diff,2.0)+eps_);
158 diff = (*zp)[i]-(*zp)[i+1];
159 (*gp)[i] = alpha_ * diff/std::sqrt(std::pow(diff,2.0)+eps_);
160 diff = (*zp)[i-1]-(*zp)[i];
161 (*gp)[i] -= alpha_ * diff/std::sqrt(std::pow(diff,2.0)+eps_);
171 if ( reg_type_ == 2 ) {
173 hv.
scale(alpha_*hz_);
175 else if ( reg_type_ == 1 ) {
180 for (uint i = 0; i <
nz_; i++) {
181 (*hvp)[i] = alpha_*hz_*(*vp)[i]*eps_/std::pow(std::pow((*zp)[i],2.0)+eps_,3.0/2.0);
184 else if ( reg_type_ == 0 ) {
191 for (uint i = 0; i <
nz_; i++) {
193 diff1 = (*zp)[i]-(*zp)[i+1];
194 diff1 = eps_/std::pow(std::pow(diff1,2.0)+eps_,3.0/2.0);
195 (*hvp)[i] = alpha_* ((*vp)[i]*diff1 - (*vp)[i+1]*diff1);
197 else if ( i == nz_-1 ) {
198 diff2 = (*zp)[i-1]-(*zp)[i];
199 diff2 = eps_/std::pow(std::pow(diff2,2.0)+eps_,3.0/2.0);
200 (*hvp)[i] = alpha_* (-(*vp)[i-1]*diff2 + (*vp)[i]*diff2);
203 diff1 = (*zp)[i]-(*zp)[i+1];
204 diff1 = eps_/std::pow(std::pow(diff1,2.0)+eps_,3.0/2.0);
205 diff2 = (*zp)[i-1]-(*zp)[i];
206 diff2 = eps_/std::pow(std::pow(diff2,2.0)+eps_,3.0/2.0);
207 (*hvp)[i] = alpha_* (-(*vp)[i-1]*diff2 + (*vp)[i]*(diff1 + diff2) - (*vp)[i+1]*diff1);
220 for (uint i = 0; i <
nu_; i++) {
222 (*Mfp)[i] = hu_/6.0*(4.0*(*fp)[i] + (*fp)[i+1]);
224 else if ( i == nu_-1 ) {
225 (*Mfp)[i] = hu_/6.0*((*fp)[i-1] + 4.0*(*fp)[i]);
228 (*Mfp)[i] = hu_/6.0*((*fp)[i-1] + 4.0*(*fp)[i] + (*fp)[i+1]);
243 for ( uint i = 0; i <
nu_; i++ ) {
244 d[i] = (std::exp((*zp)[i]) + std::exp((*zp)[i+1]))/hu_;
246 o[i] *= -std::exp((*zp)[i+1])/
hu_;
251 Teuchos::LAPACK<int,Real> lp;
255 lp.PTTRF(nu_,&d[0],&o[0],&info);
256 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&(*bp)[0],ldb,&info);
274 for (uint i = 0; i <
nu_; i++) {
276 (*Bdp)[i] = 1.0/hu_*( std::exp((*zp)[i])*(*up)[i]*(*dp)[i]
277 + std::exp((*zp)[i+1])*((*up)[i]-(*up)[i+1])*(*dp)[i+1] );
279 else if ( i == nu_-1 ) {
280 (*Bdp)[i] = 1.0/hu_*( std::exp((*zp)[i])*((*up)[i]-(*up)[i-1])*(*dp)[i]
281 + std::exp((*zp)[i+1])*(*up)[i]*(*dp)[i+1] );
284 (*Bdp)[i] = 1.0/hu_*( std::exp((*zp)[i])*((*up)[i]-(*up)[i-1])*(*dp)[i]
285 + std::exp((*zp)[i+1])*((*up)[i]-(*up)[i+1])*(*dp)[i+1] );
299 for (uint i = 0; i <
nz_; i++) {
301 (*Bdp)[i] = std::exp((*zp)[i])/hu_*(*up)[i]*(*dp)[i];
303 else if ( i == nz_-1 ) {
304 (*Bdp)[i] = std::exp((*zp)[i])/hu_*(*up)[i-1]*(*dp)[i-1];
307 (*Bdp)[i] = std::exp((*zp)[i])/hu_*( ((*up)[i]-(*up)[i-1])*((*dp)[i]-(*dp)[i-1]) );
321 for (uint i = 0; i <
nz_; i++) {
323 (*Bdp)[i] = (*vp)[i]*std::exp((*zp)[i])/hu_*(*up)[i]*(*dp)[i];
325 else if ( i == nz_-1 ) {
326 (*Bdp)[i] = (*vp)[i]*std::exp((*zp)[i])/hu_*(*up)[i-1]*(*dp)[i-1];
329 (*Bdp)[i] = (*vp)[i]*std::exp((*zp)[i])/hu_*( ((*up)[i]-(*up)[i-1])*((*dp)[i]-(*dp)[i-1]) );
343 RCP<vector> bp = rcp(
new vector(nu_, 0.0) );
344 for ( uint i = 0; i <
nu_; i++ ) {
345 if ( (Real)(i+1)*hu_ < 0.5 ) {
346 (*bp)[i] = 2.0*k1*
hu_;
348 else if ( std::abs((Real)(i+1)*hu_ - 0.5) <
ROL_EPSILON ) {
349 (*bp)[i] = (k1+k2)*hu_;
351 else if ( (Real)(i+1)*hu_ > 0.5 ) {
352 (*bp)[i] = 2.0*k2*
hu_;
367 RCP<vector> rp = rcp(
new vector(nu_,0.0) );
370 for ( uint i = 0; i <
nu_; i++) {
373 StdVector<Real> Mres( Teuchos::rcp(
new std::vector<Real>(nu_,0.0) ) );
382 SV b( rcp(
new vector(nu_,0.0) ) );
392 SV res( rcp(
new vector(nu_,0.0) ) );
394 SV res1( rcp(
new vector(nu_,0.0) ) );
407 RCP<vector> up = rcp(
new vector(nu_,0.0) );
413 RCP<vector> rp = rcp(
new vector(nu_,0.0) );
416 for ( uint i = 0; i <
nu_; i++) {
420 RCP<V> Mres = res.
clone();
422 return 0.5*Mres->dot(res) +
reg_value(z);
430 SV u( rcp(
new vector(nu_,0.0) ) );
434 SV p( Teuchos::rcp(
new std::vector<Real>(nu_,0.0) ) );
441 SV g_reg( rcp(
new vector(nz_,0.0) ) );
453 SV u( rcp(
new vector(nu_,0.0) ) );
457 SV p( rcp(
new vector(nu_,0.0) ) );
461 SV w( rcp(
new vector(nu_,0.0) ) );
465 SV q( rcp(
new vector(nu_,0.0) ) );
472 SV tmp( rcp(
new vector(nz_,0.0) ) );
482 SV hv_reg( rcp(
new vector(nz_,0.0) ) );
500 int dim =
static_cast<int>(vp.size());
504 Teuchos::SerialDenseMatrix<int, Real> H(dim, dim);
506 H = computeDenseHessian<Real>(obj, x);
509 std::vector<vector> eigenvals = computeEigenvalues<Real>(H);
510 std::sort((eigenvals[0]).begin(), (eigenvals[0]).end());
513 Real inertia = (eigenvals[0])[0];
514 Real correction = 0.0;
515 if ( inertia <= 0.0 ) {
516 correction = 2.0*std::abs(inertia);
517 if ( inertia == 0.0 ) {
519 while ( eigenvals[0][cnt] == 0.0 ) {
522 correction = 0.5*eigenvals[0][cnt];
523 if ( cnt == dim-1 ) {
527 for (
int i=0; i<dim; i++) {
528 H(i,i) += correction;
533 Teuchos::SerialDenseMatrix<int, Real> invH = computeInverse<Real>(H);
536 Teuchos::SerialDenseVector<int, Real> hv_teuchos(Teuchos::View, &((*hvp)[0]), dim);
537 const Teuchos::SerialDenseVector<int, Real> v_teuchos(Teuchos::View, &(vp[0]), dim);
538 hv_teuchos.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, invH, v_teuchos, 0.0);
551 Teuchos::RCP<std::vector<Real> > x0p = Teuchos::rcp(
new std::vector<Real>(n,0.0));
552 for (
int i = 0; i < n; i++ ) {
558 Teuchos::RCP<std::vector<Real> > xp = Teuchos::rcp(
new std::vector<Real>(n,0.0));
559 Real h = 1.0/((Real)n+1), pt = 0.0, k1 = 1.0, k2 = 2.0;
560 for(
int i = 0; i < n; i++ ) {
562 if ( pt >= 0.0 && pt < 0.5 ) {
563 (*xp)[i] = std::log(k1);
565 else if ( pt >= 0.5 && pt < 1.0 ) {
566 (*xp)[i] = std::log(k2);
Provides the interface to evaluate objective functions.
void apply_transposed_linearized_control_operator(Vector< Real > &Bd, const Vector< Real > &z, const Vector< Real > &d, const Vector< Real > &u)
virtual void scale(const Real alpha)=0
Compute where .
void axpy(const Real alpha, const Vector< Real > &x)
Compute where .
void solve_state_sensitivity_equation(Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z)
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
void gradient(Vector< Real > &g, const Vector< Real > &z, Real &tol)
Compute gradient.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual void zero()
Set to zero vector.
Contains definitions for helper functions in ROL.
Defines the linear algebra or vector space interface.
void solve_poisson(Vector< Real > &u, const Vector< Real > &z, Vector< Real > &b)
void solve_adjoint_sensitivity_equation(Vector< Real > &q, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &u, const Vector< Real > &z)
void apply_transposed_linearized_control_operator_2(Vector< Real > &Bd, const Vector< Real > &z, const Vector< Real > &v, const Vector< Real > &d, const Vector< Real > &u)
Real evaluate_target(Real x)
void invHessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply inverse Hessian approximation to vector.
void solve_state_equation(Vector< Real > &u, const Vector< Real > &z)
void reg_hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &z)
void getPoissonInversion(Teuchos::RCP< Objective< Real > > &obj, Teuchos::RCP< Vector< Real > > &x0, Teuchos::RCP< Vector< Real > > &x)
void apply_mass(Vector< Real > &Mf, const Vector< Real > &f)
void solve_adjoint_equation(Vector< Real > &p, const Vector< Real > &u, const Vector< Real > &z)
virtual Teuchos::RCP< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Teuchos::RCP< const vector > getVector(const V &x)
void apply_linearized_control_operator(Vector< Real > &Bd, const Vector< Real > &z, const Vector< Real > &d, const Vector< Real > &u)
Real value(const Vector< Real > &z, Real &tol)
Compute value.
Objective_PoissonInversion(uint nz=32, Real alpha=1.e-4)
virtual void set(const Vector &x)
Set where .
Real reg_value(const Vector< Real > &z)
void reg_gradient(Vector< Real > &g, const Vector< Real > &z)
Teuchos::RCP< vector > getVector(V &x)
Poisson material inversion.
std::vector< Real > vector
static const double ROL_EPSILON
Platform-dependent machine epsilon.