53 #ifndef ROL_POISSONINVERSION_HPP 54 #define ROL_POISSONINVERSION_HPP 60 #include "Teuchos_LAPACK.hpp" 74 typedef typename vector::size_type
uint;
89 using Teuchos::dyn_cast;
94 using Teuchos::dyn_cast;
112 for (
uint i = 0; i <
nz_; i++) {
114 val +=
alpha_/2.0 *
hz_ * (*zp)[i]*(*zp)[i];
121 val +=
alpha_ * std::sqrt(std::pow((*zp)[i]-(*zp)[i+1],2.0)+
eps_);
139 for (
uint i = 0; i <
nz_; i++) {
140 (*gp)[i] =
alpha_ *
hz_ * (*zp)[i]/std::sqrt(std::pow((*zp)[i],2.0)+
eps_);
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_);
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);
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<Real>() ) {
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++) {
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);
434 SV p( Teuchos::rcp(
new std::vector<Real>(
nu_,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.
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)
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)
virtual Teuchos::RCP< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
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