Homec4science
Diffusion tamaas (petsc) Browse

tamaas/src/core/conic.hh
c23abaf90bd2petsc

Last Change
Raw File

conic.hh
View Options

/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: 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 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef CONIC_HH
#define CONIC_HH
/* -------------------------------------------------------------------------- */
#include "static_types.hh"
#include <cmath>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Conic algebra operations
namespace conic {
template <typename DT, typename ST, UInt m>
__device__ __host__ std::remove_cv_t<DT> det(const StaticVector<DT, ST, m>& v) {
static_assert(m >= 1, "degenerate cone");
VectorProxy<const DT, m - 1> tail(&v(1));
return v(0) * v(0) - tail.l2squared();
}
template <typename DT, typename ST, UInt m>
__device__ __host__ Matrix<std::remove_cv_t<DT>, m, m>
mat(const StaticVector<DT, ST, m>& v) {
Matrix<std::remove_cv_t<DT>, m, m> matrix;
for (UInt i = 0; i < m; ++i) {
matrix(0, m) = matrix(m, 0) = v(m);
matrix(m, m) = v(0);
}
return matrix;
}
template <typename DT, typename ST, UInt m>
__device__ __host__ Vector<std::remove_cv_t<DT>, m>
hat(const StaticVector<DT, ST, m>& v) {
Vector<std::remove_cv_t<DT>, m> res(v);
VectorProxy<std::remove_cv_t<DT>, m - 1> tail(&res(1));
tail *= -1;
return res;
}
/// theta in Nesterov-Todd scaling
template <typename DT, typename ST1, typename ST2, UInt m>
__device__ __host__ std::remove_cv_t<DT>
scaling_factor(const StaticVector<DT, ST1, m>& x,
const StaticVector<DT, ST2, m>& s) {
return std::pow(det(s) / det(x), 0.25);
}
/// w in Nesterov-Todd scaling
template <typename DT, typename ST1, typename ST2, UInt m>
__device__ __host__ std::remove_cv_t<DT>
scaling_vector(const StaticVector<DT, ST1, m>& x,
const StaticVector<DT, ST2, m>& s) {
auto theta = scaling_factor(x, s);
auto res = (1. / theta) * s + theta * hat(x);
auto scal = x.dot(s);
auto z = std::sqrt(det(s) * det(x));
res *= 1 / std::sqrt(2 * (scal + z));
return res;
}
/// F in Nesterov-Todd scaling
template <typename DT, typename ST1, typename ST2, UInt m>
__device__ __host__ std::remove_cv_t<DT>
scaling_matrix(const StaticVector<DT, ST1, m>& x,
const StaticVector<DT, ST2, m>& s) {
auto theta = scaling_factor(x, s);
auto w = scaling_vector(x, s);
Matrix<std::remove_cv_t<DT>, m, m> F;
for (UInt i = 1; i < m; ++i) {
F(0, i) = F(i, 0) = w(i);
for (UInt j = 1; j < m; ++i)
F(i, j) = 1 + (w(i) * w(j)) / (1 + w(0));
}
F *= theta;
return F;
}
/// F² in Nesterov-Todd scaling
template <typename DT, typename ST1, typename ST2, UInt m>
__device__ __host__ std::remove_cv_t<DT>
scaling_matrix_squared(const StaticVector<DT, ST1, m>& x,
const StaticVector<DT, ST2, m>& s) {
auto theta = scaling_factor(x, s);
auto w = scaling_factor(x, s);
Matrix<std::remove_cv_t<DT>, m, m> F;
for (UInt i = 0; i < m; ++i)
for (UInt j = 0; j < m; ++j)
F(i, j) = (i == j) + 2 * w(i) * w(j);
F(0, 0) -= 2;
F *= theta * theta;
return F;
}
} // namespace conic
} // namespace tamaas
#endif // CONIC_HH
c4science · Help