diff --git a/CMakeLists.txt b/CMakeLists.txt index 056fb5fca2695400b8fab1124e29af9ded8eae37..ea6a5abfa57d826191645e2702d8511cd9c1e021 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 4e19bef74ece96b015ebfcff98371e6c57aa3ad6..1397294925769f6a15985a12bf14c5ef685ffc6e 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 412087d1a8118952f31eb7acd59ed4bc1305f1fc..b7faf073b6551a0ff29e2be8b475f26a81df2190 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 1246daf2519867fc5a45579c4b57083ff91a5658..0000000000000000000000000000000000000000 --- 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 cadabb741d6b73dc4e152aee3f1aa5d70010ca5e..455452e7fc87cef5958b0ab05a28961a33866a06 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 54cf68edc1131009e533002d69db5c5c7e501092..7e377dc8726a4a9072015a2b1b84fc0cd2bac9ca 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 8f5ed1daf075b365c5220fd279a9d6332b05e2aa..de881f7551de44e7544bd3978dd945770088dfe2 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 e3fd39fb8e53a2ef5f135c44c76f3c83b6487984..acf7b3953fad009d12965fc72b974ac2972b59d8 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 75458d9942905c4e2dfcec2f5e4c72831f97a803..e19f123804286a2108f464f17518f0386c014975 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 83bf44b3c2d1c72f04b2cd16ace01296fc908664..d641ce414db36c7dede3f5797777a5b8ffa57b1e 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 f369ef135d0368edee8c6277123eca861dc5ecca..ff6ebe4b67f2f9879faa1eabb865d265450555c9 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 0000000000000000000000000000000000000000..166aee3a4ab91ef4c5cc828021e5c782c2087f5c --- /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 7d38bdabc69d75efb595d53a92cb974ffeb9b9ea..dcbe92723919d02cd9d1acf033d5c930256afe34 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 f0166d7290635bed89bfe62adc08aef92ec03f10..ca0c0dca7d8c07c96b56db0ac9b3c3f0f788b06e 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 618614bdecca9e13e2e914c1e6c68950e5912ff0..1a5073b8c4aa364f53a3fcf4bce2e99b0790ddeb 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 174346543473a1dc4653db049ab3660cf29a8ea1..adb4730932b80139308790bb03518c8fab7e52aa 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 c6b762567017cb716683b0031722d05bb66b6402..b443dd945579a9a47dc2ec6bed8219f2d445c085 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 c6f5f98f3e300f3213cf323c34e255bcfdba389a..e21f91f14b0c9233787667542dda7b02eee65126 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 7ac9ec439635cb6223f74ecd7e4b8b4f04a657de..a300a596215b2c36999146b90f6c91d3ac673985 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 ec47e71f287d79a5287e8d264d93f0885015c6d0..8340323c55168f91b7f9506c72e1cc2cee7220dc 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 1519dc46cf8b8888d73280ff72125fc6040b9e3a..cc9d823f6204fbb688fc87706825c17ed2f8ac99 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 48b3373b7a9dad1b846ac548230673808e5e5c67..8997df43d39fb0d158e66aca564c3ef03af43377 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)