Skip to content

GitLab

Commit ae0bf6b1 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Fixes and updates 

Including issue #15
parent 5cbb6419
Expand all Hide whitespace changes
Inline Side-by-side
Examples/ADFVDiscr/main.cpp
View file @ ae0bf6b1
......@@ -149,6 +149,7 @@ int main(int argc,char ** argv)
{
Automatizator aut(m);
Automatizator::MakeCurrent(&aut);
INMOST_DATA_ENUM_TYPE iphi = aut.RegisterDynamicTag(phi,CELL);
aut.EnumerateDynamicTags();
......@@ -164,7 +165,7 @@ int main(int argc,char ** argv)
#pragma omp parallel
#endif
{
enhanced_variable flux(aut); //should be more efficient to define here to avoid multiple memory allocations if storage for variations should be expanded
variable flux; //should be more efficient to define here to avoid multiple memory allocations if storage for variations should be expanded
#if defined(USE_OMP)
#pragma omp for
#endif
......
Examples/ADMFD/main.cpp
View file @ ae0bf6b1
......@@ -454,9 +454,10 @@ int main(int argc,char ** argv)
{ //Main loop for problem solution
Automatizator aut(m); // declare class to help manage unknowns
Automatizator::MakeCurrent(&aut);
Sparse::RowMerger & merger = aut.GetMerger(); //get structure that helps matrix assembly
dynamic_variable P(aut,aut.RegisterDynamicTag(tag_P,CELL|FACE)); //register pressure as primary unknown
enhanced_variable calc(aut); //declare variable that helps calculating the value with variations
variable calc; //declare variable that helps calculating the value with variations
aut.EnumerateDynamicTags(); //enumerate all primary variables
......
Source/Autodiff/autodiff.cpp
View file @ ae0bf6b1
......@@ -13,9 +13,40 @@
namespace INMOST
{
Automatizator * Automatizator::CurrentAutomatizator = NULL;
bool print_ad_ctor = false;
bool GetAutodiffPrint() {return print_ad_ctor;}
void SetAutodiffPrint(bool set) {print_ad_ctor = set;}
bool CheckCurrentAutomatizator() {return Automatizator::HaveCurrent();}
void multivar_expression::FromBasicExpression(const basic_expression & expr)
{
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
expr.GetDerivative(1.0,merger);
merger.RetriveRow(entries);
merger.Clear();
}
void multivar_expression::AddBasicExpression(INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr)
{
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
merger.PushRow(multme,entries);
expr.GetDerivative(multit,merger);
merger.RetriveRow(entries);
merger.Clear();
}
void multivar_expression::FromGetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
{
Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
GetDerivative(mult,merger);
merger.AddRow(1.0,r);
merger.RetriveRow(r);
merger.Clear();
}
#if defined(NEW_VERSION)
//! returns offset from the end of precomputed z
//! returns offset from the end of precomputed z
void Automatizator::DerivativeFill(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, Sparse::RowMerger & entries, INMOST_DATA_REAL_TYPE multval, void * user_data)
{
INMOST_DATA_ENUM_TYPE voffset = var.values_offset(element), doffset = var.derivatives_offset(element);
......
Source/Data/tag.cpp
View file @ ae0bf6b1
......@@ -400,6 +400,8 @@ namespace INMOST
{
if( !tag.isSparse(mask) )
{
//this was already done in Mesh::DeleteTag()
//ReallocateData(tag,ElementNum(mask),0); //this should clean up the structures
dense_data[tpos].clear(); //here all data should be deleted
empty_dense_data.push_back(tpos);
}
......@@ -517,6 +519,17 @@ namespace INMOST
}
}
}
#if defined(USE_AUTODIFF)
else if( data_type == DATA_VARIABLE ) //Have to perform in-memory deallocation to correctly remove class inheritance
{
for(INMOST_DATA_ENUM_TYPE it = new_size; it < old_size; ++it)
{
variable * p = static_cast<variable *>(static_cast<void *>(&arr[it]));
for(INMOST_DATA_ENUM_TYPE jt = 0; jt < data_size; ++jt)
p[jt].~variable();
}
}
#endif
#if defined(USE_OMP)
#pragma omp critical
#endif
......@@ -538,6 +551,17 @@ namespace INMOST
#endif
}
}
#if defined(USE_AUTODIFF)
else if( data_type == DATA_VARIABLE ) //Have to perform in-memory allocation to correctly setup class inheritance
{
for(INMOST_DATA_ENUM_TYPE it = old_size; it < new_size; ++it)
{
variable * p = static_cast<variable *>(static_cast<void *>(&arr[it]));
for(INMOST_DATA_ENUM_TYPE jt = 0; jt < data_size; ++jt)
new (p+jt) variable();
}
}
#endif
}
}
......
Source/Headers/container.hpp
View file @ ae0bf6b1
......@@ -18,18 +18,19 @@
//TODO
// 1. change to uniform size_type instead of size_t, make it INMOST_DATA_ENUM_TYPE
/*
template<class element, class T1> struct isInputRandomIterators
{
static void constraints(T1 a, T1 b) { /*element x = static_cast<element>(*a); (void)x;*/ ++a; (void)a++; a==a; a!=a; a-b; }
static void constraints(T1 a, T1 b) { ++a; (void)a++; a==a; a!=a; a-b; }
isInputRandomIterators() { void(*p)(T1,T1) = constraints; (void)p; }
};
template<class element, class T1> struct isInputForwardIterators
{
static void constraints(T1 a) { /*element x = static_cast<element>(*a); (void)x;*/ ++a; (void)a++; }
static void constraints(T1 a) {++a; (void)a++; }
isInputForwardIterators() { void(*p)(T1) = constraints; (void)p; }
};
*/
namespace INMOST
{
......
Source/Headers/inmost_autodiff.h
View file @ ae0bf6b1
......@@ -92,6 +92,7 @@ namespace INMOST
typedef dynarray<stencil_pair, 64> stencil_pairs;
typedef void(*stencil_callback)(const Storage & current_element, stencil_pairs & out_stencil, void * user_data);
private:
static Automatizator * CurrentAutomatizator;
#if defined(USE_OMP)
std::vector<Sparse::RowMerger> merger;
#else
......@@ -187,6 +188,14 @@ namespace INMOST
return merger;
#endif
}
/// Remove global current automatizator.
static void RemoveCurrent() {CurrentAutomatizator = NULL;}
/// Set current global automatizator, so that variable will be optimized with row merger.
static void MakeCurrent(Automatizator * aut) {CurrentAutomatizator = aut;}
/// Check that there is an automatizator.
static bool HaveCurrent() {return CurrentAutomatizator != NULL;}
/// Retrive the automatizator.
static Automatizator * GetCurrent() {return CurrentAutomatizator;}
};
#if defined(NEW_VERSION)
......
Source/Headers/inmost_data.h
View file @ ae0bf6b1
......@@ -240,24 +240,38 @@ namespace INMOST
typedef tag_array_type::iterator iteratorTag;
public:
virtual ~TagManager();
/// Check existance of a data tag by it's name.
bool HaveTag(std::string name) const;
/// Retrive a data tag by it's name.
Tag GetTag(std::string name) const;
/// Retrive names for all the tags present on the mesh.
void ListTagNames(std::vector<std::string> & list) const;
/// Create tag with prescribed attributes.
Tag CreateTag(Mesh * m, std::string name, DataType dtype, ElementType etype, ElementType sparse, INMOST_DATA_ENUM_TYPE size = ENUMUNDEF);
/// Delete tag from certain elements.
virtual Tag DeleteTag(Tag tag, ElementType mask);
/// Check that the tag was defined on certain elements.
bool ElementDefined(Tag const & tag, ElementType etype) const;
protected:
/// Shrink or enlarge arrays for a dense data.
void ReallocateData(const Tag & t, INMOST_DATA_INTEGER_TYPE etypenum,INMOST_DATA_ENUM_TYPE new_size);
/// Reallocate all the data for all the tags.
void ReallocateData(INMOST_DATA_INTEGER_TYPE etypenum, INMOST_DATA_ENUM_TYPE new_size);
///Retrive substructure for representation of the sparse data without permission for modification.
__INLINE sparse_sub_type const & GetSparseData(int etypenum, int local_id) const {return sparse_data[etypenum][local_id];}
///Retrive substructure for representation of the sparse data.
__INLINE sparse_sub_type & GetSparseData(int etypenum, int local_id) {return sparse_data[etypenum][local_id];}
///Retrive substructure for representation of the dense data without permission for modification.
__INLINE dense_sub_type const & GetDenseData(int pos) const {return dense_data[pos];}
///Retrive substructure for representation of the dense data.
__INLINE dense_sub_type & GetDenseData(int pos) {return dense_data[pos];}
///Copy data from one element to another.
static void CopyData(const Tag & t, void * adata, const void * bdata);
///Destroy data that represents array of variable size.
static void DestroyVariableData(const Tag & t, void * adata);
protected:
typedef tag_array_type::iterator tag_iterator;
typedef tag_array_type::const_iterator tag_const_iterator;
typedef tag_array_type::iterator tag_iterator; //< Use this type to iterate over tags of the mesh.
typedef tag_array_type::const_iterator tag_const_iterator; //< Use this type to iterate over tags of the mesh without right for modification.
protected:
tag_array_type tags;
empty_data empty_dense_data;
......@@ -432,16 +446,16 @@ namespace INMOST
/// If there is a link to handle provided (automatically by ElementArray and reference_array),
/// then remote handle value will be modified
Storage & operator = (Storage const & other);
__INLINE bool operator < (const Storage & other) const;
__INLINE bool operator > (const Storage & other) const;
__INLINE bool operator <= (const Storage & other) const;
__INLINE bool operator >= (const Storage & other) const;
__INLINE bool operator == (const Storage & other) const;
__INLINE bool operator != (const Storage & other) const;
__INLINE Storage * operator-> ();
__INLINE const Storage * operator-> () const;
__INLINE Storage & self ();
__INLINE const Storage & self () const;
__INLINE bool operator < (const Storage & other) const {return handle < other.handle;}
__INLINE bool operator > (const Storage & other) const {return handle > other.handle;}
__INLINE bool operator <= (const Storage & other) const {return handle <= other.handle;}
__INLINE bool operator >= (const Storage & other) const {return handle >= other.handle;}
__INLINE bool operator == (const Storage & other) const {return handle == other.handle;}
__INLINE bool operator != (const Storage & other) const {return handle != other.handle;}
__INLINE Storage * operator-> () {return this;}
__INLINE const Storage * operator-> () const {return this;}
__INLINE Storage & self () {return *this;}
__INLINE const Storage & self () const {return *this;}
virtual ~Storage() {}
public:
/// Retrieve real value associated with Tag. Implemented in inmost_mesh.h.
......
Source/Headers/inmost_dense.h
View file @ ae0bf6b1
......@@ -4,6 +4,7 @@
#if defined(USE_AUTODIFF)
#include "inmost_expression.h"
#endif
#include <iomanip>
// Matrix with n columns and m rows
// __m__
......@@ -41,36 +42,64 @@ namespace INMOST
enumerator n, m;
static Var sign_func(const Var & a, const Var & b) {return (b >= 0.0 ? 1.0 : -1.0)*fabs(a);}
static Var max_func(const Var & x, const Var & y) { return x > y ? x : y; }
static Var sign_func(const Var & a, const Var & b) {return (b >= 0.0 ? fabs(a) : -fabs(a));}
static INMOST_DATA_REAL_TYPE max_func(INMOST_DATA_REAL_TYPE x, INMOST_DATA_REAL_TYPE y) { return x > y ? x : y; }
static Var pythag(const Var & a, const Var & b)
{
Var at = fabs(a), bt = fabs(b), ct, result;
if (at > bt) { ct = bt / at; result = sqrt(at) * sqrt(at + ct * bt); }
else if (bt > 0.0) { ct = at / bt; result = sqrt(bt) * sqrt(bt + ct * at); }
if (at > bt) { ct = bt / at; result = at * sqrt(1.0 + ct * ct); }
else if (bt > 0.0) { ct = at / bt; result = bt * sqrt(1.0 + ct * ct); }
else result = 0.0;
return result;
}
public:
void RemoveRow(enumerator row)
{
for(enumerator k = row+1; k < n; ++k)
{
for(enumerator l = 0; l < m; ++l)
(*this)(k-1,l) = (*this)(k,l);
}
space.resize((n-1)*m);
--n;
}
void RemoveColumn(enumerator col)
{
array<Var> tmp(n*(m-1));
for(enumerator k = 0; k < n; ++k)
{
for(enumerator l = col+1; l < m; ++l)
(*this)(k,l-1) = (*this)(k,l);
}
space.swap(tmp);
--m;
}
void Swap(Matrix & b)
{
space.swap(b.space);
std::swap(n,b.n);
std::swap(m,b.m);
}
bool SVD(Matrix & U, Matrix & Sigma, Matrix & VT)
/// Singular value decomposition.
/// Reconstruct matrix: A = U*Sigma*V.Transpose().
/// @param U Left unitary matrix, U^T U = I.
/// @param Sigma Diagonal matrix with singular values.
/// @param V Right unitary matrix, not transposed.
/// @param order_singular_values
bool SVD(Matrix & U, Matrix & Sigma, Matrix & V, bool order_singular_values = true)
{
int flag, i, its, j, jj, k, l, nm;
Var c, f, h, s, x, y, z;
Var anorm = 0.0, g = 0.0, scale = 0.0;
Var g = 0.0, scale = 0.0;
INMOST_DATA_REAL_TYPE anorm = 0.0;
if (n < m)
{
bool success = Transpose().SVD(U,Sigma,VT);
bool success = Transpose().SVD(U,Sigma,V);
if( success )
{
U.Swap(VT);
U.Swap(V);
U = U.Transpose();
VT = VT.Transpose();
V = V.Transpose();
return true;
}
else return false;
......@@ -80,92 +109,92 @@ namespace INMOST
U = (*this);
Sigma.Resize(m,m);
Sigma.Zero();
VT.Resize(m,m);
V.Resize(m,m);
std::swap(n,m); //this how original algorithm takes it
// Householder reduction to bidiagonal form
for (i = 0; i < n; i++)
for (i = 0; i < (int)n; i++)
{
// left-hand reduction
l = i + 1;
rv1[i] = scale * g;
g = s = scale = 0.0;
if (i < m)
if (i < (int)m)
{
for (k = i; k < m; k++) scale += fabs(U[k][i]);
for (k = i; k < (int)m; k++) scale += fabs(U(k,i));
if (get_value(scale))
{
for (k = i; k < m; k++)
for (k = i; k < (int)m; k++)
{
U[k][i] /= scale;
s += U[k][i] * U[k][i];
U(k,i) /= scale;
s += U(k,i) * U(k,i);
}
f = U[i][i];
f = U(i,i);
g = -sign_func(sqrt(s), f);
h = f * g - s;
U[i][i] = f - g;
U(i,i) = f - g;
if (i != n - 1)
{
for (j = l; j < n; j++)
for (j = l; j < (int)n; j++)
{
for (s = 0.0, k = i; k < m; k++) s += U[k][i] * U[k][j];
for (s = 0.0, k = i; k < (int)m; k++) s += U(k,i) * U(k,j);
f = s / h;
for (k = i; k < m; k++) U[k][j] += f * U[k][i];
for (k = i; k < (int)m; k++) U(k,j) += f * U(k,i);
}
}
for (k = i; k < m; k++) U[k][i] *= scale;
for (k = i; k < (int)m; k++) U(k,i) *= scale;
}
}
Sigma[i][i] = scale * g;
Sigma(i,i) = scale * g;
// right-hand reduction
g = s = scale = 0.0;
if (i < m && i != n - 1)
if (i < (int)m && i != n - 1)
{
for (k = l; k < n; k++) scale += fabs(U[i][k]);
for (k = l; k < (int)n; k++) scale += fabs(U(i,k));
if (get_value(scale))
{
for (k = l; k < n; k++)
for (k = l; k < (int)n; k++)
{
U[i][k] = U[i][k]/scale;
s += U[i][k] * U[i][k];
U(i,k) = U(i,k)/scale;
s += U(i,k) * U(i,k);
}
f = U[i][l];
f = U(i,l);
g = -sign_func(sqrt(s), f);
h = f * g - s;
U[i][l] = f - g;
for (k = l; k < n; k++) rv1[k] = U[i][k] / h;
U(i,l) = f - g;
for (k = l; k < (int)n; k++) rv1[k] = U(i,k) / h;
if (i != m - 1)
{
for (j = l; j < m; j++)
for (j = l; j < (int)m; j++)
{
for (s = 0.0, k = l; k < n; k++) s += U[j][k] * U[i][k];
for (k = l; k < n; k++) U[j][k] += s * rv1[k];
for (s = 0.0, k = l; k < (int)n; k++) s += U(j,k) * U(i,k);
for (k = l; k < (int)n; k++) U(j,k) += s * rv1[k];
}
}
for (k = l; k < n; k++) U[i][k] *= scale;
for (k = l; k < (int)n; k++) U(i,k) *= scale;
}
}
anorm = max_func(anorm,fabs(Sigma[i][i]) + fabs(rv1[i]));
anorm = max_func(anorm,fabs(get_value(Sigma(i,i))) + fabs(get_value(rv1[i])));
}
// accumulate the right-hand transformation
for (i = n - 1; i >= 0; i--)
{
if (i < n - 1)
if (i < (int)(n - 1))
{
if (get_value(g))
{
for (j = l; j < n; j++) VT[j][i] = ((U[i][j] / U[i][l]) / g);
for (j = l; j < (int)n; j++) V(j,i) = ((U(i,j) / U(i,l)) / g);
// double division to avoid underflow
for (j = l; j < n; j++)
for (j = l; j < (int)n; j++)
{
for (s = 0.0, k = l; k < n; k++) s += U[i][k] * VT[k][j];
for (k = l; k < n; k++) VT[k][j] += s * VT[k][i];
for (s = 0.0, k = l; k < (int)n; k++) s += U(i,k) * V(k,j);
for (k = l; k < (int)n; k++) V(k,j) += s * V(k,i);
}
}
for (j = l; j < n; j++) VT[i][j] = VT[j][i] = 0.0;
for (j = l; j < (int)n; j++) V(i,j) = V(j,i) = 0.0;
}
VT[i][i] = 1.0;
V(i,i) = 1.0;
g = rv1[i];
l = i;
}
......@@ -174,26 +203,26 @@ namespace INMOST
for (i = n - 1; i >= 0; i--)
{
l = i + 1;
g = Sigma[i][i];
if (i < n - 1)
for (j = l; j < n; j++)
U[i][j] = 0.0;
g = Sigma(i,i);
if (i < (int)(n - 1))
for (j = l; j < (int)n; j++)
U(i,j) = 0.0;
if (get_value(g))
{
g = 1.0 / g;
if (i != n - 1)
{
for (j = l; j < n; j++)
for (j = l; j < (int)n; j++)
{
for (s = 0.0, k = l; k < m; k++) s += (U[k][i] * U[k][j]);
f = (s / U[i][i]) * g;
for (k = i; k < m; k++) U[k][j] += f * U[k][i];
for (s = 0.0, k = l; k < (int)m; k++) s += (U(k,i) * U(k,j));
f = (s / U(i,i)) * g;
for (k = i; k < (int)m; k++) U(k,j) += f * U(k,i);
}
}
for (j = i; j < m; j++) U[j][i] = U[j][i]*g;
for (j = i; j < (int)m; j++) U(j,i) = U(j,i)*g;
}
else for (j = i; j < m; j++) U[j][i] = 0.0;
U[i][i] += 1;
else for (j = i; j < (int)m; j++) U(j,i) = 0.0;
U(i,i) += 1;
}
// diagonalize the bidiagonal form
......@@ -205,12 +234,12 @@ namespace INMOST
for (l = k; l >= 0; l--)
{// test for splitting
nm = l - 1;
if (fabs(rv1[l]) + anorm == anorm)
if (fabs(get_value(rv1[l])) + anorm == anorm)
{
flag = 0;
break;
}
if (fabs(Sigma[nm][nm]) + anorm == anorm)
if (fabs(get_value(Sigma(nm,nm))) + anorm == anorm)
break;
}
if (flag)
......@@ -220,31 +249,31 @@ namespace INMOST
for (i = l; i <= k; i++)
{
f = s * rv1[i];
if (fabs(f) + anorm != anorm)
if (fabs(get_value(f)) + anorm != anorm)
{
g = Sigma[i][i];
g = Sigma(i,i);
h = pythag(f, g);
Sigma[i][i] = h;
Sigma(i,i) = h;
h = 1.0 / h;
c = g * h;
s = (- f * h);
for (j = 0; j < m; j++)
for (j = 0; j < (int)m; j++)
{
y = U[j][nm];
z = U[j][i];
U[j][nm] = (y * c + z * s);
U[j][i] = (z * c - y * s);
y = U(j,nm);
z = U(j,i);
U(j,nm) = (y * c + z * s);
U(j,i) = (z * c - y * s);
}
}
}
}
z = Sigma[k][k];
z = Sigma(k,k);
if (l == k)
{// convergence
if (z < 0.0)
{// make singular value nonnegative
Sigma[k][k] = -z;
for (j = 0; j < n; j++) VT[j][k] = -VT[j][k];
Sigma(k,k) = -z;
for (j = 0; j < (int)n; j++) V(j,k) = -V(j,k);
}
break;
}
......@@ -255,9 +284,9 @@ namespace INMOST
return false;
}
// shift from bottom 2 x 2 minor
x = Sigma[l][l];
x = Sigma(l,l);
nm = k - 1;
y = Sigma[nm][nm];
y = Sigma(nm,nm);
g = rv1[nm];
h = rv1[k];
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
......@@ -269,7 +298,7 @@ namespace INMOST
{
i = j + 1;
g = rv1[i];
y = Sigma[i][i];
y = Sigma(i,i);
h = s * g;
g = c * g;
z = pythag(f, h);
......@@ -280,15 +309,15 @@ namespace INMOST
g = g * c - x * s;
h = y * s;