110 lines
3.5 KiB
C++
110 lines
3.5 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2023 Yi Zhang (yizhang-geo@zju.edu.cn)
|
|
*
|
|
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
|
* it and/or modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation, either version 2
|
|
* of the License, or (at your option) any later version. You should have
|
|
* received a copy of the GNU Lesser General Public License along with this
|
|
* program. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
|
* the GCTL, please consider the option to obtain a commercial license for a
|
|
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
|
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
|
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
|
* to include some description of your company and the realm of its activities.
|
|
* Also add information on how to contact you by electronic and paper mail.
|
|
******************************************************/
|
|
|
|
#include "cmn_grad.h"
|
|
|
|
gctl::common_gradient::common_gradient()
|
|
{
|
|
set_lcg_message(LCG_THROW);
|
|
}
|
|
|
|
gctl::common_gradient::common_gradient(size_t Ln, size_t Mn)
|
|
{
|
|
init(Ln, Mn);
|
|
}
|
|
|
|
gctl::common_gradient::~common_gradient(){}
|
|
|
|
void gctl::common_gradient::LCG_Ax(const array<double> &x, array<double> &ax)
|
|
{
|
|
matvec(t_, G_, x, NoTrans);
|
|
matvec(ax, G_, t_, Trans);
|
|
return;
|
|
}
|
|
|
|
void gctl::common_gradient::set_solver(const lcg_para ¶)
|
|
{
|
|
set_lcg_para(para);
|
|
return;
|
|
}
|
|
|
|
void gctl::common_gradient::init(size_t Ln, size_t Mn)
|
|
{
|
|
Ln_ = Ln;
|
|
Mn_ = Mn;
|
|
g_.resize(Mn_);
|
|
B_.resize(Mn_);
|
|
G_.resize(Ln_, Mn_);
|
|
t_.resize(Ln_);
|
|
gm_.resize(Ln_);
|
|
x_.resize(Ln_, 1.0);
|
|
filled_.resize(Ln_, false);
|
|
return;
|
|
}
|
|
|
|
void gctl::common_gradient::fill_model_gradient(size_t id, const _1d_array &g)
|
|
{
|
|
if (id >= Ln_) throw std::runtime_error("[gctl::common_gradient] Invalid index.");
|
|
if (g.size() != Mn_) throw std::runtime_error("[gctl::common_gradient] Invalid array size.");
|
|
|
|
G_.fill_row(id, g);
|
|
filled_[id] = true;
|
|
gm_[id] = g.module();
|
|
return;
|
|
}
|
|
|
|
const gctl::_1d_array &gctl::common_gradient::get_common_gradient(bool normalized)
|
|
{
|
|
for (size_t i = 0; i < Ln_; i++)
|
|
{
|
|
if (!filled_[i]) throw std::runtime_error("[gctl::common_gradient] Unfilled model gradient.");
|
|
}
|
|
filled_.assign(false);
|
|
|
|
G_.normalize(RowMajor);
|
|
matvec(B_, G_, x_, Trans);
|
|
|
|
g_.random_float(-1.0, 1.0, RdUniform);
|
|
g_.normalize();
|
|
LCG_Minimize(g_, B_, LCG_CG);
|
|
|
|
if (normalized) g_.normalize();
|
|
else
|
|
{
|
|
g_.normalize();
|
|
matvec(t_, G_, g_, NoTrans);
|
|
|
|
double tmod = t_.module(L1);
|
|
double gmod = 0.0;
|
|
for (size_t i = 0; i < Ln_; i++)
|
|
{
|
|
gmod += gm_[i]*abs(t_[i])/tmod;
|
|
}
|
|
g_.scale(gmod);
|
|
}
|
|
return g_;
|
|
} |