From 50d10e74001525c92c754270e58b6290b1ecda0a Mon Sep 17 00:00:00 2001 From: Kirill Terekhov Date: Tue, 1 Dec 2015 10:57:03 -0800 Subject: [PATCH] Added functionality and fixed issues. Fixes for GCC. Moved functionality to operate with dense matrices as INMOST service. Added possibility to exchange data of type DATA_VARIABLE in parallel. Added possibility to include external content in xml format. Moved overloaded operators for variables from namespace INMOST into global scope to prevent possible ambiguities in compilers. Added perturbation algorithm for BCGS solver. --- CMakeLists.txt | 2 + Examples/ADMFD/CMakeLists.txt | 2 +- Examples/ADMFD/main.cpp | 27 +- Examples/ADMFD/matrix.hpp | 453 ------------------ Source/Autodiff/autodiff.cpp | 9 +- Source/Data/tag.cpp | 23 +- Source/Headers/CMakeLists.txt | 1 + Source/Headers/inmost.h | 1 + Source/Headers/inmost_autodiff.h | 107 +++-- Source/Headers/inmost_common.h | 7 +- Source/Headers/inmost_data.h | 9 + Source/Headers/inmost_dense.h | 710 +++++++++++++++++++++++++++++ Source/Headers/inmost_expression.h | 180 +++----- Source/Headers/inmost_mesh.h | 13 +- Source/Headers/inmost_solver.h | 3 +- Source/Headers/inmost_sparse.h | 5 + Source/Headers/inmost_variable.h | 79 ++-- Source/IO/mesh_xml_file.cpp | 110 +++-- Source/Mesh/geometry.cpp | 34 +- Source/Mesh/mesh.cpp | 129 ++++++ Source/Mesh/parallel.cpp | 21 +- Source/Solver/solver.cpp | 55 +-- Source/Solver/solver_bcgsl.hpp | 78 +++- Source/Solver/sparse.cpp | 48 ++ 24 files changed, 1309 insertions(+), 797 deletions(-) delete mode 100644 Examples/ADMFD/matrix.hpp create mode 100644 Source/Headers/inmost_dense.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 056fb5f..ea6a5ab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -240,6 +240,7 @@ set(INMOST_INSTALL_HEADERS Source/Headers/inmost.h Source/Headers/inmost_autodiff.h Source/Headers/inmost_common.h Source/Headers/inmost_data.h + Source/Headers/inmost_dense.h Source/Headers/inmost_expression.h Source/Headers/inmost_mesh.h Source/Headers/inmost_nonlinear.h @@ -268,6 +269,7 @@ set_property(TARGET inmost PROPERTY PUBLIC_HEADER "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_autodiff.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_common.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_data.h" + "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_dense.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_expression.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_mesh.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_nonlinear.h" diff --git a/Examples/ADMFD/CMakeLists.txt b/Examples/ADMFD/CMakeLists.txt index 4e19bef..1397294 100644 --- a/Examples/ADMFD/CMakeLists.txt +++ b/Examples/ADMFD/CMakeLists.txt @@ -1,5 +1,5 @@ project(ADMFD) -set(SOURCE main.cpp matrix.hpp) +set(SOURCE main.cpp) add_executable(ADMFD ${SOURCE}) diff --git a/Examples/ADMFD/main.cpp b/Examples/ADMFD/main.cpp index 412087d..b7faf07 100644 --- a/Examples/ADMFD/main.cpp +++ b/Examples/ADMFD/main.cpp @@ -4,7 +4,6 @@ #include #include "inmost.h" -#include "matrix.hpp" using namespace INMOST; #ifndef M_PI @@ -184,8 +183,8 @@ int main(int argc,char ** argv) rMatrix nKL(3,1); //normal to face rMatrix lKs(3,1); //reminder of the co-normal vector Kn in cell K when projection to normal is substracted, i.e. (Kn - n n.Kn) rMatrix lLs(3,1); //reminder of the co-normal vector Kn in cell L when projection to normal is substracted, i.e. (Kn - n n.Kn) - rMatrix KK = FromTensor(cK->RealArray(tag_K)).Transpose(); //diffusion tensor in cell K - rMatrix KL = FromTensor(cL->RealArray(tag_K)).Transpose(); //diffusion tensor in cell L + rMatrix KK = rMatrix::FromTensor(cK->RealArray(tag_K).data(),cK->RealArray(tag_K).size()).Transpose(); //diffusion tensor in cell K + rMatrix KL = rMatrix::FromTensor(cL->RealArray(tag_K).data(),cL->RealArray(tag_K).size()).Transpose(); //diffusion tensor in cell L cK->Centroid(xK.data()); //retrive center of cell K cL->Centroid(xL.data()); //retrive center of cell L @@ -237,7 +236,7 @@ int main(int argc,char ** argv) real aF, vP = cell->Volume(); //area of the face cell->Centroid(xP); //obtain cell center ElementArray faces = cell->getFaces(); //obtain faces of the cell - enumerator NF = faces.size(), k; //number of faces; + enumerator NF = (enumerator)faces.size(), k; //number of faces; rMatrix IP(NF,NF), N(NF,3), R(NF,3); //matrices for inner product @@ -261,8 +260,8 @@ int main(int argc,char ** argv) } //end of loop over faces //inner product definition from Corollary 4.2, "Mimetic scalar products of differential forms" by Brezzi et al - IP = R*(R.Transpose()*N).Invert()*R.Transpose(); // Consistency part - IP += (rMatrix::Unit(NF) - N*(N.Transpose()*N).Invert()*N.Transpose())*(2.0/static_cast(NF)*IP.Trace()); //Stability part + IP = R*(R.Transpose()*N).Invert().first*R.Transpose(); // Consistency part + IP += (rMatrix::Unit(NF) - N*(N.Transpose()*N).Invert().first*N.Transpose())*(2.0/static_cast(NF)*IP.Trace()); //Stability part //assert(IP.isSymmetric()); //test positive definitiness as well! /* @@ -288,9 +287,9 @@ int main(int argc,char ** argv) real nF[3]; //normal to the face cell->Centroid(xP); ElementArray faces = cell->getFaces(); //obtain faces of the cell - enumerator NF = faces.size(), k; //number of faces; + enumerator NF = (enumerator)faces.size(), k; //number of faces; dynarray NG(NF,0), G(NF,0); // count how many times gradient was calculated for face - rMatrix K = FromTensor(cell->RealArrayDF(tag_K)); //get permeability for the cell + rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),cell->RealArrayDF(tag_K).size()); //get permeability for the cell rMatrix GRAD(NF,NF), NK(NF,3), R(NF,3); //big gradient matrix, co-normals, directions rMatrix GRADF(NF,NF), NKF(NF,3), RF(NF,3); //corresponding matrices with masked negative values when multiplied by individual co-normals rMatrix GRADFstab(NF,NF); @@ -305,18 +304,18 @@ int main(int argc,char ** argv) R(k,1) = (yF[1]-xP[1]); R(k,2) = (yF[2]-xP[2]); // assemble matrix of co-normals - rMatrix nK = FromVector(nF).Transpose()*K; + rMatrix nK = rMatrix::FromVector(nF,3).Transpose()*K; NK(k,0) = nK(0,0); NK(k,1) = nK(0,1); NK(k,2) = nK(0,2); k++; } //end of loop over faces - GRAD = NK*(NK.Transpose()*R).Invert()*NK.Transpose(); //stability part + GRAD = NK*(NK.Transpose()*R).Invert().first*NK.Transpose(); //stability part //std::cout << "consistency" << std::endl; //GRADF.Print(); - GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert()*R.Transpose())*(2.0/NF*GRADF.Trace()); + GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert().first*R.Transpose())*(2.0/NF*GRADF.Trace()); /* std::cout << "stability" << std::endl; @@ -400,7 +399,7 @@ int main(int argc,char ** argv) for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell) { ElementArray faces = cell->getFaces(); //obtain faces of the cell - enumerator NF = faces.size(), k = 0; + enumerator NF = (enumerator)faces.size(), k = 0; Cell cK = cell->self(); real gK = cK->Real(tag_F0), gL; //force on current cell and on adjacent cell rMatrix M(cK->RealArrayDF(tag_M).data(),NF,NF); //inner product matrix @@ -475,7 +474,7 @@ int main(int argc,char ** argv) for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell) //loop over cells { ElementArray faces = cell->getFaces(); //obtain faces of the cell - enumerator NF = faces.size(), k = 0; + enumerator NF = (enumerator)faces.size(), k = 0; Cell cK = cell->self(); variable pK = P(cK), pL; //pressure for back and front cells rMatrix GRAD(cK->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient @@ -509,7 +508,7 @@ int main(int argc,char ** argv) for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell) //loop over cells { ElementArray facesK = cell->getFaces(), facesL; //obtain faces of the cell - enumerator NF = facesK.size(), k, l; + enumerator NF = (enumerator)facesK.size(), k, l; Cell cK = cell->self(); variable pK = P(cK), pL; //pressure for back and front cells rMatrix GRAD(cK->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient diff --git a/Examples/ADMFD/matrix.hpp b/Examples/ADMFD/matrix.hpp deleted file mode 100644 index 1246daf..0000000 --- a/Examples/ADMFD/matrix.hpp +++ /dev/null @@ -1,453 +0,0 @@ -#ifndef __MATRIX_H -#define __MATRIX_H -#include "inmost.h" - -using namespace INMOST; -// Matrix with n columns and m rows -// __m__ -// | | -// n| | -// |_____| -// -// todo: -// 1. expression templates for operations -// 2. (ok) template matrix type for AD variables -template -class Matrix -{ -public: - typedef unsigned enumerator; -protected: - array space; - enumerator n, m; -public: - static Matrix Unit(enumerator pn) - { - Matrix ret(pn,pn); - for(enumerator i = 0; i < pn; ++i) - { - for(enumerator j = 0; j < pn; ++j) - { - if( i == j ) - ret(i,j) = 1; - else - ret(i,j) = 0; - } - } - return ret; - } - Matrix(Var * pspace, enumerator pn, enumerator pm) : space(pspace,pspace+pn*pm), n(pn), m(pm) {} - Matrix(enumerator pn, enumerator pm) : space(pn*pm), n(pn), m(pm) {} - Matrix(const Matrix & other) : space(other.n*other.m), n(other.n), m(other.m) - { - for(enumerator i = 0; i < n*m; ++i) - space[i] = other.space[i]; - } - ~Matrix() {remove();} - Matrix & operator =(Matrix const & other) - { - remove(); - space.resize(other.n*other.m); - for(enumerator i = 0; i < other.n*other.m; ++i) - space[i] = other.space[i]; - n = other.n; - m = other.m; - return *this; - } - void remove() {} - // i is in [0,n] - row index - // j is in [0,m] - column index - Var & operator()(enumerator i, enumerator j) - { - assert(i >= 0 && i < n); - assert(j >= 0 && j < m); - assert(i*m+j < n*m); //overflow check? - return space[i*m+j]; - } - Var operator()(enumerator i, enumerator j) const - { - assert(i >= 0 && i < n); - assert(j >= 0 && j < m); - assert(i*m+j < n*m); //overflow check? - return space[i*m+j]; - } - Matrix operator-(const Matrix & other) const - { - assert(n == other.n); - assert(m == other.m); - Matrix ret(n,m); //check RVO - for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]-other.space[k]; - return ret; - } - Matrix & operator-=(const Matrix & other) - { - assert(n == other.n); - assert(m == other.m); - for(enumerator k = 0; k < n*m; ++k) space[k] -= other.space[k]; - return *this; - } - Matrix operator+(const Matrix & other) const - { - assert(n == other.n); - assert(m == other.m); - Matrix ret(n,m); //check RVO - for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]+other.space[k]; - return ret; - } - Matrix & operator+=(const Matrix & other) - { - assert(n == other.n); - assert(m == other.m); - for(enumerator k = 0; k < n*m; ++k) space[k] += other.space[k]; - return *this; - } - Matrix operator*(Var coef) const - { - Matrix ret(n,m); //check RVO - for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]*coef; - return ret; - } - Matrix & operator*=(Var coef) - { - for(enumerator k = 0; k < n*m; ++k) space[k] *= coef; - return *this; - } - Matrix operator/(Var coef) const - { - Matrix ret(n,m); //check RVO - for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]/coef; - return ret; - } - Matrix & operator/=(Var coef) - { - for(enumerator k = 0; k < n*m; ++k) space[k] /= coef; - return *this; - } - Matrix operator*(const Matrix & other) const - { - assert(m == other.n); - Matrix ret(n,other.m); //check RVO - for(enumerator i = 0; i < n; ++i) //loop rows - { - for(enumerator j = 0; j < other.m; ++j) //loop columns - { - Var tmp = 0.0; - for(enumerator k = 0; k < m; ++k) - tmp += (*this)(i,k)*other(k,j); - ret(i,j) = tmp; - } - } - return ret; - } - /* - template - Matrix operator*(const Matrix & other) - { - assert(m == other.n); - Matrix ret(n,other.m); //check RVO - for(enumerator i = 0; i < n; ++i) //loop rows - { - for(enumerator j = 0; j < other.m; ++j) //loop columns - { - InType tmp = 0.0; - for(enumerator k = 0; k < m; ++k) - tmp += (*this)(i,k)*other(k,j); - ret(i,j) = tmp; - } - } - return ret; - } - template - Matrix operator*(const Matrix & other) - { - assert(m == other.n); - Matrix ret(n,other.m); //check RVO - for(enumerator i = 0; i < n; ++i) //loop rows - { - for(enumerator j = 0; j < other.m; ++j) //loop columns - { - Var tmp = 0.0; - for(enumerator k = 0; k < m; ++k) - tmp += (*this)(i,k)*other(k,j); - ret(i,j) = tmp; - } - } - return ret; - } - */ - Matrix Transpose() const - { - Matrix ret(m,n); - for(enumerator i = 0; i < n; ++i) - { - for(enumerator j = 0; j < m; ++j) - { - ret(j,i) = (*this)(i,j); - } - } - return ret; - } - Matrix Invert() const - { - Matrix ret(m,n); - Matrix At = Transpose(); //m by n matrix - Matrix AtB = At*Unit(n); //m by n matrix - Matrix AtA = At*(*this); //m by m matrix - enumerator * order = new enumerator [m]; - for(enumerator i = 0; i < m; ++i) order[i] = i; - for(enumerator i = 0; i < m; i++) - { - enumerator maxk = i, maxq = i, temp2; - Var max, temp; - max = fabs(AtA(maxk,maxq)); - //Find best pivot - for(enumerator k = i; k < m; k++) // over rows - { - for(enumerator q = i; q < m; q++) // over columns - { - if( fabs(AtA(k,q)) > max ) - { - max = fabs(AtA(k,q)); - maxk = k; - maxq = q; - } - } - } - //Exchange rows - if( maxk != i ) - { - for(enumerator q = 0; q < m; q++) // over columns of A - { - temp = AtA(maxk,q); - AtA(maxk,q) = AtA(i,q); - AtA(i,q) = temp; - } - //exchange rhs - for(enumerator q = 0; q < n; q++) // over columns of B - { - temp = AtB(maxk,q); - AtB(maxk,q) = AtB(i,q); - AtB(i,q) = temp; - } - } - //Exchange columns - if( maxq != i ) - { - for(enumerator k = 0; k < m; k++) //over rows - { - temp = AtA(k,maxq); - AtA(k,maxq) = AtA(k,i); - AtA(k,i) = temp; - } - //remember order in sol - { - temp2 = order[maxq]; - order[maxq] = order[i]; - order[i] = temp2; - } - } - - if( fabs(AtA(i,i)) < 1.0e-54 ) - { - bool ok = true; - for(enumerator k = 0; k < n; k++) // over columns of B - { - if( fabs(AtB(i,k)/1.0e-54) > 1 ) - { - ok = false; - break; - } - } - if( ok ) AtA(i,i) = AtA(i,i) < 0.0 ? - 1.0e-12 : 1.0e-12; - else assert(false && "Matrix is singular"); - } - for(enumerator k = i+1; k < m; k++) - { - AtA(i,k) /= AtA(i,i); - AtA(k,i) /= AtA(i,i); - } - for(enumerator k = i+1; k < m; k++) - for(enumerator q = i+1; q < m; q++) - { - AtA(k,q) -= AtA(k,i) * AtA(i,i) * AtA(i,q); - } - for(enumerator k = 0; k < n; k++) - { - for(enumerator j = i+1; j < m; j++) //iterate over columns of L - { - AtB(j,k) -= AtB(i,k) * AtA(j,i); - } - AtB(i,k) /= AtA(i,i); - } - } - - for(enumerator k = 0; k < n; k++) - { - for(enumerator i = m; i-- > 0; ) //iterate over rows of U - for(enumerator j = i+1; j < m; j++) - { - AtB(i,k) -= AtB(j,k) * AtA(i,j); - } - for(enumerator i = 0; i < m; i++) - ret(order[i],k) = AtB(i,k); - } - delete [] order; - return ret; - } - void Zero() - { - for(enumerator i = 0; i < n*m; ++i) space[i] = 0.0; - } - Var Trace() const - { - assert(n == m); - Var ret = 0.0; - for(enumerator i = 0; i < n; ++i) ret += (*this)(i,i); - return ret; - } - Var * data() {return space.data();} - const Var * data() const {return space.data();} - enumerator Rows() const {return n;} - enumerator Cols() const {return m;} - void Print(Storage::real threshold = 1.0e-10) const - { - for(enumerator k = 0; k < n; ++k) - { - for(enumerator l = 0; l < m; ++l) - { - if( fabs((*this)(k,l)) > threshold ) - std::cout << (*this)(k,l); - else - std::cout << 0; - std::cout << " "; - } - std::cout << std::endl; - } - } - bool isSymmetric() const - { - if( n != m ) return false; - for(enumerator k = 0; k < n; ++k) - { - for(enumerator l = k+1; l < n; ++l) - if( fabs((*this)(k,l)-(*this)(l,k)) > 1.0e-7 ) - return false; - } - return true; - } - Var DotProduct(const Matrix & other) const - { - assert(n == other.n); - assert(m == other.m); - Var ret = 0.0; - for(enumerator i = 0; i < n*m; ++i) ret += space[i]*other.space[i]; - return ret; - } -}; - - -typedef Matrix rMatrix; //shortcut for real matrix -typedef Matrix vMatrix; //shortcut for matrix with variations - - - - -// Multiplication of matrices with promotion -__INLINE vMatrix operator *(const rMatrix & A, const vMatrix & B) -{ - assert(A.Cols() == B.Rows()); - vMatrix ret(A.Rows(),B.Cols()); //check RVO - for(rMatrix::enumerator i = 0; i < A.Rows(); ++i) //loop rows - { - for(rMatrix::enumerator j = 0; j < B.Cols(); ++j) //loop columns - { - variable tmp = 0.0; - for(rMatrix::enumerator k = 0; k < A.Cols(); ++k) - tmp += A(i,k)*B(k,j); - ret(i,j) = tmp; - } - } - return ret; -} - -// Multiplication of matrices with promotion -__INLINE vMatrix operator *(const vMatrix & A, const rMatrix & B) -{ - assert(A.Cols() == B.Rows()); - vMatrix ret(A.Rows(),B.Cols()); //check RVO - for(rMatrix::enumerator i = 0; i < A.Rows(); ++i) //loop rows - { - for(rMatrix::enumerator j = 0; j < B.Cols(); ++j) //loop columns - { - variable tmp = 0.0; - for(rMatrix::enumerator k = 0; k < A.Cols(); ++k) - tmp += A(i,k)*B(k,j); - ret(i,j) = tmp; - } - } - return ret; -} - -//retrive permeability matrix from stored data -__INLINE rMatrix FromTensor(Storage::real_array & K) -{ - rMatrix Kc(3,3); - switch(K.size()) - { - case 1: //scalar permeability tensor - Kc.Zero(); - Kc(0,0) = Kc(1,1) = Kc(2,2) = K[0]; - case 3: - Kc.Zero(); //diagonal permeability tensor - Kc(0,0) = K[0]; //KXX - Kc(1,1) = K[1]; //KYY - Kc(2,2) = K[2]; //KZZ - break; - case 6: //symmetric permeability tensor - Kc(0,0) = K[0]; //KXX - Kc(0,1) = Kc(1,0) = K[1]; //KXY - Kc(0,2) = Kc(2,0) = K[2]; //KXZ - Kc(1,1) = K[3]; //KYY - Kc(1,2) = Kc(2,1) = K[4]; //KYZ - Kc(2,2) = K[5]; //KZZ - break; - case 9: //full permeability tensor - Kc(0,0) = K[0]; //KXX - Kc(0,1) = K[1]; //KXY - Kc(0,2) = K[2]; //KXZ - Kc(1,0) = K[3]; //KYX - Kc(1,1) = K[4]; //KYY - Kc(1,2) = K[5]; //KYZ - Kc(2,0) = K[6]; //KZX - Kc(2,1) = K[7]; //KZY - Kc(2,2) = K[8]; //KZZ - break; - } - return Kc; -} - -//retrive vector in matrix form from array -__INLINE rMatrix FromVector(Storage::real n[3]) -{ - return rMatrix(n,3,1); -} - -//create diagonal matrix from array -__INLINE rMatrix FromDiagonal(Storage::real * r, rMatrix::enumerator size) -{ - rMatrix ret(size,size); - ret.Zero(); - for(rMatrix::enumerator k = 0; k < size; ++k) ret(k,k) = r[k]; - return ret; -} - -//create diagonal matrix from array of inversed values -__INLINE rMatrix FromDiagonalInverse(Storage::real * r, rMatrix::enumerator size) -{ - rMatrix ret(size,size); - ret.Zero(); - for(rMatrix::enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k]; - return ret; -} - -#endif \ No newline at end of file diff --git a/Source/Autodiff/autodiff.cpp b/Source/Autodiff/autodiff.cpp index cadabb7..455452e 100644 --- a/Source/Autodiff/autodiff.cpp +++ b/Source/Autodiff/autodiff.cpp @@ -201,15 +201,16 @@ namespace INMOST INMOST_DATA_REAL_TYPE Automatizator::Derivative(expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply, void * user_data) { + Sparse::RowMerger & m = GetMerger(); INMOST_DATA_REAL_TYPE ret; var.current_stencil.resize(1); var.current_stencil[0] = stencil_pair(e->GetHandle(),1.0); var.resize_for_stencil(); ret = EvaluateSub(var,0,ENUMUNDEF,user_data); - merger.PushRow(1.0,out); - DerivativeFill(var, 0, ENUMUNDEF, merger, multiply, user_data); - merger.RetriveRow(out); - merger.Clear(); + m.PushRow(1.0,out); + DerivativeFill(var, 0, ENUMUNDEF, m, multiply, user_data); + m.RetriveRow(out); + m.Clear(); return ret*multiply; } #else diff --git a/Source/Data/tag.cpp b/Source/Data/tag.cpp index 54cf68e..7e377dc 100644 --- a/Source/Data/tag.cpp +++ b/Source/Data/tag.cpp @@ -11,7 +11,9 @@ namespace INMOST case DATA_INTEGER: return "INTEGER"; case DATA_BULK: return "BULK"; case DATA_REFERENCE: return "REFERENCE"; +#if defined(USE_AUTODIFF) case DATA_VARIABLE: return "VARIABLE"; +#endif } return "UNKNOWN"; } @@ -25,7 +27,9 @@ namespace INMOST case DATA_INTEGER: return sizeof(INMOST_DATA_INTEGER_TYPE); case DATA_REAL: return sizeof(INMOST_DATA_REAL_TYPE); case DATA_REFERENCE: return sizeof(HandleType); +#if defined(USE_AUTODIFF) case DATA_VARIABLE: return sizeof(variable); +#endif } return 0; } @@ -39,7 +43,9 @@ namespace INMOST case DATA_INTEGER: return sizeof(inner_integer_array); case DATA_BULK: return sizeof(inner_bulk_array); case DATA_REFERENCE: return sizeof(inner_reference_array); +#if defined(USE_AUTODIFF) case DATA_VARIABLE: return sizeof(inner_variable_array); +#endif } return 0; } @@ -73,13 +79,17 @@ namespace INMOST else if( type == DATA_INTEGER ) new (adata) inner_integer_array (*static_cast(bdata)); else if( type == DATA_BULK ) new (adata) inner_bulk_array (*static_cast(bdata)); else if( type == DATA_REFERENCE ) new (adata) inner_reference_array(*static_cast(bdata)); +#if defined(USE_AUTODIFF) else if( type == DATA_VARIABLE ) new (adata) inner_variable_array (*static_cast(bdata)); +#endif } +#if defined(USE_AUTODIFF) else if( type == DATA_VARIABLE ) //have to call constructor { for(INMOST_DATA_ENUM_TYPE k = 0; k < data_size; ++k) new (static_cast(adata)+k) variable(*(static_cast(bdata)+k)); } +#endif else // fixed size array memcpy(adata,bdata,data_size*bytes); } @@ -109,11 +119,13 @@ namespace INMOST (*static_cast (adata)).~inner_reference_array(); new (adata) inner_reference_array(); } +#if defined(USE_AUTODIFF) else if( type == DATA_VARIABLE ) { (*static_cast (adata)).~inner_variable_array(); new (adata) inner_variable_array(); } +#endif } } @@ -177,7 +189,12 @@ namespace INMOST case DATA_REAL: mem->bulk_data_type = INMOST_MPI_DATA_REAL_TYPE; break; case DATA_INTEGER: mem->bulk_data_type = INMOST_MPI_DATA_INTEGER_TYPE; break; case DATA_REFERENCE: mem->bulk_data_type = INMOST_MPI_DATA_ENUM_TYPE; break; - case DATA_VARIABLE: mem->bulk_data_type = INMOST_MPI_DATA_REAL_TYPE; break; //mixture of types, should not get here +#if defined(USE_AUTODIFF) + case DATA_VARIABLE: + if( !Sparse::HaveRowEntryType() ) Sparse::CreateRowEntryType(); + mem->bulk_data_type = Sparse::GetRowEntryType(); + break; +#endif } mem->bytes_size = DataTypeBytesSize(mem->dtype); if(mem->size == ENUMUNDEF ) @@ -481,7 +498,9 @@ namespace INMOST case DATA_INTEGER: for(INMOST_DATA_ENUM_TYPE it = new_size; it < old_size; ++it) {void * p = static_cast(&arr[it]); if( p != NULL ) (*static_cast( p )).~inner_integer_array(); } break; case DATA_BULK: for(INMOST_DATA_ENUM_TYPE it = new_size; it < old_size; ++it) {void * p = static_cast(&arr[it]); if( p != NULL ) (*static_cast( p )).~inner_bulk_array(); } break; case DATA_REFERENCE: for(INMOST_DATA_ENUM_TYPE it = new_size; it < old_size; ++it) {void * p = static_cast(&arr[it]); if( p != NULL ) (*static_cast( p )).~inner_reference_array();} break; +#if defined(USE_AUTODIFF) case DATA_VARIABLE: for(INMOST_DATA_ENUM_TYPE it = new_size; it < old_size; ++it) {void * p = static_cast(&arr[it]); if( p != NULL ) (*static_cast( p )).~inner_variable_array(); } break; +#endif } } } @@ -499,7 +518,9 @@ namespace INMOST case DATA_INTEGER: for(INMOST_DATA_ENUM_TYPE it = old_size; it < new_size; ++it) new ( &arr[it] ) inner_integer_array(); break; case DATA_BULK: for(INMOST_DATA_ENUM_TYPE it = old_size; it < new_size; ++it) new ( &arr[it] ) inner_bulk_array(); break; case DATA_REFERENCE: for(INMOST_DATA_ENUM_TYPE it = old_size; it < new_size; ++it) new ( &arr[it] ) inner_reference_array(); break; +#if defined(USE_AUTODIFF) case DATA_VARIABLE: for(INMOST_DATA_ENUM_TYPE it = old_size; it < new_size; ++it) new ( &arr[it] ) inner_variable_array(); break; +#endif } } } diff --git a/Source/Headers/CMakeLists.txt b/Source/Headers/CMakeLists.txt index 8f5ed1d..de881f7 100644 --- a/Source/Headers/CMakeLists.txt +++ b/Source/Headers/CMakeLists.txt @@ -9,6 +9,7 @@ set(HEADER ${CMAKE_CURRENT_SOURCE_DIR}/inmost_options_cmake.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_common.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_data.h + ${CMAKE_CURRENT_SOURCE_DIR}/inmost_dense.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_mesh.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_partitioner.h diff --git a/Source/Headers/inmost.h b/Source/Headers/inmost.h index e3fd39f..acf7b39 100644 --- a/Source/Headers/inmost.h +++ b/Source/Headers/inmost.h @@ -4,6 +4,7 @@ #include "inmost_common.h" #include "inmost_mesh.h" +#include "inmost_dense.h" #include "inmost_solver.h" #include "inmost_partitioner.h" #include "inmost_variable.h" diff --git a/Source/Headers/inmost_autodiff.h b/Source/Headers/inmost_autodiff.h index 75458d9..e19f123 100644 --- a/Source/Headers/inmost_autodiff.h +++ b/Source/Headers/inmost_autodiff.h @@ -602,32 +602,42 @@ namespace INMOST friend class Automatizator; }; - +} - __INLINE expr operator +(INMOST_DATA_REAL_TYPE left, const expr & right) { return expr(left) + right; } - __INLINE expr operator -(INMOST_DATA_REAL_TYPE left, const expr & right) { return expr(left) - right; } - __INLINE expr operator *(INMOST_DATA_REAL_TYPE left, const expr & right) { return expr(left) * right; } - __INLINE expr operator /(INMOST_DATA_REAL_TYPE left, const expr & right) { return expr(left) / right; } - __INLINE expr ad_pow(const expr & v, const expr & n) { return expr(AD_POW, v, n); } - //__INLINE expr ad_pow(const expr & v, INMOST_DATA_REAL_TYPE n) { return expr(AD_POW, v, expr(n)); } - //__INLINE expr ad_pow(const expr & v, INMOST_DATA_ENUM_TYPE n) { return expr(AD_POW, v, expr(static_cast(n))); } - __INLINE expr ad_abs(const expr & v) { return expr(AD_ABS, v); } - __INLINE expr ad_exp(const expr & v) { return expr(AD_EXP, v); } - __INLINE expr ad_log(const expr & v) { return expr(AD_LOG, v); } - __INLINE expr ad_sin(const expr & v) { return expr(AD_SIN, v); } - __INLINE expr ad_cos(const expr & v) { return expr(AD_COS, v); } - __INLINE expr ad_sqrt(const expr & v) {return expr(AD_SQRT, v);} - __INLINE expr ad_val(const expr & v, const expr & multiplyer = expr(0.0)) {return expr(v, multiplyer);} - __INLINE expr measure() { return expr(AD_MES,ENUMUNDEF); } - __INLINE expr condition(const expr & cond, const expr & if_gt_zero, const expr & if_le_zero) { return expr(cond, if_gt_zero, if_le_zero); } - __INLINE expr condition_etype(ElementType etypes, const expr & if_true, const expr & if_false) { return expr(expr(AD_COND_TYPE, etypes), if_true, if_false); } - __INLINE expr condition_marker(MarkerType marker, const expr & if_true, const expr & if_false) { return expr(expr(AD_COND_MARK, marker), if_true, if_false); } - __INLINE expr stencil(INMOST_DATA_ENUM_TYPE stncl, const expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return expr(stncl, &v); } - __INLINE expr tabval(INMOST_DATA_ENUM_TYPE tabl, const expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return expr(tabl, v); } - __INLINE expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return expr(reg_tag, comp); } - __INLINE expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return expr(reg_func, ENUMUNDEF); } +__INLINE INMOST::expr operator +(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) + right; } +__INLINE INMOST::expr operator -(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) - right; } +__INLINE INMOST::expr operator *(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) * right; } +__INLINE INMOST::expr operator /(INMOST_DATA_REAL_TYPE left, const INMOST::expr & right) { return INMOST::expr(left) / right; } +__INLINE INMOST::expr ad_pow(const INMOST::expr & v, const INMOST::expr & n) { return INMOST::expr(AD_POW, v, n); } +//__INLINE INMOST::expr ad_pow(const INMOST::expr & v, INMOST_DATA_REAL_TYPE n) { return INMOST::expr(AD_POW, v, INMOST::expr(n)); } +//__INLINE INMOST::expr ad_pow(const INMOST::expr & v, INMOST_DATA_ENUM_TYPE n) { return INMOST::expr(AD_POW, v, INMOST::expr(static_cast(n))); } +__INLINE INMOST::expr ad_abs(const INMOST::expr & v) { return INMOST::expr(AD_ABS, v); } +__INLINE INMOST::expr ad_exp(const INMOST::expr & v) { return INMOST::expr(AD_EXP, v); } +__INLINE INMOST::expr ad_log(const INMOST::expr & v) { return INMOST::expr(AD_LOG, v); } +__INLINE INMOST::expr ad_sin(const INMOST::expr & v) { return INMOST::expr(AD_SIN, v); } +__INLINE INMOST::expr ad_cos(const INMOST::expr & v) { return INMOST::expr(AD_COS, v); } +__INLINE INMOST::expr ad_sqrt(const INMOST::expr & v) {return INMOST::expr(AD_SQRT, v);} +__INLINE INMOST::expr ad_val(const INMOST::expr & v, const INMOST::expr & multiplyer = INMOST::expr(0.0)) {return INMOST::expr(v, multiplyer);} +__INLINE INMOST::expr measure() { return INMOST::expr(AD_MES,ENUMUNDEF); } +__INLINE INMOST::expr condition(const INMOST::expr & cond, const INMOST::expr & if_gt_zero, const INMOST::expr & if_le_zero) { return INMOST::expr(cond, if_gt_zero, if_le_zero); } +__INLINE INMOST::expr condition_etype(INMOST::ElementType etypes, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(INMOST::expr(AD_COND_TYPE, etypes), if_true, if_false); } +__INLINE INMOST::expr condition_marker(INMOST::MarkerType marker, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(INMOST::expr(AD_COND_MARK, marker), if_true, if_false); } +__INLINE INMOST::expr stencil(INMOST_DATA_ENUM_TYPE stncl, const INMOST::expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return INMOST::expr(stncl, &v); } +__INLINE INMOST::expr tabval(INMOST_DATA_ENUM_TYPE tabl, const INMOST::expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return INMOST::expr(tabl, v); } +__INLINE INMOST::expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return INMOST::expr(reg_tag, comp); } +__INLINE INMOST::expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return INMOST::expr(reg_func, ENUMUNDEF); } +namespace INMOST +{ #else +} +__INLINE INMOST::expr operator +(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right); +__INLINE INMOST::expr operator -(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right); +__INLINE INMOST::expr operator *(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right); +__INLINE INMOST::expr operator /(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right); + +namespace INMOST +{ class expr { INMOST_DATA_ENUM_TYPE op; @@ -683,36 +693,39 @@ namespace INMOST expr operator -(const INMOST_DATA_REAL_TYPE & other) const { return expr(AD_MINUS, new expr(*this), new expr(other)); } expr operator *(const INMOST_DATA_REAL_TYPE & other) const { expr var(*this); var.coef *= other; return var; } expr operator /(const INMOST_DATA_REAL_TYPE & other) const { expr var(*this); var.coef /= other; return var; } - __INLINE friend expr operator +(const INMOST_DATA_REAL_TYPE & left, const expr & right); - __INLINE friend expr operator -(const INMOST_DATA_REAL_TYPE & left, const expr & right); - __INLINE friend expr operator *(const INMOST_DATA_REAL_TYPE & left, const expr & right); - __INLINE friend expr operator /(const INMOST_DATA_REAL_TYPE & left, const expr & right); + __INLINE friend expr (::operator +)(const INMOST_DATA_REAL_TYPE & left, const expr & right); + __INLINE friend expr (::operator -)(const INMOST_DATA_REAL_TYPE & left, const expr & right); + __INLINE friend expr (::operator *)(const INMOST_DATA_REAL_TYPE & left, const expr & right); + __INLINE friend expr (::operator /)(const INMOST_DATA_REAL_TYPE & left, const expr & right); bool operator ==(const expr & other) {return this == &other || (op == other.op && left == other.left && right == other.right);} bool is_endpoint() { return op >= AD_TAG && op < AD_STNCL; } friend class Automatizator; }; +} - __INLINE expr operator +(const INMOST_DATA_REAL_TYPE & left, const expr & right) { return expr(AD_PLUS, new expr(left), new expr(right)); } - __INLINE expr operator -(const INMOST_DATA_REAL_TYPE & left, const expr & right) { return expr(AD_MINUS, new expr(left), new expr(right)); } - __INLINE expr operator *(const INMOST_DATA_REAL_TYPE & left, const expr & right) { expr var(right); var.coef *= left; return var; } - __INLINE expr operator /(const INMOST_DATA_REAL_TYPE & left, const expr & right) { expr var; var.op = AD_INV; var.coef = left; var.left = new expr(right); var.right = NULL; return var; } - __INLINE expr ad_pow(const expr & v, const expr n) { return expr(AD_POW, new expr(v), new expr(n)); } - __INLINE expr ad_abs(const expr & v) { return expr(AD_ABS, new expr(v), NULL); } - __INLINE expr ad_exp(const expr & v) { return expr(AD_EXP, new expr(v), NULL); } - __INLINE expr ad_log(const expr & v) { return expr(AD_LOG, new expr(v), NULL); } - __INLINE expr ad_sin(const expr & v) { return expr(AD_SIN, new expr(v), NULL); } - __INLINE expr ad_cos(const expr & v) { return expr(AD_COS, new expr(v), NULL); } - __INLINE expr ad_sqrt(const expr & v) {return expr(AD_SQRT, new expr(v), NULL);} - __INLINE expr ad_val(const expr & v, const expr & multiplyer = expr(0.0)) {return expr(AD_VAL,new expr(v), new expr(multiplyer));} - __INLINE expr measure() { return expr(AD_MES, NULL, NULL); } - __INLINE expr condition_etype(ElementType etype, const expr & if_true, const expr & if_false) { return expr(AD_COND, new expr(AD_COND_TYPE,etype), new expr(AD_ALTR, new expr(if_true), new expr(if_false))); } - __INLINE expr condition_marker(MarkerType marker, const expr & if_true, const expr & if_false) { return expr(AD_COND, new expr(AD_COND_MARK,marker), new expr(AD_ALTR, new expr(if_true), new expr(if_false))); } - __INLINE expr condition(const expr & cond, const expr & if_gt_zero, const expr & if_le_zero) { return expr(AD_COND, new expr(cond), new expr(AD_ALTR, new expr(if_gt_zero), new expr(if_le_zero))); } - __INLINE expr stencil(INMOST_DATA_ENUM_TYPE stncl, const expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return expr(stncl, new expr(v), NULL); } - __INLINE expr tabval(INMOST_DATA_ENUM_TYPE tabl, const expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return expr(tabl, new expr(v), NULL); } - __INLINE expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return expr(reg_tag, comp); } - __INLINE expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return expr(reg_func, NULL, NULL); } +__INLINE INMOST::expr operator +(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { return INMOST::expr(AD_PLUS, new INMOST::expr(left), new INMOST::expr(right)); } +__INLINE INMOST::expr operator -(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { return INMOST::expr(AD_MINUS, new INMOST::expr(left), new INMOST::expr(right)); } +__INLINE INMOST::expr operator *(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { INMOST::expr var(right); var.coef *= left; return var; } +__INLINE INMOST::expr operator /(const INMOST_DATA_REAL_TYPE & left, const INMOST::expr & right) { INMOST::expr var; var.op = AD_INV; var.coef = left; var.left = new INMOST::expr(right); var.right = NULL; return var; } +__INLINE INMOST::expr ad_pow(const INMOST::expr & v, const INMOST::expr n) { return INMOST::expr(AD_POW, new INMOST::expr(v), new INMOST::expr(n)); } +__INLINE INMOST::expr ad_abs(const INMOST::expr & v) { return INMOST::expr(AD_ABS, new INMOST::expr(v), NULL); } +__INLINE INMOST::expr ad_exp(const INMOST::expr & v) { return INMOST::expr(AD_EXP, new INMOST::expr(v), NULL); } +__INLINE INMOST::expr ad_log(const INMOST::expr & v) { return INMOST::expr(AD_LOG, new INMOST::expr(v), NULL); } +__INLINE INMOST::expr ad_sin(const INMOST::expr & v) { return INMOST::expr(AD_SIN, new INMOST::expr(v), NULL); } +__INLINE INMOST::expr ad_cos(const INMOST::expr & v) { return INMOST::expr(AD_COS, new INMOST::expr(v), NULL); } +__INLINE INMOST::expr ad_sqrt(const INMOST::expr & v) {return INMOST::expr(AD_SQRT, new INMOST::expr(v), NULL);} +__INLINE INMOST::expr ad_val(const INMOST::expr & v, const INMOST::expr & multiplyer = INMOST::expr(0.0)) {return INMOST::expr(AD_VAL,new INMOST::expr(v), new INMOST::expr(multiplyer));} +__INLINE INMOST::expr measure() { return INMOST::expr(AD_MES, NULL, NULL); } +__INLINE INMOST::expr condition_etype(INMOST::ElementType etype, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(AD_COND, new INMOST::expr(AD_COND_TYPE,etype), new INMOST::expr(AD_ALTR, new INMOST::expr(if_true), new INMOST::expr(if_false))); } +__INLINE INMOST::expr condition_marker(INMOST::MarkerType marker, const INMOST::expr & if_true, const INMOST::expr & if_false) { return INMOST::expr(AD_COND, new INMOST::expr(AD_COND_MARK,marker), new INMOST::expr(AD_ALTR, new INMOST::expr(if_true), new INMOST::expr(if_false))); } +__INLINE INMOST::expr condition(const INMOST::expr & cond, const INMOST::expr & if_gt_zero, const INMOST::expr & if_le_zero) { return INMOST::expr(AD_COND, new INMOST::expr(cond), new INMOST::expr(AD_ALTR, new INMOST::expr(if_gt_zero), new INMOST::expr(if_le_zero))); } +__INLINE INMOST::expr stencil(INMOST_DATA_ENUM_TYPE stncl, const INMOST::expr & v) { assert(stncl >= AD_STNCL && stncl < AD_TABLE); return INMOST::expr(stncl, new INMOST::expr(v), NULL); } +__INLINE INMOST::expr tabval(INMOST_DATA_ENUM_TYPE tabl, const INMOST::expr & v) { assert(tabl >= AD_TABLE && tabl < AD_FUNC); return INMOST::expr(tabl, new INMOST::expr(v), NULL); } +__INLINE INMOST::expr tagval(INMOST_DATA_ENUM_TYPE reg_tag, INMOST_DATA_ENUM_TYPE comp = 0) { assert(reg_tag >= AD_TAG && reg_tag < AD_STNCL); return INMOST::expr(reg_tag, comp); } +__INLINE INMOST::expr funcval(INMOST_DATA_ENUM_TYPE reg_func) { assert(reg_func >= AD_FUNC); return INMOST::expr(reg_func, NULL, NULL); } +namespace INMOST +{ #endif } diff --git a/Source/Headers/inmost_common.h b/Source/Headers/inmost_common.h index 83bf44b..d641ce4 100644 --- a/Source/Headers/inmost_common.h +++ b/Source/Headers/inmost_common.h @@ -99,7 +99,9 @@ #if defined(USE_OMP) #include #endif - +#if defined(USE_MPI) +#include +#endif #if !defined(USE_MPI) #define INMOST_MPI_Request int @@ -111,6 +113,7 @@ #define INMOST_MPI_DOUBLE 0 #define INMOST_MPI_UNSIGNED 0 #define INMOST_MPI_Win int +#define INMOST_MPI_DATATYPE_NULL 0 #else #define INMOST_MPI_Request MPI_Request #define INMOST_MPI_Type MPI_Datatype @@ -121,8 +124,10 @@ #define INMOST_MPI_DOUBLE MPI_DOUBLE #define INMOST_MPI_UNSIGNED MPI_UNSIGNED #define INMOST_MPI_Win MPI_Win +#define INMOST_MPI_DATATYPE_NULL MPI_DATATYPE_NULL #endif + #define INMOST_MPI_SIZE int //in case MPI standard changes and compiler gives tons of warning #define INMOST_DATA_INTEGER_TYPE int diff --git a/Source/Headers/inmost_data.h b/Source/Headers/inmost_data.h index f369ef1..ff6ebe4 100644 --- a/Source/Headers/inmost_data.h +++ b/Source/Headers/inmost_data.h @@ -435,12 +435,21 @@ namespace INMOST /// Return the data length associated with Tag. /// For abstract data return the number of bytes, otherwise return the length of associated array. + /// @param tag tag that represents the data /// @see Storage::SetDataSize + /// @see Mesh::GetDataSize __INLINE INMOST_DATA_ENUM_TYPE GetDataSize (const Tag & tag) const; + /// Return the size of the structure required to represent the data on current element. + /// This is equal to GetDataSize times Tag::GetBytesSize for all the data types, + /// except for DATA_VARIABLE, that requires a larger structure to accomodate derivatives. + /// @param tag tag that represents the data + /// @see Mesh::GetDataCapacity + __INLINE INMOST_DATA_ENUM_TYPE GetDataCapacity (const Tag & tag) const; /// Set the length of data associated with Tag. /// @param tag Identifying Tag. /// @param new_size The number of bytes for abstract data, otherwise the length of the array. /// @see Storage::GetDataSize + /// @see Mesh::SetDataSize __INLINE void SetDataSize (const Tag & tag, INMOST_DATA_ENUM_TYPE new_size) const; /// Extract part of the data associated with Tag. diff --git a/Source/Headers/inmost_dense.h b/Source/Headers/inmost_dense.h new file mode 100644 index 0000000..166aee3 --- /dev/null +++ b/Source/Headers/inmost_dense.h @@ -0,0 +1,710 @@ +#ifndef INMOST_DENSE_INCLUDED +#define INMOST_DENSE_INCLUDED +#include "inmost_common.h" +#if defined(USE_AUTODIFF) +#include "inmost_expression.h" +#endif + +// Matrix with n columns and m rows +// __m__ +// | | +// n| | +// |_____| +// +// todo: +// 1. expression templates for operations +// (???) how to for multiplication? +// 2. (ok) template matrix type for AD variables + + +namespace INMOST +{ + + + template struct Promote; + template<> struct Promote {typedef INMOST_DATA_REAL_TYPE type;}; +#if defined(USE_AUTODIFF) + template<> struct Promote {typedef variable type;}; + template<> struct Promote {typedef variable type;}; + template<> struct Promote {typedef variable type;}; +#else + INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE x) {return x;} +#endif + + template + class Matrix + { + public: + typedef unsigned enumerator; + protected: + array space; + 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 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); } + else result = 0.0; + return result; + } + public: + 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) + { + 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; + if (n < m) + { + bool success = Transpose().SVD(U,Sigma,VT); + if( success ) + { + U.Swap(VT); + U = U.Transpose(); + VT = VT.Transpose(); + return true; + } + else return false; + } //m <= n + array _rv1(m); + shell rv1(_rv1); + U = (*this); + Sigma.Resize(m,m); + Sigma.Zero(); + VT.Resize(m,m); + + std::swap(n,m); //this how original algorithm takes it + // Householder reduction to bidiagonal form + for (i = 0; i < n; i++) + { + // left-hand reduction + l = i + 1; + rv1[i] = scale * g; + g = s = scale = 0.0; + if (i < m) + { + for (k = i; k < m; k++) scale += fabs(U[k][i]); + if (get_value(scale)) + { + for (k = i; k < m; k++) + { + U[k][i] /= scale; + s += U[k][i] * U[k][i]; + } + f = U[i][i]; + g = -sign_func(sqrt(s), f); + h = f * g - s; + U[i][i] = f - g; + if (i != n - 1) + { + for (j = l; j < n; j++) + { + for (s = 0.0, k = i; k < 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 < m; k++) U[k][i] *= scale; + } + } + Sigma[i][i] = scale * g; + // right-hand reduction + g = s = scale = 0.0; + if (i < m && i != n - 1) + { + for (k = l; k < n; k++) scale += fabs(U[i][k]); + if (get_value(scale)) + { + for (k = l; k < n; k++) + { + U[i][k] = U[i][k]/scale; + s += U[i][k] * U[i][k]; + } + 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; + if (i != m - 1) + { + for (j = l; j < 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 (k = l; k < n; k++) U[i][k] *= scale; + } + } + anorm = max_func(anorm,fabs(Sigma[i][i]) + fabs(rv1[i])); + } + + // accumulate the right-hand transformation + for (i = n - 1; i >= 0; i--) + { + if (i < n - 1) + { + if (get_value(g)) + { + for (j = l; j < n; j++) VT[j][i] = ((U[i][j] / U[i][l]) / g); + // double division to avoid underflow + for (j = l; j < 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 (j = l; j < n; j++) VT[i][j] = VT[j][i] = 0.0; + } + VT[i][i] = 1.0; + g = rv1[i]; + l = i; + } + + // accumulate the left-hand transformation + 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; + if (get_value(g)) + { + g = 1.0 / g; + if (i != n - 1) + { + for (j = l; j < 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 (j = i; j < 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; + } + + // diagonalize the bidiagonal form + for (k = n - 1; k >= 0; k--) + {// loop over singular values + for (its = 0; its < 30; its++) + {// loop over allowed iterations + flag = 1; + for (l = k; l >= 0; l--) + {// test for splitting + nm = l - 1; + if (fabs(rv1[l]) + anorm == anorm) + { + flag = 0; + break; + } + if (fabs(Sigma[nm][nm]) + anorm == anorm) + break; + } + if (flag) + { + c = 0.0; + s = 1.0; + for (i = l; i <= k; i++) + { + f = s * rv1[i]; + if (fabs(f) + anorm != anorm) + { + g = Sigma[i][i]; + h = pythag(f, g); + Sigma[i][i] = h; + h = 1.0 / h; + c = g * h; + s = (- f * h); + for (j = 0; j < m; j++) + { + 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]; + 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]; + } + break; + } + if (its >= 30) + { + std::cout << "No convergence after " << its << " iterations" << std::endl; + std::swap(n,m); + return false; + } + // shift from bottom 2 x 2 minor + x = Sigma[l][l]; + nm = k - 1; + y = Sigma[nm][nm]; + g = rv1[nm]; + h = rv1[k]; + f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y); + g = pythag(f, 1.0); + f = ((x - z) * (x + z) + h * ((y / (f + sign_func(g, f))) - h)) / x; + // next QR transformation + c = s = 1.0; + for (j = l; j <= nm; j++) + { + i = j + 1; + g = rv1[i]; + y = Sigma[i][i]; + h = s * g; + g = c * g; + z = pythag(f, h); + rv1[j] = z; + c = f / z; + s = h / z; + f = x * c + g * s; + g = g * c - x * s; + h = y * s; + y = y * c; + for (jj = 0; jj < n; jj++) + { + x = VT[jj][j]; + z = VT[jj][i]; + VT[jj][j] = (x * c + z * s); + VT[jj][i] = (z * c - x * s); + } + z = pythag(f, h); + Sigma[j][j] = z; + if (z) + { + z = 1.0 / z; + c = f * z; + s = h * z; + } + f = (c * g) + (s * y); + x = (c * y) - (s * g); + for (jj = 0; jj < m; jj++) + { + y = U[jj][j]; + z = U[jj][i]; + U[jj][j] = (y * c + z * s); + U[jj][i] = (z * c - y * s); + } + } + rv1[l] = 0.0; + rv1[k] = f; + Sigma[k][k] = x; + } + } + std::swap(n,m); + return true; + } + + Matrix(Var * pspace, enumerator pn, enumerator pm) : space(pspace,pspace+pn*pm), n(pn), m(pm) {} + Matrix(enumerator pn, enumerator pm) : space(pn*pm), n(pn), m(pm) {} + Matrix(const Matrix & other) : space(other.n*other.m), n(other.n), m(other.m) + { + for(enumerator i = 0; i < n*m; ++i) + space[i] = other.space[i]; + } + ~Matrix() {} + void Resize(enumerator nrows, enumerator mcols) + { + if( space.size() != mcols*nrows ) + space.resize(mcols*nrows); + n = nrows; + m = mcols; + } + Matrix & operator =(Matrix const & other) + { + space.resize(other.n*other.m); + for(enumerator i = 0; i < other.n*other.m; ++i) + space[i] = other.space[i]; + n = other.n; + m = other.m; + return *this; + } + // i is in [0,n] - row index + // j is in [0,m] - column index + Var & operator()(enumerator i, enumerator j) + { + assert(i >= 0 && i < n); + assert(j >= 0 && j < m); + assert(i*m+j < n*m); //overflow check? + return space[i*m+j]; + } + Var operator()(enumerator i, enumerator j) const + { + assert(i >= 0 && i < n); + assert(j >= 0 && j < m); + assert(i*m+j < n*m); //overflow check? + return space[i*m+j]; + } + Matrix operator-() const + { + Matrix ret(n,m); + for(enumerator k = 0; k < n*m; ++k) ret.space[k] = -space[k]; + return ret; + } + template + Matrix::type> operator-(const Matrix & other) const + { + assert(Rows() == other.Rows()); + assert(Cols() == other.Cols()); + Matrix::type> ret(n,m); //check RVO + for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]-other.space[k]; + return ret; + } + Matrix & operator-=(const Matrix & other) + { + assert(n == other.n); + assert(m == other.m); + for(enumerator k = 0; k < n*m; ++k) space[k] -= other.space[k]; + return *this; + } + template + Matrix::type> operator+(const Matrix & other) const + { + assert(Rows() == other.Rows()); + assert(Cols() == other.Cols()); + Matrix::type> ret(n,m); //check RVO + for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]+other.space[k]; + return ret; + } + Matrix & operator+=(const Matrix & other) + { + assert(n == other.n); + assert(m == other.m); + for(enumerator k = 0; k < n*m; ++k) space[k] += other.space[k]; + return *this; + } + template + Matrix::type> operator*(typeB coef) const + { + Matrix::type> ret(n,m); //check RVO + for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]*coef; + return ret; + } + Matrix & operator*=(Var coef) + { + for(enumerator k = 0; k < n*m; ++k) space[k] *= coef; + return *this; + } + template + Matrix::type> operator/(typeB coef) const + { + Matrix::type> ret(n,m); //check RVO + for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]/coef; + return ret; + } + Matrix & operator/=(Var coef) + { + for(enumerator k = 0; k < n*m; ++k) space[k] /= coef; + return *this; + } + template + Matrix::type> operator*(const Matrix & other) const + { + assert(Cols() == other.Rows()); + Matrix::type> ret(Rows(),other.Cols()); //check RVO + for(enumerator i = 0; i < Rows(); ++i) //loop rows + { + for(enumerator j = 0; j < other.Cols(); ++j) //loop columns + { + typename Promote::type tmp = 0.0; + for(enumerator k = 0; k < Cols(); ++k) + tmp += (*this)(i,k)*other(k,j); + ret(i,j) = tmp; + } + } + return ret; + } + /// performs A*B^{-1} + /// checks existence of B^{-1} only in debug mode. + template + Matrix::type> operator/(const Matrix & other) const + { + std::pair,bool> other_inv = other.Invert(); + assert(other_inv.second); + assert(Cols() == other_inv.Rows()); + Matrix::type> ret(n,other.m); //check RVO + for(enumerator i = 0; i < Rows(); ++i) //loop rows + { + for(enumerator j = 0; j < other_inv.Cols(); ++j) //loop columns + { + typename Promote::type tmp = 0.0; + for(enumerator k = 0; k < Cols(); ++k) + tmp += (*this)(i,k)*other_inv.first(k,j); + ret(i,j) = tmp; + } + } + return ret; + } + Matrix Transpose() const + { + Matrix ret(m,n); + for(enumerator i = 0; i < n; ++i) + { + for(enumerator j = 0; j < m; ++j) + { + ret(j,i) = (*this)(i,j); + } + } + return ret; + } + std::pair Invert() const + { + std::pair ret = std::make_pair(Matrix(m,n),true); + Matrix At = Transpose(); //m by n matrix + Matrix AtB = At*Unit(n); //m by n matrix + Matrix AtA = At*(*this); //m by m matrix + enumerator * order = new enumerator [m]; + for(enumerator i = 0; i < m; ++i) order[i] = i; + for(enumerator i = 0; i < m; i++) + { + enumerator maxk = i, maxq = i, temp2; + Var max, temp; + max = fabs(AtA(maxk,maxq)); + //Find best pivot + for(enumerator k = i; k < m; k++) // over rows + { + for(enumerator q = i; q < m; q++) // over columns + { + if( fabs(AtA(k,q)) > max ) + { + max = fabs(AtA(k,q)); + maxk = k; + maxq = q; + } + } + } + //Exchange rows + if( maxk != i ) + { + for(enumerator q = 0; q < m; q++) // over columns of A + { + temp = AtA(maxk,q); + AtA(maxk,q) = AtA(i,q); + AtA(i,q) = temp; + } + //exchange rhs + for(enumerator q = 0; q < n; q++) // over columns of B + { + temp = AtB(maxk,q); + AtB(maxk,q) = AtB(i,q); + AtB(i,q) = temp; + } + } + //Exchange columns + if( maxq != i ) + { + for(enumerator k = 0; k < m; k++) //over rows + { + temp = AtA(k,maxq); + AtA(k,maxq) = AtA(k,i); + AtA(k,i) = temp; + } + //remember order in sol + { + temp2 = order[maxq]; + order[maxq] = order[i]; + order[i] = temp2; + } + } + + if( fabs(AtA(i,i)) < 1.0e-54 ) + { + bool ok = true; + for(enumerator k = 0; k < n; k++) // over columns of B + { + if( fabs(AtB(i,k)/1.0e-54) > 1 ) + { + ok = false; + break; + } + } + if( ok ) AtA(i,i) = AtA(i,i) < 0.0 ? - 1.0e-12 : 1.0e-12; + else + { + ret.second = false; + return ret; + } + } + for(enumerator k = i+1; k < m; k++) + { + AtA(i,k) /= AtA(i,i); + AtA(k,i) /= AtA(i,i); + } + for(enumerator k = i+1; k < m; k++) + for(enumerator q = i+1; q < m; q++) + { + AtA(k,q) -= AtA(k,i) * AtA(i,i) * AtA(i,q); + } + for(enumerator k = 0; k < n; k++) + { + for(enumerator j = i+1; j < m; j++) //iterate over columns of L + { + AtB(j,k) -= AtB(i,k) * AtA(j,i); + } + AtB(i,k) /= AtA(i,i); + } + } + + for(enumerator k = 0; k < n; k++) + { + for(enumerator i = m; i-- > 0; ) //iterate over rows of U + for(enumerator j = i+1; j < m; j++) + { + AtB(i,k) -= AtB(j,k) * AtA(i,j); + } + for(enumerator i = 0; i < m; i++) + ret.first(order[i],k) = AtB(i,k); + } + delete [] order; + return ret; + } + void Zero() + { + for(enumerator i = 0; i < n*m; ++i) space[i] = 0.0; + } + Var Trace() const + { + assert(n == m); + Var ret = 0.0; + for(enumerator i = 0; i < n; ++i) ret += (*this)(i,i); + return ret; + } + Var * data() {return space.data();} + const Var * data() const {return space.data();} + enumerator Rows() const {return n;} + enumerator Cols() const {return m;} + void Print(INMOST_DATA_REAL_TYPE threshold = 1.0e-10) const + { + for(enumerator k = 0; k < n; ++k) + { + for(enumerator l = 0; l < m; ++l) + { + if( fabs(get_value((*this)(k,l))) > threshold ) + std::cout << (*this)(k,l); + else + std::cout << 0; + std::cout << " "; + } + std::cout << std::endl; + } + } + bool isSymmetric() const + { + if( n != m ) return false; + for(enumerator k = 0; k < n; ++k) + { + for(enumerator l = k+1; l < n; ++l) + if( fabs((*this)(k,l)-(*this)(l,k)) > 1.0e-7 ) + return false; + } + return true; + } + Var DotProduct(const Matrix & other) const + { + assert(n == other.n); + assert(m == other.m); + Var ret = 0.0; + for(enumerator i = 0; i < n*m; ++i) ret += space[i]*other.space[i]; + return ret; + } + Var FrobeniusNorm() + { + Var ret = 0; + for(enumerator i = 0; i < n*m; ++i) ret += space[i]*space[i]; + return sqrt(ret); + } + static Matrix FromTensor(Var * K, enumerator size) + { + Matrix Kc(3,3); + switch(size) + { + case 1: //scalar permeability tensor + Kc.Zero(); + Kc(0,0) = Kc(1,1) = Kc(2,2) = K[0]; + case 3: + Kc.Zero(); //diagonal permeability tensor + Kc(0,0) = K[0]; //KXX + Kc(1,1) = K[1]; //KYY + Kc(2,2) = K[2]; //KZZ + break; + case 6: //symmetric permeability tensor + Kc(0,0) = K[0]; //KXX + Kc(0,1) = Kc(1,0) = K[1]; //KXY + Kc(0,2) = Kc(2,0) = K[2]; //KXZ + Kc(1,1) = K[3]; //KYY + Kc(1,2) = Kc(2,1) = K[4]; //KYZ + Kc(2,2) = K[5]; //KZZ + break; + case 9: //full permeability tensor + Kc(0,0) = K[0]; //KXX + Kc(0,1) = K[1]; //KXY + Kc(0,2) = K[2]; //KXZ + Kc(1,0) = K[3]; //KYX + Kc(1,1) = K[4]; //KYY + Kc(1,2) = K[5]; //KYZ + Kc(2,0) = K[6]; //KZX + Kc(2,1) = K[7]; //KZY + Kc(2,2) = K[8]; //KZZ + break; + } + return Kc; + } + //retrive vector in matrix form from array + static Matrix FromVector(Var * n, enumerator size) + { + return Matrix(n,size,1); + } + //create diagonal matrix from array + static Matrix FromDiagonal(Var * r, enumerator size) + { + Matrix ret(size,size); + ret.Zero(); + for(enumerator k = 0; k < size; ++k) ret(k,k) = r[k]; + return ret; + } + //create diagonal matrix from array of inversed values + static Matrix FromDiagonalInverse(Var * r, enumerator size) + { + Matrix ret(size,size); + ret.Zero(); + for(enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k]; + return ret; + } + //Unit matrix + static Matrix Unit(enumerator pn) + { + Matrix ret(pn,pn); + ret.Zero(); + for(enumerator i = 0; i < pn; ++i) ret(i,i) = 1.0; + return ret; + } + }; + + + typedef Matrix rMatrix; //shortcut for real matrix +#if defined(USE_AUTODIFF) + typedef Matrix vMatrix; //shortcut for matrix with variations +#endif +} + + +#endif //INMOST_DENSE_INCLUDED \ No newline at end of file diff --git a/Source/Headers/inmost_expression.h b/Source/Headers/inmost_expression.h index 7d38bda..dcbe927 100644 --- a/Source/Headers/inmost_expression.h +++ b/Source/Headers/inmost_expression.h @@ -34,9 +34,9 @@ namespace INMOST public: basic_expression() {}//std::cout << this << " Created" << std::endl;} basic_expression(const basic_expression & other) {};//std::cout << this << " Created from " << &other << std::endl;} - __INLINE virtual INMOST_DATA_REAL_TYPE GetValue() const = 0; - __INLINE virtual void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const = 0; - __INLINE virtual void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const = 0; + virtual INMOST_DATA_REAL_TYPE GetValue() const = 0; + virtual void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const = 0; + virtual void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const = 0; //virtual ~basic_expression() {std::cout << this << " Destroied" << std::endl;} }; @@ -393,7 +393,7 @@ namespace INMOST value = arg.GetValue(); value = ::exp(value); } - exp_expression(const exp_expression & b) : unary_expression(b) {} + exp_expression(const exp_expression & b) : arg(b.arg), value(b.value) {} __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }; __INLINE void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const { @@ -797,8 +797,8 @@ namespace INMOST { INMOST_DATA_REAL_TYPE rval = right.GetValue(); value = ::pow(pleft,rval); - if( lval != 0 ) - rdmult = value * ::log(lval); + if( pleft != 0 ) + rdmult = value * ::log(pleft); else rdmult = 0; } @@ -857,131 +857,87 @@ namespace INMOST stencil_expression(const dynarray< std::pair, 64 > & parg) : arg(parg) { value = 0.0; - for(dynarray< std::pair, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) + for(typename dynarray< std::pair, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) value += it->first * it->second.GetValue(); } stencil_expression(const stencil_expression & other) : arg(other.arg), value(other.value) {} __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; } __INLINE void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const { - for(dynarray< std::pair, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) + for(typename dynarray< std::pair, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) it->second.GetDerivative(it->first*mult,r); } __INLINE void GetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const { - for(dynarray< std::pair, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) + for(typename dynarray< std::pair, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) it->second.GetDerivative(it->first*mult,r); } }; - - - - template __INLINE condition_expression condition(shell_expression const & control, shell_expression const & if_ge_zero, shell_expression const & if_lt_zero) { return condition_expression(control,if_ge_zero,if_lt_zero); } - __INLINE INMOST_DATA_REAL_TYPE condition(INMOST_DATA_REAL_TYPE control, INMOST_DATA_REAL_TYPE if_ge_zero, INMOST_DATA_REAL_TYPE if_lt_zero) {return control >= 0.0 ? if_ge_zero : if_lt_zero;} - template __INLINE unary_minus_expression operator-(shell_expression const & Arg) { return unary_minus_expression(Arg); } - template __INLINE abs_expression abs(shell_expression const & Arg) { return abs_expression(Arg); } - template __INLINE exp_expression exp(shell_expression const & Arg) { return exp_expression (Arg); } - template __INLINE log_expression log(shell_expression const & Arg) { return log_expression (Arg); } - template __INLINE sin_expression sin(shell_expression const & Arg) { return sin_expression (Arg ); } - template __INLINE cos_expression cos(shell_expression const & Arg) { return cos_expression (Arg); } - template __INLINE sqrt_expression sqrt(shell_expression const & Arg) { return sqrt_expression (Arg); } - template __INLINE variation_multiplication_expression variation(shell_expression const & Arg,INMOST_DATA_REAL_TYPE Mult) {return variation_multiplication_expression(Arg,Mult);} - __INLINE INMOST_DATA_REAL_TYPE variation(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE) {return Arg;} - template __INLINE INMOST_DATA_REAL_TYPE get_value(shell_expression const & Arg) { return Arg.GetValue(); } - __INLINE INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE Arg) {return Arg;} - template __INLINE soft_abs_expression soft_abs(shell_expression const & Arg, INMOST_DATA_REAL_TYPE tol) { return soft_abs_expression(Arg,tol); } - __INLINE INMOST_DATA_REAL_TYPE soft_abs(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE tol) {return ::sqrt(Arg*Arg+tol*tol);} - template __INLINE soft_sign_expression soft_sign(shell_expression const & Arg, INMOST_DATA_REAL_TYPE tol) { return soft_sign_expression(Arg,tol); } - __INLINE INMOST_DATA_REAL_TYPE soft_sign(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE tol) {return Arg/::sqrt(Arg*Arg+tol*tol);} - - template __INLINE multiplication_expression operator*(shell_expression const & Left, shell_expression const & Right) { return multiplication_expression (Left, Right); } - template __INLINE division_expression operator/(shell_expression const & Left, shell_expression const & Right) { return division_expression (Left, Right); } - template __INLINE addition_expression operator+(shell_expression const & Left, shell_expression const & Right) { return addition_expression (Left, Right); } - template __INLINE subtraction_expression operator-(shell_expression const & Left, shell_expression const & Right) { return subtraction_expression (Left, Right); } - template __INLINE pow_expression pow(shell_expression const & Left, shell_expression const & Right) { return pow_expression (Left, Right); } - template __INLINE soft_max_expression soft_max(shell_expression const & Left, shell_expression const & Right ,INMOST_DATA_REAL_TYPE tol) { return soft_max_expression (Left, Right,tol); } - __INLINE INMOST_DATA_REAL_TYPE soft_max(INMOST_DATA_REAL_TYPE Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol) {return 0.5*(Left+Right+::sqrt((Left-Right)*(Left-Right)+tol*tol));} - template __INLINE soft_min_expression soft_min(shell_expression const & Left, shell_expression const & Right ,INMOST_DATA_REAL_TYPE tol) { return soft_min_expression (Left, Right,tol); } - __INLINE INMOST_DATA_REAL_TYPE soft_min(INMOST_DATA_REAL_TYPE Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol) {return 0.5*(Left+Right-::sqrt((Left-Right)*(Left-Right)+tol*tol));} - - template __INLINE const_pow_expression pow(INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) { return const_pow_expression (Left, Right); } - template __INLINE pow_const_expression pow(shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return pow_const_expression (Left, Right); } - template __INLINE const_multiplication_expression operator*(INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) { return const_multiplication_expression(Right,Left); } - template __INLINE const_multiplication_expression operator*(shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return const_multiplication_expression(Left,Right); } - template __INLINE reciprocal_expression operator/(INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) { return reciprocal_expression(Right,Left); } - template __INLINE const_division_expression operator/(shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return const_division_expression(Left, Right); } - template __INLINE const_addition_expression operator+(INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) { return const_addition_expression(Right,Left); } - template __INLINE const_addition_expression operator+(shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return const_addition_expression(Left,Right); } - template __INLINE const_subtraction_expression operator-(INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) { return const_subtraction_expression(Right, Left); } - template __INLINE const_addition_expression operator-(shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return const_addition_expression(Left, -Right); } - - - template __INLINE bool operator == (shell_expression const & Left, shell_expression const & Right) {return Left.GetValue() == Right.GetValue();} - template __INLINE bool operator != (shell_expression const & Left, shell_expression const & Right) {return Left.GetValue() != Right.GetValue();} - template __INLINE bool operator < (shell_expression const & Left, shell_expression const & Right) {return Left.GetValue() < Right.GetValue();} - template __INLINE bool operator > (shell_expression const & Left, shell_expression const & Right) {return Left.GetValue() > Right.GetValue();} - template __INLINE bool operator <= (shell_expression const & Left, shell_expression const & Right) {return Left.GetValue() <= Right.GetValue();} - template __INLINE bool operator >= (shell_expression const & Left, shell_expression const & Right) {return Left.GetValue() >= Right.GetValue();} - - template __INLINE bool operator == (shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() == Right;} - template __INLINE bool operator != (shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() != Right;} - template __INLINE bool operator < (shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() < Right;} - template __INLINE bool operator > (shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() > Right;} - template __INLINE bool operator <= (shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() <= Right;} - template __INLINE bool operator >= (shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() >= Right;} - - - template __INLINE bool operator == (INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) {return Left == Right.GetValue();} - template __INLINE bool operator != (INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) {return Left != Right.GetValue();} - template __INLINE bool operator < (INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) {return Left < Right.GetValue();} - template __INLINE bool operator > (INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) {return Left > Right.GetValue();} - template __INLINE bool operator <= (INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) {return Left <= Right.GetValue();} - template __INLINE bool operator >= (INMOST_DATA_REAL_TYPE Left, shell_expression const & Right) {return Left >= Right.GetValue();} - - /* - template - class get_value - { - operator INMOST_DATA_REAL_TYPE() = 0; - }; - - template<> - class get_value - { - INMOST_DATA_REAL_TYPE val; - public: - get_value(const INMOST_DATA_REAL_TYPE & v) :val(v) {} - operator INMOST_DATA_REAL_TYPE() {return val;} - }; - - template<> - class get_value - { - INMOST_DATA_REAL_TYPE val; - public: - get_value(const var_expression & v) :val(v.GetValue()) {} - operator INMOST_DATA_REAL_TYPE() {return val;} - }; - - template<> - class get_value - { - INMOST_DATA_REAL_TYPE val; - public: - get_value(const multivar_expression & v) :val(v.GetValue()) {} - operator INMOST_DATA_REAL_TYPE() {return val;} - }; - */ - __INLINE bool check_nans(INMOST_DATA_REAL_TYPE val) {return val != val;} __INLINE bool check_nans(var_expression const & e) {return e.check_nans();} __INLINE bool check_nans(multivar_expression const & e) {return e.check_nans();} - typedef multivar_expression variable; - typedef var_expression unknown; -}; + typedef var_expression unknown; +} + +template __INLINE INMOST::condition_expression condition(INMOST::shell_expression const & control, INMOST::shell_expression const & if_ge_zero, INMOST::shell_expression const & if_lt_zero) { return INMOST::condition_expression(control,if_ge_zero,if_lt_zero); } + __INLINE INMOST_DATA_REAL_TYPE condition(INMOST_DATA_REAL_TYPE control, INMOST_DATA_REAL_TYPE if_ge_zero, INMOST_DATA_REAL_TYPE if_lt_zero) {return control >= 0.0 ? if_ge_zero : if_lt_zero;} +template __INLINE INMOST::unary_minus_expression operator-(INMOST::shell_expression const & Arg) { return INMOST::unary_minus_expression(Arg); } +template __INLINE INMOST::abs_expression abs(INMOST::shell_expression const & Arg) { return INMOST::abs_expression(Arg); } +template __INLINE INMOST::exp_expression exp(INMOST::shell_expression const & Arg) { return INMOST::exp_expression (Arg); } +template __INLINE INMOST::log_expression log(INMOST::shell_expression const & Arg) { return INMOST::log_expression (Arg); } +template __INLINE INMOST::sin_expression sin(INMOST::shell_expression const & Arg) { return INMOST::sin_expression (Arg ); } +template __INLINE INMOST::cos_expression cos(INMOST::shell_expression const & Arg) { return INMOST::cos_expression (Arg); } +template __INLINE INMOST::sqrt_expression sqrt(INMOST::shell_expression const & Arg) { return INMOST::sqrt_expression (Arg); } +template __INLINE INMOST::variation_multiplication_expression variation(INMOST::shell_expression const & Arg,INMOST_DATA_REAL_TYPE Mult) {return INMOST::variation_multiplication_expression(Arg,Mult);} + __INLINE INMOST_DATA_REAL_TYPE variation(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE) {return Arg;} +template __INLINE INMOST_DATA_REAL_TYPE get_value(INMOST::shell_expression const & Arg) { return Arg.GetValue(); } + __INLINE INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE Arg) {return Arg;} +template __INLINE INMOST::soft_abs_expression soft_abs(INMOST::shell_expression const & Arg, INMOST_DATA_REAL_TYPE tol) { return INMOST::soft_abs_expression(Arg,tol); } + __INLINE INMOST_DATA_REAL_TYPE soft_abs(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE tol) {return ::sqrt(Arg*Arg+tol*tol);} +template __INLINE INMOST::soft_sign_expression soft_sign(INMOST::shell_expression const & Arg, INMOST_DATA_REAL_TYPE tol) { return INMOST::soft_sign_expression(Arg,tol); } + __INLINE INMOST_DATA_REAL_TYPE soft_sign(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE tol) {return Arg/::sqrt(Arg*Arg+tol*tol);} +template __INLINE INMOST::multiplication_expression operator*(INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) { return INMOST::multiplication_expression (Left, Right); } +template __INLINE INMOST::division_expression operator/(INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) { return INMOST::division_expression (Left, Right); } +template __INLINE INMOST::addition_expression operator+(INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) { return INMOST::addition_expression (Left, Right); } +template __INLINE INMOST::subtraction_expression operator-(INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) { return INMOST::subtraction_expression (Left, Right); } +template __INLINE INMOST::pow_expression pow(INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) { return INMOST::pow_expression (Left, Right); } +template __INLINE INMOST::soft_max_expression soft_max(INMOST::shell_expression const & Left, INMOST::shell_expression const & Right ,INMOST_DATA_REAL_TYPE tol) { return INMOST::soft_max_expression (Left, Right,tol); } + __INLINE INMOST_DATA_REAL_TYPE soft_max(INMOST_DATA_REAL_TYPE Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol) {return 0.5*(Left+Right+::sqrt((Left-Right)*(Left-Right)+tol*tol));} +template __INLINE INMOST::soft_min_expression soft_min(INMOST::shell_expression const & Left, INMOST::shell_expression const & Right ,INMOST_DATA_REAL_TYPE tol) { return INMOST::soft_min_expression (Left, Right,tol); } + __INLINE INMOST_DATA_REAL_TYPE soft_min(INMOST_DATA_REAL_TYPE Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol) {return 0.5*(Left+Right-::sqrt((Left-Right)*(Left-Right)+tol*tol));} +template __INLINE INMOST::const_pow_expression pow(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) { return INMOST::const_pow_expression (Left, Right); } +template __INLINE INMOST::pow_const_expression pow(INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::pow_const_expression (Left, Right); } +template __INLINE INMOST::const_multiplication_expression operator*(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) { return INMOST::const_multiplication_expression(Right,Left); } +template __INLINE INMOST::const_multiplication_expression operator*(INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::const_multiplication_expression(Left,Right); } +template __INLINE INMOST::reciprocal_expression operator/(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) { return INMOST::reciprocal_expression(Right,Left); } +template __INLINE INMOST::const_division_expression operator/(INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::const_division_expression(Left, Right); } +template __INLINE INMOST::const_addition_expression operator+(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) { return INMOST::const_addition_expression(Right,Left); } +template __INLINE INMOST::const_addition_expression operator+(INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::const_addition_expression(Left,Right); } +template __INLINE INMOST::const_subtraction_expression operator-(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) { return INMOST::const_subtraction_expression(Right, Left); } +template __INLINE INMOST::const_addition_expression operator-(INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::const_addition_expression(Left, -Right); } +template __INLINE bool operator == (INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) {return Left.GetValue() == Right.GetValue();} +template __INLINE bool operator != (INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) {return Left.GetValue() != Right.GetValue();} +template __INLINE bool operator < (INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) {return Left.GetValue() < Right.GetValue();} +template __INLINE bool operator > (INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) {return Left.GetValue() > Right.GetValue();} +template __INLINE bool operator <= (INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) {return Left.GetValue() <= Right.GetValue();} +template __INLINE bool operator >= (INMOST::shell_expression const & Left, INMOST::shell_expression const & Right) {return Left.GetValue() >= Right.GetValue();} +template __INLINE bool operator == (INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() == Right;} +template __INLINE bool operator != (INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() != Right;} +template __INLINE bool operator < (INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() < Right;} +template __INLINE bool operator > (INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() > Right;} +template __INLINE bool operator <= (INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() <= Right;} +template __INLINE bool operator >= (INMOST::shell_expression const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() >= Right;} +template __INLINE bool operator == (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) {return Left == Right.GetValue();} +template __INLINE bool operator != (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) {return Left != Right.GetValue();} +template __INLINE bool operator < (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) {return Left < Right.GetValue();} +template __INLINE bool operator > (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) {return Left > Right.GetValue();} +template __INLINE bool operator <= (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) {return Left <= Right.GetValue();} +template __INLINE bool operator >= (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression const & Right) {return Left >= Right.GetValue();} + + #endif #endif diff --git a/Source/Headers/inmost_mesh.h b/Source/Headers/inmost_mesh.h index f0166d7..ca0c0dc 100644 --- a/Source/Headers/inmost_mesh.h +++ b/Source/Headers/inmost_mesh.h @@ -1851,8 +1851,14 @@ namespace INMOST /// For constant size array returns the same value that may be obtained through GetSize. /// @param h handle of element /// @param tag tag that represents the data - /// @see Mesh::GetSize + /// @see Tag::GetSize INMOST_DATA_ENUM_TYPE GetDataSize (HandleType h,const Tag & tag) const; //For DATA_BULK return number of bytes, otherwise return the length of array + /// Return the size of the structure required to represent the data on current element. + /// This is equal to GetDataSize times Tag::GetBytesSize for all the data types, + /// except for DATA_VARIABLE, that requires a larger structure to accomodate derivatives. + /// @param h handle of element + /// @param tag tag that represents the data + INMOST_DATA_ENUM_TYPE GetDataCapacity (HandleType h,const Tag & tag) const; /// Sets the size of the array for data of variable size. /// If you try to change size of data of constant size then if size is /// different from current then in debug mode (NDEBUG not set) assertion will fail, @@ -2845,6 +2851,7 @@ namespace INMOST bool remember[5][3]; private: void RestoreGeometricTags(); + void RepairGeometricTags(); bool HideGeometricData (GeometricData type, ElementType mask) {return remember[type][ElementNum(mask)-1] = false;} bool ShowGeometricData (GeometricData type, ElementType mask) {return remember[type][ElementNum(mask)-1] = true;} public: @@ -3217,6 +3224,10 @@ namespace INMOST { return GetMeshLink()->GetDataSize(GetHandle(),tag); } + __INLINE INMOST_DATA_ENUM_TYPE Storage::GetDataCapacity(const Tag & tag) const + { + return GetMeshLink()->GetDataCapacity(GetHandle(),tag); + } __INLINE void Storage::SetDataSize(const Tag & tag,INMOST_DATA_ENUM_TYPE new_size) const { GetMeshLink()->SetDataSize(GetHandle(),tag,new_size); diff --git a/Source/Headers/inmost_solver.h b/Source/Headers/inmost_solver.h index 618614b..1a5073b 100644 --- a/Source/Headers/inmost_solver.h +++ b/Source/Headers/inmost_solver.h @@ -58,7 +58,7 @@ namespace INMOST }; static std::string TypeName(Type t); - static INMOST_MPI_Type & GetRowEntryType() {return RowEntryType;} + //solver.cpp:::::::::::::::::::::::::::::::::::::::::::::::::::: public: @@ -347,7 +347,6 @@ namespace INMOST static void Finalize(); static bool isInitialized() {return is_initialized;} static bool isFinalized() {return is_finalized;} - /// Clear all internal data of the current solver including matrix, preconditioner etc. void Clear(); }; diff --git a/Source/Headers/inmost_sparse.h b/Source/Headers/inmost_sparse.h index 1743465..adb4730 100644 --- a/Source/Headers/inmost_sparse.h +++ b/Source/Headers/inmost_sparse.h @@ -14,6 +14,11 @@ namespace INMOST { namespace Sparse { + INMOST_MPI_Type GetRowEntryType(); + void CreateRowEntryType(); + void DestroyRowEntryType(); + bool HaveRowEntryType(); + /// Distributed vector class. /// This class can be used to store both local and distributed dense data of real type. /// For example, to form the right-hand side or initial guess to the solution. diff --git a/Source/Headers/inmost_variable.h b/Source/Headers/inmost_variable.h index c6b7625..b443dd9 100644 --- a/Source/Headers/inmost_variable.h +++ b/Source/Headers/inmost_variable.h @@ -184,17 +184,7 @@ namespace INMOST __INLINE bool check_nans(enhanced_multivar_expression const & e) {return e.check_nans();} - /* - template<> - class get_value - { - INMOST_DATA_REAL_TYPE val; - public: - get_value(const enhanced_multivar_expression & v) :val(v.GetValue()) {} - operator INMOST_DATA_REAL_TYPE() {return val;} - }; - */ - + template class enhance { public: @@ -519,7 +509,7 @@ namespace INMOST { dynarray, 64> tmp; Automatizator::stencil_pairs stncl; - aut.GetStencil(stnclind, e, puser_data, stncl); + aut.GetStencil(stnclind, e, user_data, stncl); tmp.resize(stncl.size()); for(INMOST_DATA_ENUM_TYPE k = 0; k < stncl.size(); ++k) tmp[k] = std::make_pair(stncl[k].first, Arg(stncl[k].second)); @@ -553,7 +543,7 @@ namespace INMOST } unary_pool_expression operator [](const Storage & e) const { - binary_pool pool(Arg(e)); + unary_pool pool(Arg(e)); return unary_pool_expression(pool); } void GetVariation(const Storage & e, Sparse::Row & r) const { (*this)[e].GetDerivative(1.0,r); } @@ -657,42 +647,39 @@ namespace INMOST void GetVariation(const Storage & e, Sparse::Row & r) const { (*this)[e].GetDerivative(1.0,r); } void GetVariation(const Storage & e, Sparse::RowMerger & r) const { (*this)[e].GetDerivative(1.0,r); } }; - - template __INLINE ternary_custom_variable,A,B,C> condition(shell_dynamic_variable const & control, shell_dynamic_variable const & if_ge_zero, shell_dynamic_variable const & if_lt_zero) { return ternary_custom_variable,A,B,C>(control,if_ge_zero,if_lt_zero); } - template __INLINE unary_custom_variable,A> operator-(shell_dynamic_variable const & Arg) { return unary_custom_variable,A>(Arg); } - template __INLINE unary_custom_variable,A> fabs(shell_dynamic_variable const & Arg) { return unary_custom_variable,A>(Arg); } - template __INLINE unary_custom_variable,A> exp(shell_dynamic_variable const & Arg) { return unary_custom_variable,A>(Arg); } - template __INLINE unary_custom_variable,A> log(shell_dynamic_variable const & Arg) { return unary_custom_variable,A>(Arg); } - template __INLINE unary_custom_variable,A> sin(shell_dynamic_variable const & Arg) { return unary_custom_variable,A>(Arg ); } - template __INLINE unary_custom_variable,A> cos(shell_dynamic_variable const & Arg) { return unary_custom_variable,A>(Arg); } - template __INLINE unary_custom_variable,A> sqrt(shell_dynamic_variable const & Arg) { return unary_custom_variable,A>(Arg); } - template __INLINE unary_const_custom_variable,A> variation(shell_dynamic_variable const & Arg, INMOST_DATA_REAL_TYPE Mult) {return unary_const_custom_variable,A>(Arg,Mult);} - - - - template __INLINE binary_custom_variable,A, B> operator+(shell_dynamic_variable const & Left, shell_dynamic_variable const & Right) { return binary_custom_variable,A, B> (Left, Right); } - template __INLINE binary_custom_variable,A, B> operator-(shell_dynamic_variable const & Left, shell_dynamic_variable const & Right) { return binary_custom_variable, A, B> (Left, Right); } - template __INLINE binary_custom_variable,A, B> operator*(shell_dynamic_variable const & Left, shell_dynamic_variable const & Right) { return binary_custom_variable, A, B> (Left, Right); } - template __INLINE binary_custom_variable,A, B> operator/(shell_dynamic_variable const & Left, shell_dynamic_variable const & Right) { return binary_custom_variable, A, B> (Left, Right); } - template __INLINE binary_custom_variable,A, B> pow(shell_dynamic_variable const & Left, shell_dynamic_variable const & Right) { return binary_custom_variable,A, B>(Left, Right); } - - template __INLINE unary_const_custom_variable,B> pow(INMOST_DATA_REAL_TYPE Left, shell_dynamic_variable const & Right) { return unary_const_custom_variable,B>(Left, Right); } - template __INLINE unary_const_custom_variable,A> pow(shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return unary_const_custom_variable,A>(Left, Right); } - template __INLINE unary_const_custom_variable,B> operator*(INMOST_DATA_REAL_TYPE Left, shell_dynamic_variable const & Right) { return unary_const_custom_variable,B>(Right,Left); } - template __INLINE unary_const_custom_variable,A> operator*(shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return unary_const_custom_variable,A>(Left,Right); } - template __INLINE unary_const_custom_variable,B> operator/(INMOST_DATA_REAL_TYPE Left, shell_dynamic_variable const & Right) { return unary_const_custom_variable,B>(Right,Left); } - template __INLINE unary_const_custom_variable,A> operator/(shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return unary_const_custom_variable,A>(Left, Right); } - template __INLINE unary_const_custom_variable,B> operator+(INMOST_DATA_REAL_TYPE Left, shell_dynamic_variable const & Right) { return unary_const_custom_variable,B>(Right,Left); } - template __INLINE unary_const_custom_variable,A> operator+(shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return unary_const_custom_variable,A>(Left,Right); } - template __INLINE unary_const_custom_variable,B> operator-(INMOST_DATA_REAL_TYPE Left, shell_dynamic_variable const & Right) { return unary_const_custom_variable,B>(Right, Left); } - template __INLINE unary_const_custom_variable,A> operator-(shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return unary_const_custom_variable,A>(Left, -Right); } - - template __INLINE stencil_variable stencil(Automatizator & aut, INMOST_DATA_ENUM_TYPE stncl, shell_dynamic_variable const & Arg, void * user_data = NULL) { return stencil_variable(aut,stncl,Arg,user_data); } - typedef abstract_dynamic_variable abstract_variable; typedef enhanced_multivar_expression enhanced_variable; +} + +template __INLINE INMOST::ternary_custom_variable,A,B,C> condition(INMOST::shell_dynamic_variable const & control, INMOST::shell_dynamic_variable const & if_ge_zero, INMOST::shell_dynamic_variable const & if_lt_zero) { return INMOST::ternary_custom_variable,A,B,C>(control,if_ge_zero,if_lt_zero); } +template __INLINE INMOST::unary_custom_variable,A> operator-(INMOST::shell_dynamic_variable const & Arg) { return INMOST::unary_custom_variable,A>(Arg); } +template __INLINE INMOST::unary_custom_variable,A> fabs(INMOST::shell_dynamic_variable const & Arg) { return INMOST::unary_custom_variable,A>(Arg); } +template __INLINE INMOST::unary_custom_variable,A> exp(INMOST::shell_dynamic_variable const & Arg) { return INMOST::unary_custom_variable,A>(Arg); } +template __INLINE INMOST::unary_custom_variable,A> log(INMOST::shell_dynamic_variable const & Arg) { return INMOST::unary_custom_variable,A>(Arg); } +template __INLINE INMOST::unary_custom_variable,A> sin(INMOST::shell_dynamic_variable const & Arg) { return INMOST::unary_custom_variable,A>(Arg ); } +template __INLINE INMOST::unary_custom_variable,A> cos(INMOST::shell_dynamic_variable const & Arg) { return INMOST::unary_custom_variable,A>(Arg); } +template __INLINE INMOST::unary_custom_variable,A> sqrt(INMOST::shell_dynamic_variable const & Arg) { return INMOST::unary_custom_variable,A>(Arg); } +template __INLINE INMOST::unary_const_custom_variable,A> variation(INMOST::shell_dynamic_variable const & Arg, INMOST_DATA_REAL_TYPE Mult) {return INMOST::unary_const_custom_variable,A>(Arg,Mult);} +template __INLINE INMOST::binary_custom_variable,A, B> operator+(INMOST::shell_dynamic_variable const & Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::binary_custom_variable,A, B> (Left, Right); } +template __INLINE INMOST::binary_custom_variable,A, B> operator-(INMOST::shell_dynamic_variable const & Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::binary_custom_variable, A, B> (Left, Right); } +template __INLINE INMOST::binary_custom_variable,A, B> operator*(INMOST::shell_dynamic_variable const & Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::binary_custom_variable, A, B> (Left, Right); } +template __INLINE INMOST::binary_custom_variable,A, B> operator/(INMOST::shell_dynamic_variable const & Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::binary_custom_variable, A, B> (Left, Right); } +template __INLINE INMOST::binary_custom_variable,A, B> pow(INMOST::shell_dynamic_variable const & Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::binary_custom_variable,A, B>(Left, Right); } +template __INLINE INMOST::unary_const_custom_variable,B> pow(INMOST_DATA_REAL_TYPE Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::unary_const_custom_variable,B>(Left, Right); } +template __INLINE INMOST::unary_const_custom_variable,A> pow(INMOST::shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::unary_const_custom_variable,A>(Left, Right); } +template __INLINE INMOST::unary_const_custom_variable,B> operator*(INMOST_DATA_REAL_TYPE Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::unary_const_custom_variable,B>(Right,Left); } +template __INLINE INMOST::unary_const_custom_variable,A> operator*(INMOST::shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::unary_const_custom_variable,A>(Left,Right); } +template __INLINE INMOST::unary_const_custom_variable,B> operator/(INMOST_DATA_REAL_TYPE Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::unary_const_custom_variable,B>(Right,Left); } +template __INLINE INMOST::unary_const_custom_variable,A> operator/(INMOST::shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::unary_const_custom_variable,A>(Left, Right); } +template __INLINE INMOST::unary_const_custom_variable,B> operator+(INMOST_DATA_REAL_TYPE Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::unary_const_custom_variable,B>(Right,Left); } +template __INLINE INMOST::unary_const_custom_variable,A> operator+(INMOST::shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::unary_const_custom_variable,A>(Left,Right); } +template __INLINE INMOST::unary_const_custom_variable,B> operator-(INMOST_DATA_REAL_TYPE Left, INMOST::shell_dynamic_variable const & Right) { return INMOST::unary_const_custom_variable,B>(Right, Left); } +template __INLINE INMOST::unary_const_custom_variable,A> operator-(INMOST::shell_dynamic_variable const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::unary_const_custom_variable,A>(Left, -Right); } +template __INLINE INMOST::stencil_variable stencil(INMOST::Automatizator & aut, INMOST_DATA_ENUM_TYPE stncl, INMOST::shell_dynamic_variable const & Arg, void * user_data = NULL) { return INMOST::stencil_variable(aut,stncl,Arg,user_data); } + -}; + + #endif #endif diff --git a/Source/IO/mesh_xml_file.cpp b/Source/IO/mesh_xml_file.cpp index c6f5f98..e21f91f 100644 --- a/Source/IO/mesh_xml_file.cpp +++ b/Source/IO/mesh_xml_file.cpp @@ -261,7 +261,7 @@ class XMLReader }; Interpreter intrp; std::string src; - std::istream & inp; + std::vector inp; int linebreak, linechar; int hadlinebreak, hadlinechar; int verbose; @@ -295,13 +295,14 @@ class XMLReader EndOfFile, // end of file reached Failure //unexpected error } _state; + std::istream & get_stream() {return *inp.back();} //should not share the reference to the stream with another reader - XMLReader(const XMLReader & other) :src(other.src), inp(other.inp) {} + XMLReader(const XMLReader & other) {} XMLReader & operator =(XMLReader & other) {return *this;} char GetChar() { char c = '\0'; - inp.get(c); + get_stream().get(c); hadlinebreak = linebreak; hadlinechar = linechar; if( c == '\n' ) @@ -318,8 +319,8 @@ class XMLReader { linebreak = hadlinebreak; linechar = hadlinechar; - inp.unget(); - if( inp.fail() ) + get_stream().unget(); + if( get_stream().fail() ) { Report("Stream failed while ungetting the char"); _state = Failure; @@ -417,7 +418,14 @@ public: } std::cout << std::endl; } - XMLReader(std::string sourcename, std::istream & input) :intrp(),src(sourcename),inp(input),linebreak(0),linechar(0),_state(Intro){verbose = 0;} + XMLReader(std::string sourcename, std::istream & input) :intrp(),src(sourcename),linebreak(0),linechar(0),_state(Intro){inp.push_back(&input); verbose = 0;} + void PushStream(std::istream & stream) {inp.push_back(&stream);} + void PopStream() + { + inp.pop_back(); + if( _state == EndOfFile && !inp.empty() ) + _state = Intro; + } //read in or /> skipping for attributes - bool ReadCloseTag() + int ReadCloseTag() { char tmp[2]; tmp[0] = GetChar(); @@ -507,7 +515,7 @@ public: { _state = Intro; if( verbose ) Report("info: closed tag"); - return true; + return 1; //tag was finished with > } else if( tmp[0] == '/' ) //close single stage tag { @@ -516,13 +524,13 @@ public: { _state = Intro; if( verbose ) Report("info: closed tag"); - return true; + return 2; //tag was halted with /> } Report("Encountered '%c%c' while expecting '/>' for tag closing",tmp[0],tmp[1]); } Report("Encountered '%c' while expecting '>' for tag closing",tmp[0]); _state = Failure; - return false; + return 0; } bool isTagFinish() {return _state == ReadCloseTagSlash;} //read or fail @@ -839,8 +847,7 @@ public: { if( ret[ret.size()-2] == ']' && ret[ret.size()-1] == ']' ) { - ret.pop_back(); - ret.pop_back(); + ret.resize(ret.size()-2); RetChar(); //return '>' _state = EndContents; break; @@ -896,13 +903,13 @@ public: case ReadContentsMultiplier: if( isspace(c) ) { - if( ret.back() != '*' ) //maybe user have put a space after the multiplier + if( ret[ret.size()-1] != '*' ) //maybe user have put a space after the multiplier { done = true; _state = WaitContents; } } - else if( c == '(' && ret.back() == '*' ) //expression for the skope + else if( c == '(' && ret[ret.size()-1] == '*' ) //expression for the skope { ret.push_back(c); _state = ReadContentsMultiplierSkopes; @@ -921,9 +928,9 @@ public: case ReadContentsQuotes: if( c == '"' ) { - if( ret.back() == '\\' ) //this was a guarded quote + if( ret[ret.size()-1] == '\\' ) //this was a guarded quote { - ret.pop_back(); //remove guard + ret.resize(ret.size()-1); //remove guard ret.push_back(c); //put skope } else //may have a multiplier after the closing quote @@ -1249,8 +1256,8 @@ namespace INMOST std::vector tags; std::vector new_nodes, new_edges, new_faces, new_cells, new_sets; std::vector set_comparators; - std::string tag, attr, val, save; - int nmeshes = 1; + std::string tag, attr, val, save, incl; + int nmeshes = 1, fintag; tag = reader.ReadOpenTag(); //ParallelMesh if( tag != "ParallelMesh" ) { @@ -1374,6 +1381,7 @@ namespace INMOST } { //Nodes + incl = ""; int nnodes = 0, ndims = 3; tag = reader.ReadOpenTag(); if( tag != "Nodes" ) @@ -1391,45 +1399,59 @@ namespace INMOST ndims = atoi(val.c_str()); if( ndims != 3 || GetDimensions() != ndims ) SetDimensions(ndims); } + else if( attr == "IncludeContents" ) + incl = attr; else reader.Report("Unused attribute for Tags %s='%s'",attr.c_str(),val.c_str()); attr = reader.AttributeName(); } - reader.ReadCloseTag(); + fintag = reader.ReadCloseTag(); new_nodes.reserve(nnodes); - if( reader.ReadOpenContents() ) { - std::vector Vector; - int Repeat; - dynarray xyz; - val = reader.GetContentsWord(); - while(!reader.isContentsEnded() ) + std::fstream incl_stream; + if( !incl.empty() ) { - bool vector = false; - reader.ParseReal(val,Vector,Repeat,nnodes); - for(int l = 0; l < Repeat; ++l) + incl_stream.open(incl.c_str(),std::ios::in); + reader.PushStream(incl_stream); + } + + if( reader.ReadOpenContents() ) + { + std::vector Vector; + int Repeat; + dynarray xyz; + val = reader.GetContentsWord(); + while(!reader.isContentsEnded() ) { - for(int q = 0; q < (int)Vector.size(); ++q) + bool vector = false; + reader.ParseReal(val,Vector,Repeat,nnodes); + for(int l = 0; l < Repeat; ++l) { - xyz.push_back(Vector[q]); - if( xyz.size() == ndims ) + for(int q = 0; q < (int)Vector.size(); ++q) { - new_nodes.push_back(CreateNode(xyz.data())->GetHandle()); - xyz.clear(); + xyz.push_back(Vector[q]); + if( xyz.size() == ndims ) + { + new_nodes.push_back(CreateNode(xyz.data())->GetHandle()); + xyz.clear(); + } } } + val = reader.GetContentsWord(); } - val = reader.GetContentsWord(); + reader.ReadCloseContents(); } - reader.ReadCloseContents(); - } - else - { - reader.Report("Cannot find contents of XML tag"); - throw BadFile; + else + { + reader.Report("Cannot find contents of XML tag"); + throw BadFile; + } + + if( !incl.empty() ) reader.PopStream(); + + if( fintag != 2 ) reader.ReadFinishTag("Nodes"); } - reader.ReadFinishTag("Nodes"); if( new_nodes.size() != nnodes ) { @@ -2048,6 +2070,8 @@ namespace INMOST reader.ReadFinishTag("Mesh"); + RepairGeometricTags(); + if( repair_orientation ) { int numfixed = 0; @@ -2088,8 +2112,8 @@ namespace INMOST if( tags[k].isDefined(etype) ) definition += names[ElementNum(etype)] + ","; if( tags[k].isSparse(etype) ) sparse += names[ElementNum(etype)] + ","; } - if( !definition.empty() ) definition.pop_back(); //remove trailing comma - if( !sparse.empty() ) sparse.pop_back(); //remove trailing comma + if( !definition.empty() ) definition.resize(definition.size()-1); //remove trailing comma + if( !sparse.empty() ) sparse.resize(sparse.size()-1); //remove trailing comma if( sparse == "" ) sparse = "None"; fout << "\t\t\t(static_cast(adata)->size()); case DATA_BULK: return static_cast(static_cast(adata)->size()); case DATA_REFERENCE:return static_cast(static_cast(adata)->size()); +#if defined(USE_AUTODIFF) + case DATA_VARIABLE: return static_cast(static_cast(adata)->size()); +#endif } throw BadTag; } return tag.GetSize(); } + INMOST_DATA_ENUM_TYPE Mesh::GetDataCapacity(HandleType h,const Tag & tag) const + { + assert( tag.GetMeshLink() == this ); + if( tag.GetSize() == ENUMUNDEF ) + { + const void * adata = MGetLink(h,tag); + assert( adata != NULL ); + switch(tag.GetDataType()) + { + case DATA_REAL: return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); + case DATA_INTEGER: return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); + case DATA_BULK: return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); + case DATA_REFERENCE:return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); +#if defined(USE_AUTODIFF) + case DATA_VARIABLE: + { + INMOST_DATA_ENUM_TYPE ret = 0; + const inner_variable_array * arr = static_cast(adata); + ret += arr->size(); + for(inner_variable_array::size_type k = 0; k < arr->size(); ++k) + ret += (*arr)[k].GetRow().Size(); + return ret*sizeof(Sparse::Row::entry_s); + } +#endif + } + throw BadTag; + } +#if defined(USE_AUTODIFF) + if( tag.GetDataType() == DATA_VARIABLE ) + { + INMOST_DATA_ENUM_TYPE ret = 0; + const var * v = static_cast(MGetLink(h,tag)); + return (1+v->GetRow().Size())*sizeof(Sparse::Row::entry_s); + } +#endif + return tag.GetSize()*tag.GetBytesSize(); + } void Mesh::SetDataSize(HandleType h,const Tag & tag,INMOST_DATA_ENUM_TYPE new_size) { assert( tag.GetMeshLink() == this ); @@ -1990,6 +2030,9 @@ namespace INMOST case DATA_INTEGER: static_cast(adata)->resize(new_size); break; case DATA_BULK: static_cast(adata)->resize(new_size); break; case DATA_REFERENCE:static_cast(adata)->resize(new_size); break; +#if defined(USE_AUTODIFF) + case DATA_VARIABLE: static_cast(adata)->resize(new_size); break; +#endif } return; } @@ -2012,8 +2055,51 @@ namespace INMOST case DATA_INTEGER: memcpy(data_out,&(*static_cast(adata))[shift],bytes*size); break; case DATA_BULK: memcpy(data_out,&(*static_cast(adata))[shift],bytes*size); break; case DATA_REFERENCE:memcpy(data_out,&(*static_cast(adata))[shift],bytes*size); break; +#if defined(USE_AUTODIFF) + case DATA_VARIABLE: + { + const inner_variable_array * arr = static_cast(adata); + Sparse::Row::entry_s * data = static_cast(data_out); + int k = 0; + for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) + { + const Sparse::Row & row = (*arr)[r+shift].GetRow(); + data[k].first = row.Size(); + data[k].second = (*arr)[r+shift].GetValue(); + ++k; + for(INMOST_DATA_ENUM_TYPE q = 0; q < row.Size(); ++q) + { + data[k].first = row.GetIndex(q); + data[k].second = row.GetValue(q); + ++k; + } + } + } + break; +#endif } } +#if defined(USE_AUTODIFF) + else if( tag.GetDataType() == DATA_VARIABLE ) + { + Sparse::Row::entry_s * data = static_cast(data_out); + const var * v = static_cast(MGetLink(h,tag)); + int k = 0; + for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) + { + const Sparse::Row & row = v[shift+r].GetRow(); + data[k].first = row.Size(); + data[k].second = v[shift+r].GetValue(); + ++k; + for(INMOST_DATA_ENUM_TYPE r = 0; r < row.Size(); ++r) + { + data[k].first = row.GetIndex(r); + data[k].second = row.GetValue(r); + ++k; + } + } + } +#endif else memcpy(data_out,static_cast(adata)+shift*bytes,size*bytes); return; } @@ -2032,8 +2118,51 @@ namespace INMOST case DATA_INTEGER: memcpy(&(*static_cast(adata))[shift],data_in,bytes*size); break; case DATA_BULK: memcpy(&(*static_cast(adata))[shift],data_in,bytes*size); break; case DATA_REFERENCE:memcpy(&(*static_cast(adata))[shift],data_in,bytes*size); break; +#if defined(USE_AUTODIFF) + case DATA_VARIABLE: + { + inner_variable_array * arr = static_cast(adata); + const Sparse::Row::entry_s * data = static_cast(data_in); + int k = 0; + for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) + { + Sparse::Row & row = (*arr)[r+shift].GetRow(); + row.Resize(data[k].first); + (*arr)[r+shift].SetValue(data[k].second); + ++k; + for(INMOST_DATA_ENUM_TYPE q = 0; q < row.Size(); ++q) + { + row.GetIndex(q) = data[k].first; + row.GetValue(q) = data[k].second; + ++k; + } + } + } + break; +#endif } } +#if defined(USE_AUTODIFF) + else if( tag.GetDataType() == DATA_VARIABLE ) + { + const Sparse::Row::entry_s * data = static_cast(data_in); + var * v = static_cast(MGetLink(h,tag)); + int k = 0; + for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) + { + Sparse::Row & row = v[shift+r].GetRow(); + row.Resize(data[k].first); + v[shift+r].SetValue(data[k].second); + ++k; + for(INMOST_DATA_ENUM_TYPE r = 0; r < row.Size(); ++r) + { + row.GetIndex(r) = data[k].first; + row.GetValue(r) = data[k].second; + ++k; + } + } + } +#endif else memcpy(static_cast(adata)+shift*bytes,data_in,size*bytes); } diff --git a/Source/Mesh/parallel.cpp b/Source/Mesh/parallel.cpp index 7ac9ec4..a300a59 100644 --- a/Source/Mesh/parallel.cpp +++ b/Source/Mesh/parallel.cpp @@ -2111,7 +2111,8 @@ namespace INMOST array_size_send[count]++; INMOST_DATA_ENUM_TYPE s = GetDataSize(*eit,tag); INMOST_DATA_ENUM_TYPE had_s = static_cast(array_data_send.size()); - array_data_send.resize(had_s+s*tag.GetBytesSize()); + //array_data_send.resize(had_s+s*tag.GetBytesSize()); + array_data_send.resize(had_s+GetDataCapacity(*eit,tag)); GetData(*eit,tag,0,s,&array_data_send[had_s]); //REPORT_VAL("size",s); //for(int qq = 0; qq < s; ++qq) @@ -2129,7 +2130,8 @@ namespace INMOST // REPORT_STR("element type " << ElementTypeName(*eit) << " global id " << Integer(*eit,GlobalIDTag())); INMOST_DATA_ENUM_TYPE s = GetDataSize(*eit,tag); INMOST_DATA_ENUM_TYPE had_s = static_cast(array_data_send.size()); - array_data_send.resize(had_s+s*tag.GetBytesSize()); + //array_data_send.resize(had_s+s*tag.GetBytesSize()); + array_data_send.resize(had_s+GetDataCapacity(*eit,tag)); GetData(*eit,tag,0,s,&array_data_send[had_s]); if( size == ENUMUNDEF ) array_size_send.push_back(s); ++total_packed; @@ -2214,7 +2216,7 @@ namespace INMOST // REPORT_VAL("value " << qq, (*(Storage::integer *)&array_data_recv[pos+qq*tag.GetBytesSize()])); //} op(tag,Element(this,*eit),&array_data_recv[pos],array_size_recv[k]); - pos += array_size_recv[k]*tag.GetBytesSize(); + pos += GetDataCapacity(*eit,tag);// array_size_recv[k]*tag.GetBytesSize(); ++k; ++total_unpacked; } @@ -2226,7 +2228,7 @@ namespace INMOST eit = elements[i].begin() + array_size_recv[k++]; assert( !select || GetMarker(*eit,select) ); //if fires then very likely that marker was not synchronized op(tag,Element(this,*eit),&array_data_recv[pos],size); - pos += size*tag.GetBytesSize(); + pos += GetDataCapacity(*eit,tag);//size*tag.GetBytesSize(); ++total_unpacked; } } @@ -2238,7 +2240,7 @@ namespace INMOST for(eit = elements[i].begin(); eit != elements[i].end(); eit++) if( !select || GetMarker(*eit,select) ) { op(tag,Element(this,*eit),&array_data_recv[pos],array_size_recv[k]); - pos += array_size_recv[k]*tag.GetBytesSize(); + pos += GetDataCapacity(*eit,tag);//array_size_recv[k]*tag.GetBytesSize(); ++k; ++total_unpacked; } @@ -2248,7 +2250,7 @@ namespace INMOST for(eit = elements[i].begin(); eit != elements[i].end(); eit++) if( !select || GetMarker(*eit,select) ) { op(tag,Element(this,*eit),&array_data_recv[pos],size); - pos += size*tag.GetBytesSize(); + pos += GetDataCapacity(*eit,tag);//size*tag.GetBytesSize(); ++total_unpacked; } } @@ -2284,7 +2286,12 @@ namespace INMOST std::vector send_size(procs.size(),0), recv_size(procs.size(),0); bool unknown_size = false; - for(unsigned int k = 0; k < tags.size(); k++) if( tags[k].GetSize() == ENUMUNDEF ) unknown_size = true; + for(unsigned int k = 0; k < tags.size(); k++) + if( tags[k].GetSize() == ENUMUNDEF +#if defined(USE_AUTODIFF) + || tags[k].GetDataType() == DATA_VARIABLE +#endif + ) unknown_size = true; //precompute sizes for(p = procs.begin(); p != procs.end(); p++ ) diff --git a/Source/Solver/solver.cpp b/Source/Solver/solver.cpp index ec47e71..8340323 100644 --- a/Source/Solver/solver.cpp +++ b/Source/Solver/solver.cpp @@ -46,8 +46,8 @@ //#define USE_OMP -#define KSOLVER BCGSL_solver -//#define KSOLVER BCGS_solver +//#define KSOLVER BCGSL_solver +#define KSOLVER BCGS_solver //#define ACCELERATED_CONDEST //#define PRINT_CONDEST @@ -63,11 +63,6 @@ namespace INMOST static std::string k3biilu2_database_file = ""; -#if defined(USE_MPI) - INMOST_MPI_Type Solver::RowEntryType = MPI_DATATYPE_NULL; -#else - INMOST_MPI_Type Solver::RowEntryType = 0; -#endif #define GUARD_MPI(x) {ierr = x; if( ierr != MPI_SUCCESS ) {std::cout << #x << " not successfull " << std::endl; MPI_Abort(comm,-1000);}} #define HASH_TABLE_SIZE 2048 @@ -502,7 +497,7 @@ namespace INMOST INMOST_DATA_ENUM_TYPE local_size = 0; for(INMOST_DATA_ENUM_TYPE r = 0; r < vector_exchange_recv[j+1]; r++) local_size += recv_row_sizes[q+r]; - GUARD_MPI(MPI_Irecv(&recv_row_data[f],local_size,Solver::GetRowEntryType(),vector_exchange_recv[j],4*size+vector_exchange_recv[j],comm,&requests[k])); + GUARD_MPI(MPI_Irecv(&recv_row_data[f],local_size,Sparse::GetRowEntryType(),vector_exchange_recv[j],4*size+vector_exchange_recv[j],comm,&requests[k])); q += vector_exchange_recv[j+1]; j += vector_exchange_recv[j+1]+2; f += local_size; @@ -517,7 +512,7 @@ namespace INMOST INMOST_DATA_ENUM_TYPE local_size = 0; for(INMOST_DATA_ENUM_TYPE r = 0; r < vector_exchange_send[j+1]; r++) local_size += send_row_sizes[q+r]; - GUARD_MPI(MPI_Isend(&send_row_data[f],local_size,Solver::GetRowEntryType(),vector_exchange_send[j],4*size+rank,comm,&requests[k+vector_exchange_recv[0]])); + GUARD_MPI(MPI_Isend(&send_row_data[f],local_size,Sparse::GetRowEntryType(),vector_exchange_send[j],4*size+rank,comm,&requests[k+vector_exchange_recv[0]])); q += vector_exchange_send[j+1]; j += vector_exchange_send[j+1]+2; f += local_size; @@ -913,45 +908,13 @@ namespace INMOST #if defined(HAVE_SOLVER_K3BIILU2) SolverInitializeK3biilu2(argc,argv,k3biilu2_database_file.c_str()); #endif -#if defined(USE_MPI) - { - int flag = 0; - int ierr = 0; - MPI_Initialized(&flag); - if( !flag ) - { - ierr = MPI_Init(argc,argv); - if( ierr != MPI_SUCCESS ) - { - std::cout << __FILE__ << ":" << __LINE__ << "problem in MPI_Init" << std::endl; - } - } - MPI_Datatype type[3] = { INMOST_MPI_DATA_ENUM_TYPE, INMOST_MPI_DATA_REAL_TYPE, MPI_UB}; - int blocklen[3] = { 1, 1, 1 }; - MPI_Aint disp[3]; - disp[0] = offsetof(Sparse::Row::entry,first); - disp[1] = offsetof(Sparse::Row::entry,second); - disp[2] = sizeof(Sparse::Row::entry); - ierr = MPI_Type_create_struct(3, blocklen, disp, type, &RowEntryType); - if( ierr != MPI_SUCCESS ) - { - std::cout << __FILE__ << ":" << __LINE__ << "problem in MPI_Type_create_struct" << std::endl; - } - ierr = MPI_Type_commit(&RowEntryType); - if( ierr != MPI_SUCCESS ) - { - std::cout << __FILE__ << ":" << __LINE__ << "problem in MPI_Type_commit" << std::endl; - } - is_initialized = true; - } -#endif + Sparse::CreateRowEntryType(); + is_initialized = true; } void Solver::Finalize() { -#if defined(USE_MPI) - MPI_Type_free(&RowEntryType); -#endif + Sparse::DestroyRowEntryType(); #if defined(USE_SOLVER_PETSC) SolverFinalizePetsc(); #endif @@ -970,10 +933,10 @@ namespace INMOST MPI_Finalized(&flag); if( !flag ) MPI_Finalize(); - is_initialized = false; - is_finalized = true; } #endif + is_initialized = false; + is_finalized = true; } void Solver::SetParameterEnum(std::string name, INMOST_DATA_ENUM_TYPE val) { diff --git a/Source/Solver/solver_bcgsl.hpp b/Source/Solver/solver_bcgsl.hpp index 1519dc4..cc9d823 100644 --- a/Source/Solver/solver_bcgsl.hpp +++ b/Source/Solver/solver_bcgsl.hpp @@ -698,18 +698,6 @@ perturbate: for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) r_tilde[k] /= r_tilde_norm; //Recompute rho1 - /* -#if defined(USE_OMP) -#pragma omp single -#endif - rho1 = 0; -#if defined(USE_OMP) -#pragma omp for reduction(+:rho1) -#endif - for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) - rho1+= r[j][k]*r_tilde[k]; - info->Integrate(&rho1,1); - */ #else reason = "r_tilde perturbation algorithm is disabled"; halt = true; @@ -1411,10 +1399,10 @@ perturbate: bool Solve(Sparse::Vector & RHS, Sparse::Vector & SOL) { assert(isInitialized()); - INMOST_DATA_REAL_TYPE tempa = 0.0, tempb=0.0; + INMOST_DATA_REAL_TYPE tempa = 0.0, tempb=0.0, r0_norm, length; INMOST_DATA_ENUM_TYPE vbeg,vend, vlocbeg, vlocend; INMOST_DATA_INTEGER_TYPE ivbeg,ivend, ivlocbeg, ivlocend; - INMOST_DATA_REAL_TYPE rho = 1, alpha = 1, beta, omega = 1; + INMOST_DATA_REAL_TYPE rho = 1, rho1 = 0, alpha = 1, beta, omega = 1; INMOST_DATA_REAL_TYPE resid0, resid, temp[2]; info->PrepareVector(SOL); info->PrepareVector(RHS); @@ -1470,6 +1458,54 @@ perturbate: while(true) { { + if( false ) + { +perturbate: +#if defined(PERTURBATE_RTILDE) + //std::cout << "Rescue solver by perturbing orthogonal direction!" << std::endl; + //Compute length of the vector +#if defined(USE_OMP) +#pragma omp single +#endif + { + length = static_cast(ivlocend-ivlocbeg); + r0_norm = 0.0; + } + info->Integrate(&length,1); + //perform perturbation (note that rand() with openmp may give identical sequences of random values +#if defined(USE_OMP) +#pragma omp for +#endif + for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) + { + INMOST_DATA_REAL_TYPE unit = 2*static_cast(rand())/static_cast(RAND_MAX)-1.0; + r0[k] += unit/length; + } + //compute norm for orthogonal vector +#if defined(USE_OMP) +#pragma omp for reduction(+:r0_norm) +#endif + for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) + r0_norm += r0[k]*r0[k]; +#if defined(USE_OMP) +#pragma omp single +#endif + { + r0_norm = sqrt(r0_norm); + } + info->Integrate(&r0_norm,1); + //normalize orthogonal vector to unity +#if defined(USE_OMP) +#pragma omp for +#endif + for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) + r0[k] /= r0_norm; + //Recompute rho1 +#else + reason = "r_tilde perturbation algorithm is disabled"; + break; +#endif + } /* if( fabs(rho) < 1.0e-31 ) { @@ -1490,20 +1526,26 @@ perturbate: #endif { beta = 1.0 /rho * alpha / omega; - rho = 0; + rho1 = 0; } #if defined(USE_OMP) #pragma omp for reduction(+:rho) #endif for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) - rho += r0[k]*r[k]; - info->Integrate(&rho,1); + rho1 += r0[k]*r[k]; + info->Integrate(&rho1,1); + if( fabs(rho1) < 1.0e-50 ) + { + //std::cout << "Asking to perturbate r_tilde since rho1 is too small " << rho1 << std::endl; + goto perturbate; + } //info->ScalarProd(r0,r,ivlocbeg,ivlocend,rho); #if defined(USE_OMP) #pragma omp single #endif { - beta *= rho; + beta *= rho1; + rho = rho1; } if( fabs(beta) > 1.0e+100 ) diff --git a/Source/Solver/sparse.cpp b/Source/Solver/sparse.cpp index 48b3373..8997df4 100644 --- a/Source/Solver/sparse.cpp +++ b/Source/Solver/sparse.cpp @@ -7,6 +7,54 @@ namespace INMOST { namespace Sparse { + + bool _hasRowEntryType = false; + INMOST_MPI_Type RowEntryType = INMOST_MPI_DATATYPE_NULL; + + + INMOST_MPI_Type GetRowEntryType() {return RowEntryType;} + + bool HaveRowEntryType() {return _hasRowEntryType;} + + void CreateRowEntryType() + { +#if defined(USE_MPI) + if( HaveRowEntryType() ) + { + int ierr; + MPI_Datatype type[3] = { INMOST_MPI_DATA_ENUM_TYPE, INMOST_MPI_DATA_REAL_TYPE, MPI_UB}; + int blocklen[3] = { 1, 1, 1 }; + MPI_Aint disp[3]; + disp[0] = offsetof(Sparse::Row::entry,first); + disp[1] = offsetof(Sparse::Row::entry,second); + disp[2] = sizeof(Sparse::Row::entry); + ierr = MPI_Type_create_struct(3, blocklen, disp, type, &RowEntryType); + if( ierr != MPI_SUCCESS ) + { + std::cout << __FILE__ << ":" << __LINE__ << "problem in MPI_Type_create_struct" << std::endl; + } + ierr = MPI_Type_commit(&RowEntryType); + if( ierr != MPI_SUCCESS ) + { + std::cout << __FILE__ << ":" << __LINE__ << "problem in MPI_Type_commit" << std::endl; + } + } +#endif + _hasRowEntryType = true; + } + + void DestroyRowEntryType() + { +#if defined(USE_MPI) + if( HaveRowEntryType() ) + { + MPI_Type_free(&RowEntryType); + RowEntryType = INMOST_MPI_DATATYPE_NULL; + } +#endif + _hasRowEntryType = false; + } + RowMerger::RowMerger() : Sorted(true), Nonzeros(0) {} INMOST_DATA_REAL_TYPE & RowMerger::operator[] (INMOST_DATA_ENUM_TYPE pos) -- 2.26.2