diff --git a/Examples/ADMFD/main.cpp b/Examples/ADMFD/main.cpp index 0a06d04..785e47d 100644 --- a/Examples/ADMFD/main.cpp +++ b/Examples/ADMFD/main.cpp @@ -199,9 +199,9 @@ int main(int argc,char ** argv) NK(k,k+1,0,3) = n*K; Areas(k,k) = area; } //end of loop over faces - W = NK*(NK.Transpose()*R).Invert(true).first*NK.Transpose(); //stability part - W+=(rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert(true).first*R.Transpose())* - (2.0/(static_cast(NF)*volume)*(NK*K.Invert(true).first*NK.Transpose()).Trace()); + W = NK*(NK.Transpose()*R).PseudoInvert(1.0e-12)*NK.Transpose(); //stability part + W+=(rMatrix::Unit(NF) - R*(R.Transpose()*R).CholeskyInvert()*R.Transpose())* + (2.0/(static_cast(NF)*volume)*(NK*K.CholeskyInvert()*NK.Transpose()).Trace()); W = Areas*W*Areas; //access data structure for gradient matrix in mesh real_array store_W = cell->RealArrayDV(tag_W); diff --git a/Examples/ADVDIFF/stencil.cpp b/Examples/ADVDIFF/stencil.cpp index a188434..5ed7658 100644 --- a/Examples/ADVDIFF/stencil.cpp +++ b/Examples/ADVDIFF/stencil.cpp @@ -249,7 +249,7 @@ bool find_stencils(Cell cK, rK = nF.DotProduct(xF-xK); rN = nF.DotProduct(xN-xF); //correction due to tensor - iQ = -iC * iT.Invert(true).first; + iQ = -iC * iT.Invert(); rQ = nF.DotProduct(iQ); //harmonic point yS = (lambdaC*rN*(xK - iQ - rQ*nF) + lambdaN*(rK+rQ)*xN + rN*(rK+rQ)*nF*KD)/(lambdaC*rN + lambdaN*(rK+rQ)); diff --git a/Examples/GridTools/fix_tiny.cpp b/Examples/GridTools/fix_tiny.cpp index 4bde423..af132d2 100644 --- a/Examples/GridTools/fix_tiny.cpp +++ b/Examples/GridTools/fix_tiny.cpp @@ -105,9 +105,6 @@ static bool BoundingEllipse(Element e, double eps, int iters, rMatrix & Q, rMatr rMatrix M(A.Cols(), A.Cols()); - std::pair inv, iQ; - inv.first.Resize(d + 1, d + 1); - iQ.first.Resize(d, d); int done_iters = 0; double n = d+1, ep , en, kp, kn, beta; int jp,jn; @@ -115,9 +112,7 @@ static bool BoundingEllipse(Element e, double eps, int iters, rMatrix & Q, rMatr { for (int k = 0; k < p.Rows(); ++k) D(k, k) = p(k, 0); - inv = (A1*D*A1.Transpose()).Invert(true); // d+1 by d+1 - if (!inv.second) return false; - M = A1.Transpose()*inv.first*A1; // m by m + M = A1.Transpose()*(A1*D*A1.Transpose()).Invert()*A1; // m by m //std::cout << "matrix M:" << std::endl; //M.Print(); kp = -1.0e20; @@ -167,9 +162,7 @@ static bool BoundingEllipse(Element e, double eps, int iters, rMatrix & Q, rMatr } for (int k = 0; k < p.Rows(); ++k) D(k, k) = p(k, 0); - iQ = (A*(D - p*p.Transpose())*A.Transpose()).Invert(true); - if (!iQ.second) return false; - Q = iQ.first / static_cast(d); + Q = (A*(D - p*p.Transpose())*A.Transpose()).Invert() / static_cast(d); c = A*p; //check if (d == 2) @@ -686,4 +679,4 @@ int main(int argc, char ** argv) } m.Save(grid_out); return 0; -} \ No newline at end of file +} diff --git a/Examples/MFD-ES/main.cpp b/Examples/MFD-ES/main.cpp index 9d64d70..55ae12c 100644 --- a/Examples/MFD-ES/main.cpp +++ b/Examples/MFD-ES/main.cpp @@ -429,8 +429,8 @@ int main(int argc,char ** argv) if( faces[k].BackCell() == c ) tag_i[faces[k]] = k; } //end of loop over faces - W = N*(N.Transpose()*R).Invert(true).first*N.Transpose(); //stability part - W += (rMatrix::Unit(faces.size()) - R*(R.Transpose()*R).Invert(true).first*R.Transpose())*(4.0/(faces.size())*W.Trace()); + W = N*(N.Transpose()*R).Invert()*N.Transpose(); //stability part + W += (rMatrix::Unit(faces.size()) - R*(R.Transpose()*R).Invert()*R.Transpose())*(4.0/(faces.size())*W.Trace()); } //end of loop over cells //initialize normal velocity diff --git a/Source/Autodiff/autodiff.cpp b/Source/Autodiff/autodiff.cpp index ff5f229..ee8ad39 100644 --- a/Source/Autodiff/autodiff.cpp +++ b/Source/Autodiff/autodiff.cpp @@ -18,11 +18,21 @@ namespace INMOST template<> Demote::type AbstractEntry::Access (const Storage& e, INMOST_DATA_ENUM_TYPE pos) const {return Unknown(e,pos);} template<> Demote::type AbstractEntry::Access (const Storage& e, INMOST_DATA_ENUM_TYPE pos) const {return Unknown(e,pos);} - template<> Matrix::type> AbstractEntry::Access (const Storage& e) const {return Value(e);} - template<> Matrix::type> AbstractEntry::Access(const Storage& e) const {return Index(e);} - template<> Matrix::type> AbstractEntry::Access (const Storage& e) const {return Unknown(e);} - template<> Matrix::type> AbstractEntry::Access (const Storage& e) const {return Unknown(e);} - template<> Matrix::type> AbstractEntry::Access (const Storage& e) const {return Unknown(e);} + template<> + Matrix::type, pool_array_t::type> > + AbstractEntry::Access (const Storage& e) const {return Value(e);} + template<> + Matrix::type, pool_array_t::type> > + AbstractEntry::Access(const Storage& e) const {return Index(e);} + template<> + Matrix::type, pool_array_t::type> > + AbstractEntry::Access(const Storage& e) const {return Unknown(e);} + template<> + Matrix::type, pool_array_t::type> > + AbstractEntry::Access(const Storage& e) const {return Unknown(e);} + template<> + Matrix::type,pool_array_t::type> > + AbstractEntry::Access(const Storage& e) const {return Unknown(e);} #if defined(USE_MESH) @@ -392,9 +402,9 @@ namespace INMOST throw Impossible; } - rMatrix MultiEntry::Value(const Storage & e) const + rpMatrix MultiEntry::Value(const Storage & e) const { - vMatrix ret(MatrixSize(e),1); + rpMatrix ret(MatrixSize(e),1); unsigned l = 0, r, t; for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) ) { @@ -405,9 +415,9 @@ namespace INMOST return ret; } - iMatrix MultiEntry::Index(const Storage & e) const + ipMatrix MultiEntry::Index(const Storage & e) const { - iMatrix ret(MatrixSize(e),1); + ipMatrix ret(MatrixSize(e),1); unsigned l = 0, r, t; for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) ) { @@ -418,9 +428,9 @@ namespace INMOST return ret; } - uMatrix MultiEntry::operator [](const Storage & e) const + upMatrix MultiEntry::operator [](const Storage & e) const { - uMatrix ret(MatrixSize(e),1); + upMatrix ret(MatrixSize(e),1); unsigned l = 0, r, t; for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) ) { diff --git a/Source/Headers/container.hpp b/Source/Headers/container.hpp index a551401..34fc4fc 100644 --- a/Source/Headers/container.hpp +++ b/Source/Headers/container.hpp @@ -523,7 +523,7 @@ namespace INMOST replace(begin(),end(),first,last); } template friend class shell; - }; + }; #if defined(PACK_ARRAY) #pragma pack(pop,r1) #endif @@ -2837,29 +2837,366 @@ namespace INMOST const T * operator ->() const {return &item;} }; #endif //_OPENMP - /* + class memory_pool { - std::vector pool; ///< Data storage + static const unsigned pool_size_bits = 16; + //typedef char pool_type[pool_size]; + typedef std::map page_fault_type; + std::vector< char * > pool; ///< Data storage std::vector last_alloc; ///< Last allocated block position, used to track allocation + page_fault_type page_fault; public: - memory_pool() pool_size(0) {} - memory_pool(const memory_pool & b) : pool(b.pool), last_alloc(b.last_alloc) {} - memory_pool & operator = (memory_pool const & b) {pool = b.pool; last_alloc = b.last_alloc; return *this;} + memory_pool(){pool.push_back(new char[1 << pool_size_bits]); last_alloc.push_back(0); } + //memory_pool(const memory_pool & b) : pool(b.pool), last_alloc(b.last_alloc) {} + //memory_pool & operator = (memory_pool const & b) {pool = b.pool; last_alloc = b.last_alloc; return *this;} template - T * allocate(unsigned n, const T & c = T()) + void * allocate(unsigned n, const T & c = T()) { - unsigned oldsize = pool.size(); - pool.resize(pool_size() + sizeof(T)*n); - for(unsigned i = 0; i < n; ++i) - new (&static_cast(static_cast(&pool[oldsize]))[i]) T(c); - last_alloc.push_back(oldsize); - return &pool[oldsize]; + if( sizeof(T)*n < (1 << pool_size_bits) ) + { + unsigned oldpos = last_alloc.size() > 0 ? last_alloc.back() : 0; + unsigned newpos = oldpos + sizeof(T)*n; + unsigned pageold = oldpos >> pool_size_bits; + unsigned pagepos = newpos >> pool_size_bits; + unsigned datapos; + if( pagepos == pool.size() ) + { + //std::cout << "position from " << oldpos << " to " << newpos << " need new page " << pagepos << std::endl; + pool.push_back(new char[1 << pool_size_bits]); + } + + if( pagepos != pageold || last_alloc.empty() ) + { + //std::cout << "add page " << pagepos << " start marker" << std::endl; + last_alloc.push_back(pagepos << pool_size_bits); + oldpos = last_alloc.back(); + newpos = oldpos + sizeof(T)*n; + } + datapos = oldpos%(1<(data)[i]) T(c); + last_alloc.push_back(newpos); + //std::cout << "allocated " << sizeof(T)*n << " bytes from " << oldpos << " to " << newpos << " at page " << pagepos << " at " << data << std::endl; + //std::cout << this << " last_alloc[" << last_alloc.size() << "]:"; + //for(unsigned i = 0; i < last_alloc.size(); ++i) std::cout << " " << last_alloc[i]; + //std::cout << std::endl; + return data; + } + else + { + //T * data = new T[n]; + //for(unsigned i = 0; i < n; ++i) data[i] = c; + void * data = malloc(sizeof(T)*n); + for(unsigned i = 0; i < n; ++i) new (&static_cast(data)[i]) T(c); + page_fault[(void *)data] = n; + //std::cout << "page fault for " << sizeof(T)*n << " bytes allocated at " << data << std::endl; + //std::cout << this << " last_alloc[" << last_alloc.size() << "]:"; + //for(unsigned i = 0; i < last_alloc.size(); ++i) std::cout << " " << last_alloc[i]; + //std::cout << std::endl; + return (void *)data; + } } template - void deallocate + void deallocate(T * mem) + { + unsigned oldpos = last_alloc.size() > 1 ? last_alloc[last_alloc.size()-2] : 0; + unsigned newpos = last_alloc.size() > 0 ? last_alloc[last_alloc.size()-1] : 0; + unsigned pagepos = newpos >> pool_size_bits; + unsigned datapos = oldpos % (1 << pool_size_bits); + //std::cout << "deallocate " << mem << " from " << oldpos << " to " << newpos << " page " << pagepos << std::endl; + //std::cout << this << " last_alloc[" << last_alloc.size() << "]:"; + //for(unsigned i = 0; i < last_alloc.size(); ++i) std::cout << " " << last_alloc[i]; + //std::cout << std::endl; + if( last_alloc.size() > 1 && (((T*)&pool[pagepos][datapos]) == mem) ) + { + unsigned n = (newpos - oldpos)/sizeof(T); + for(unsigned i = 0; i < n; ++i) mem[i].~T(); + last_alloc.pop_back(); + //std::cout << "deallocate from " << oldpos << " to " << newpos << std::endl; + if( last_alloc.back() != 0 && last_alloc.back()%(1<> pool_size_bits) << " start marker " << std::endl; + last_alloc.pop_back(); + } + } + else + { + page_fault_type::iterator it = page_fault.find((void *)mem); + assert(it != page_fault.end() && "deallocated block does not match neither last allocated block nor page fault"); + //if(it != page_fault.end() ) + //{ + // std::cout << "deallocated block does not match neither last allocated block nor page fault"; + // throw -1; + //} + unsigned n = it->second; + //std::cout << "deallocate page fault of " << sizeof(T)*n << " bytes at " << mem << std::endl; + //delete [] mem; + for(unsigned i = 0; i < n; ++i) mem[i].~T(); + free(mem); + page_fault.erase(it); + } + } + ~memory_pool() + { + if( last_alloc.back() != 0 ) std::cout << "warning: memory pool not empty on deallocation!!!" << std::endl; + for(unsigned k = 0; k < pool.size(); ++k) + delete [] pool[k]; + if( !page_fault.empty() ) + { + std::cout << "warning: memory pool's page fault not empty on deallocation!!!" << std::endl; + for(page_fault_type::iterator it = page_fault.begin(); it != page_fault.end(); ++it) + free(it->first); + } + } + }; + + //static thread_private _pool; + memory_pool & get_pool(); + + + + template//, typename enumerator = unsigned int> + class pool_array + { + public: + typedef unsigned size_type; + typedef make_unsigned::type uenum; + template + class _iterator + { + private: + etype * e; + public: + typedef etype * pointer; + typedef etype & reference; + typedef etype value_type; + typedef ptrdiff_t difference_type; + typedef std::random_access_iterator_tag iterator_category; + _iterator():e(NULL){} + _iterator(etype * i):e(i){} + _iterator(const _iterator & other){e = other.e;} + ~_iterator() {}; + _iterator operator -(size_t n) const { return _iterator(e-n); } + _iterator & operator -=(size_t n) { e-=n; return *this; } + _iterator operator +(size_t n) const { return _iterator(e+n); } + _iterator & operator +=(size_t n) { e+=n; return *this; } + _iterator & operator ++(){ ++e; return *this;} + _iterator operator ++(int){ return _iterator(e++); } + _iterator & operator --(){ --e; return *this; } + _iterator operator --(int){ return _iterator(e--); } + ptrdiff_t operator -(const _iterator & other) const {return e-other.e;} + etype & operator *() { return *e; } + const etype & operator *() const { return *e; } + etype * operator ->() { return e; } + _iterator & operator =(_iterator const & other) { e = other.e; return *this; } + bool operator ==(const _iterator & other) const { return e == other.e;} + bool operator !=(const _iterator & other) const { return e != other.e;} + bool operator <(const _iterator & other) const { return e < other.e;} + bool operator >(const _iterator & other) const { return e > other.e;} + bool operator <=(const _iterator & other) const { return e <= other.e;} + bool operator >=(const _iterator & other) const { return e >= other.e;} + operator void *() {return static_cast (e);} + operator const void *() const {return const_cast (e);} + }; + typedef _iterator iterator; + typedef _iterator const_iterator; + template + class _reverse_iterator + { + private: + etype * e; + public: + typedef etype * pointer; + typedef etype & reference; + typedef etype value_type; + typedef ptrdiff_t difference_type; + typedef std::random_access_iterator_tag iterator_category; + _reverse_iterator():e(NULL){} + _reverse_iterator(etype * i):e(i){} + _reverse_iterator(const _reverse_iterator & other){e = other.e;} + ~_reverse_iterator() {}; + _reverse_iterator operator -(size_t n) const { return _reverse_iterator(e+n); } + _reverse_iterator & operator -=(size_t n) { e+=n; return *this; } + _reverse_iterator operator +(size_t n) const {return _reverse_iterator(e-n); } + _reverse_iterator & operator +=(size_t n) { e-=n; return *this; } + _reverse_iterator & operator ++(){ --e; return *this;} + _reverse_iterator operator ++(int){ return _reverse_iterator(e--); } + _reverse_iterator & operator --(){ ++e; return *this; } + _reverse_iterator operator --(int){ return _reverse_iterator(e++); } + ptrdiff_t operator -(const _reverse_iterator & other) const {return other.e-e;} + etype & operator *() { return *e; } + const etype & operator *() const { return *e; } + etype * operator ->() { return e; } + _reverse_iterator & operator =(_reverse_iterator const & other) { e = other.e; return *this;} + bool operator ==(const _reverse_iterator & other) const { return e == other.e;} + bool operator !=(const _reverse_iterator & other) const { return e != other.e;} + bool operator <(const _reverse_iterator & other) const { return e < other.e;} + bool operator >(const _reverse_iterator & other) const { return e > other.e;} + bool operator <=(const _reverse_iterator & other) const { return e <= other.e;} + bool operator >=(const _reverse_iterator & other) const { return e >= other.e;} + operator void *() {return static_cast (e);} + operator const void *() const {return static_cast (e);} + }; + typedef _reverse_iterator reverse_iterator; + typedef _reverse_iterator const_reverse_iterator; + private: + + element * m_arr; + size_type m_size; + public: + __INLINE element * data() {return m_arr;} + __INLINE const element * data() const {return m_arr;} + pool_array() + { + m_size = 0; + m_arr = NULL; + } + pool_array(size_type n,element c = element()) + { + m_size = n; + m_arr = (element *)get_pool().allocate(n,c); + assert(m_arr != NULL); + } + template + pool_array(InputIterator first, InputIterator last) + { + //isInputForwardIterators(); + m_size = static_cast(std::distance(first,last)); + m_arr = (element *)get_pool().allocate(m_size,element()); + assert(m_arr != NULL); + { + size_type i = 0; + InputIterator it = first; + while(it != last) new (m_arr+i++) element(*it++); + } + } + pool_array(const pool_array & other) + { + m_size = other.m_size; + if( m_size ) + { + m_arr = (element *)get_pool().allocate(m_size,element()); + assert(m_arr != NULL); + } + else m_arr = NULL; + for(size_type i = 0; i < m_size; i++) new (m_arr+i) element(other.m_arr[i]); + } + ~pool_array() + { + get_pool().deallocate(m_arr); + m_arr = NULL; + m_size = 0; + } + __INLINE const element & operator [] (size_type n) const + { + assert(n < m_size); + return m_arr[n]; + } + __INLINE element & operator [] (size_type n) + { + assert(n < m_size); + return m_arr[n]; + } + __INLINE const element & at (size_type n) const + { + assert(n < m_size); + return m_arr[n]; + } + __INLINE element & at (size_type n) + { + assert(n < m_size); + return m_arr[n]; + } + /* + //think what todo here + pool_array & operator =(pool_array const & other) + { + if( this != &other ) + { + for(size_type i = 0; i < m_size; i++) m_arr[i].~element(); + if(m_arr != NULL ) + { + free(m_arr); + m_arr = NULL; + m_size = 0; + } + if( other.m_arr != NULL ) + { + m_size = other.m_size; + m_arr = static_cast(malloc(sizeof(element)*growth_formula(m_size))); + assert(m_arr != NULL); + memcpy(m_arr,other.m_arr,sizeof(element)*m_size); + } + } + return *this; + } + */ + void resize(size_type n, element c = element()) + { + if( m_size == 0 && m_arr == NULL ) + { + m_arr = (element *)get_pool().allocate(n,c); + assert(m_arr != NULL); + m_size = n; + } + else if( n != m_size ) + { + assert(false && "resize of non-empty pool_array"); + } + } + __INLINE element & back() + { + assert(m_arr != NULL); + assert(m_size > 0); + return m_arr[m_size-1]; + } + __INLINE const element & back() const + { + assert(m_arr != NULL); + assert(m_size > 0); + return m_arr[m_size-1]; + } + __INLINE element & front() + { + assert(m_arr != NULL); + assert(m_size > 0); + return m_arr[0]; + } + __INLINE const element & front() const + { + assert(m_arr != NULL); + assert(m_size > 0); + return m_arr[0]; + } + __INLINE size_type capacity() { return m_size; } + __INLINE bool empty() const { if( m_size ) return false; return true; } + __INLINE size_type size() const {return m_size;} + __INLINE size_type capacity() const {return m_size;} + void clear() + { + get_pool().deallocate(m_arr); + m_arr = NULL; + m_size = 0; + } + void swap(pool_array & other) + { + element * t_m_arr = m_arr; + size_type t_m_size = m_size; + m_arr = other.m_arr; + m_size = other.m_size; + other.m_arr = t_m_arr; + other.m_size = t_m_size; + } + __INLINE iterator begin() { return m_arr; } + __INLINE iterator end() { return m_arr+m_size; } + __INLINE const_iterator begin() const { return m_arr; } + __INLINE const_iterator end() const { return m_arr+m_size; } + __INLINE reverse_iterator rbegin() { return reverse_iterator(m_arr+m_size-1); } + __INLINE reverse_iterator rend() { return reverse_iterator(m_arr-1); } + __INLINE const_reverse_iterator rbegin() const { return const_reverse_iterator(m_arr+m_size-1); } + __INLINE const_reverse_iterator rend() const { return const_reverse_iterator(m_arr-1); } }; - */ } #endif diff --git a/Source/Headers/inmost_autodiff.h b/Source/Headers/inmost_autodiff.h index 489eb06..acfb631 100644 --- a/Source/Headers/inmost_autodiff.h +++ b/Source/Headers/inmost_autodiff.h @@ -63,20 +63,24 @@ namespace INMOST /// @param pos Position for which to extract the unknown, should be no larger then MatrixSize. virtual unknown Unknown(const Storage & e, INMOST_DATA_ENUM_TYPE pos) const = 0; /// Return vector filled with values of unknowns of the block. - virtual rMatrix Value(const Storage & e) const = 0; + virtual rpMatrix Value(const Storage & e) const = 0; /// Return vector filled with indices of unknowns of the block. - virtual iMatrix Index(const Storage & e) const = 0; + virtual ipMatrix Index(const Storage & e) const = 0; /// Return vector filled with unknowns of the block with their derivatives. - virtual uMatrix Unknown(const Storage & e) const = 0; + virtual upMatrix Unknown(const Storage & e) const = 0; /// Return vector filled with either values or indices or unknowns of the block, /// depending on the template parameter. - template Matrix::type> Access(const Storage &e) const; + template + Matrix::type, pool_array_t::type> > + Access(const Storage &e) const; /// Return either value or index or unknown at specified position of the block, /// depending on the template parameter. /// @param pos Position in the block, should be no larger then MatrixSize. - template typename Demote::type Access(const Storage &e, INMOST_DATA_ENUM_TYPE pos) const; + template + typename Demote::type + Access(const Storage &e, INMOST_DATA_ENUM_TYPE pos) const; /// Return vector filled with unknowns of the block with their derivatives. - virtual uMatrix operator [](const Storage & e) const = 0; + virtual upMatrix operator [](const Storage & e) const = 0; /// The intended size of the matrix for this entry. virtual INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const = 0; /// Number of tags in block. @@ -131,13 +135,13 @@ namespace INMOST /// Return unknown in vector of variables of the block at certain position. unknown Unknown(const Storage & e, INMOST_DATA_ENUM_TYPE pos) const {return unknown(Value(e,pos),Index(e,pos));} /// Return vector filled with values of unknowns of the block. - rMatrix Value(const Storage & e) const {rMatrix ret(MatrixSize(e),1); for(unsigned k = 0; k < Size(); ++k) ret(k,0) = Value(e,k); return ret; } + rpMatrix Value(const Storage & e) const {rpMatrix ret(MatrixSize(e),1); for(unsigned k = 0; k < Size(); ++k) ret(k,0) = Value(e,k); return ret; } /// Return vector filled with indices of unknowns of the block. - iMatrix Index(const Storage & e) const {iMatrix ret(MatrixSize(e),1); for(unsigned k = 0; k < Size(); ++k) ret(k,0) = Index(e,k); return ret; } + ipMatrix Index(const Storage & e) const {ipMatrix ret(MatrixSize(e),1); for(unsigned k = 0; k < Size(); ++k) ret(k,0) = Index(e,k); return ret; } /// Return vector filled with unknowns of the block with their derivatives. - uMatrix Unknown(const Storage & e) const {return BlockEntry::operator [](e);} + upMatrix Unknown(const Storage & e) const {return BlockEntry::operator [](e);} /// Return vector filled with unknowns of the block with their derivatives. - uMatrix operator [](const Storage & e) const {uMatrix ret(MatrixSize(e),1); for(unsigned k = 0; k < Size(); ++k) ret(k,0) = Unknown(e,k); return ret; } + upMatrix operator [](const Storage & e) const {upMatrix ret(MatrixSize(e),1); for(unsigned k = 0; k < Size(); ++k) ret(k,0) = Unknown(e,k); return ret; } /// The intended size of the matrix for this entry. INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const {(void)e; return Size();} /// Number of tags in block. @@ -174,13 +178,13 @@ namespace INMOST /// Return unknown in vector of variables of the block at certain position. unknown Unknown(const Storage & e, INMOST_DATA_ENUM_TYPE pos) const {assert(pos==0); return unknown(Value(e,pos),Index(e,pos));} /// Return vector filled with values of unknowns of the block. - rMatrix Value(const Storage & e) const { rMatrix ret(1,1); ret(0,0) = Value(e,0); return ret; } + rpMatrix Value(const Storage & e) const { rpMatrix ret(1,1); ret(0,0) = Value(e,0); return ret; } /// Return vector filled with indices of unknowns of the block. - iMatrix Index(const Storage & e) const { iMatrix ret(1,1); ret(0,0) = Index(e,0); return ret; } + ipMatrix Index(const Storage & e) const { ipMatrix ret(1,1); ret(0,0) = Index(e,0); return ret; } /// Return vector filled with unknowns of the block with their derivatives. - uMatrix Unknown(const Storage & e) const {return SingleEntry::operator [](e);} + upMatrix Unknown(const Storage & e) const {return SingleEntry::operator [](e);} /// Return vector filled with unknowns of the block with their derivatives. - uMatrix operator [](const Storage & e) const { uMatrix ret(1,1); ret(0,0) = Unknown(e,0); return ret; } + upMatrix operator [](const Storage & e) const { upMatrix ret(1,1); ret(0,0) = Unknown(e,0); return ret; } /// The intended size of the matrix for this entry. INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const {(void)e; return 1;} /// Number of tags in block. @@ -216,13 +220,13 @@ namespace INMOST /// Return unknown in vector of variables of the block at certain position. unknown Unknown(const Storage & e, INMOST_DATA_ENUM_TYPE pos) const {assert(pos(new division_block_variable(*this));} }; @@ -694,11 +694,11 @@ namespace INMOST } /// Get value of variable expression on provided element e. rMatrix Value(const Storage & e) const - {return ArgA.Value(e).Invert(true).first;} + {return ArgA.Value(e).Invert();} /// Get value with derivatives of variable expression on provided element e. /// This function collapses associated expression tree into multivar_expression. vMatrix Variable(const Storage & e) const - {return ArgA.Variable(e).Invert(true).first;} + {return ArgA.Variable(e).Invert();} /// Make a copy of this class, used to reproduce and store a tree of variable expressions. abstract_dynamic_block_variable * Copy() const {return static_cast(new inverse_block_variable(*this));} }; @@ -725,11 +725,11 @@ namespace INMOST } /// Get value of variable expression on provided element e. rMatrix Value(const Storage & e) const - {return ArgA.Value(e).PseudoInvert(eps,true).first;} + {return ArgA.Value(e).PseudoInvert(eps);} /// Get value with derivatives of variable expression on provided element e. /// This function collapses associated expression tree into multivar_expression. vMatrix Variable(const Storage & e) const - {return ArgA.Variable(e).PseudoInvert(eps,true).first;} + {return ArgA.Variable(e).PseudoInvert(eps);} /// Make a copy of this class, used to reproduce and store a tree of variable expressions. abstract_dynamic_block_variable * Copy() const {return static_cast(new pseudo_inverse_block_variable(*this));} }; diff --git a/Source/Headers/inmost_common.h b/Source/Headers/inmost_common.h index 284730a..cb4a858 100644 --- a/Source/Headers/inmost_common.h +++ b/Source/Headers/inmost_common.h @@ -259,6 +259,12 @@ namespace INMOST UnknownWeightSize, DistributionTagWasNotFilled, + /// The list of errros that may occur in linear algebra + MatrixError = 600, + MatrixCholeskySolveFail, + MatrixSolveFail, + MatrixPseudoSolveFail, + /// The very tail of the errors list. NotImplemented = 1000, Impossible @@ -287,5 +293,11 @@ namespace INMOST class SymmetricMatrix; } +#include + +static int __isnan__(double x) { return x != x; } +static int __isinf__(double x) { return fabs(x) > DBL_MAX; } +static int __isbad(double x) { return __isnan__(x) || __isinf__(x); } + #endif //INMOST_COMMON_INCLUDED diff --git a/Source/Headers/inmost_dense.h b/Source/Headers/inmost_dense.h index 2956d0d..52b2d42 100644 --- a/Source/Headers/inmost_dense.h +++ b/Source/Headers/inmost_dense.h @@ -4,6 +4,8 @@ #include "inmost_expression.h" #include +#define pool_array_t pool_array +//#define pool_array_t array namespace INMOST { @@ -221,24 +223,35 @@ namespace INMOST Var c, f, h, s, x, y, z; Var g = 0.0, scale = 0.0; INMOST_DATA_REAL_TYPE anorm = 0.0; + if (n >= m) + { + U = (*this); + Sigma.Resize(m,m); + Sigma.Zero(); + V.Resize(m,m); + } + else // n < m + { + U.Resize(n,n); + Sigma.Resize(n,n); + Sigma.Zero(); + V.Resize(m,n); + } if (n < m) { - bool success = Transpose().SVD(U,Sigma,V); + bool success = Transpose().SVD(V,Sigma,U); if( success ) { - U.Swap(V); - U = U.Transpose(); - V = V.Transpose(); + //U.Swap(V); + //U = U.Transpose(); + //V = V.Transpose(); return true; } else return false; } //m <= n - array _rv1(m); - shell rv1(_rv1); - U = (*this); - Sigma.Resize(m,m); - Sigma.Zero(); - V.Resize(m,m); + dynarray rv1(m); + //array _rv1(m); + //shell rv1(_rv1); std::swap(n,m); //this how original algorithm takes it // Householder reduction to bidiagonal form for (i = 0; i < n; i++) @@ -512,16 +525,32 @@ namespace INMOST } /// Transpose the current matrix. /// @return Transposed matrix. - Matrix Transpose() const; + Matrix > Transpose() const; ///Exchange contents of two matrices. virtual void Swap(AbstractMatrix & b); /// Unary minus. Change sign of each element of the matrix. - Matrix operator-() const; + Matrix > operator-() const; + /// Cross-product operation for a vector. + /// Both right hand side and left hand side should be a vector + /// @param other The right hand side of cross product. + /// @return The cross product of current and right hand side vector. + /// \warning Works for 3x1 vector and 3xm m-vectors as right hand side. + template + Matrix::type, pool_array_t::type> > + CrossProduct(const AbstractMatrix & other) const; + /// Transformation matrix from current vector to provided vector. + /// @param other Vector to transform to. + /// @return A sqaure (rotation) matrix that transforms current vector into right hand side vector. + /// \warning Works only for 3x1 vectors. + template + Matrix::type, pool_array_t::type> > + Transform(const AbstractMatrix & other) const; /// Subtract a matrix. /// @param other The matrix to be subtracted. /// @return Difference of current matrix with another matrix. template - Matrix::type> operator-(const AbstractMatrix & other) const; + Matrix::type, pool_array_t::type> > + operator-(const AbstractMatrix & other) const; /// Subtract a matrix and store result in the current. /// @param other The matrix to be subtracted. /// @return Reference to the current matrix. @@ -531,7 +560,8 @@ namespace INMOST /// @param other The matrix to be added. /// @return Sum of current matrix with another matrix. template - Matrix::type> operator+(const AbstractMatrix & other) const; + Matrix::type, pool_array_t::type> > + operator+(const AbstractMatrix & other) const; /// Add a matrix and store result in the current. /// @param other The matrix to be added. /// @return Reference to the current matrix. @@ -541,7 +571,8 @@ namespace INMOST /// @param other Matrix to be multiplied from right. /// @return Matrix multipled by another matrix. template - Matrix::type> operator*(const AbstractMatrix & other) const; + Matrix::type, pool_array_t::type> > + operator*(const AbstractMatrix & other) const; /// Multiply matrix with another matrix in-place. /// @param B Another matrix to the right in multiplication. /// @return Reference to current matrix. @@ -551,7 +582,8 @@ namespace INMOST /// @param coef Coefficient. /// @return Matrix multiplied by the coefficient. template - Matrix::type> operator*(typeB coef) const; + Matrix::type, pool_array_t::type> > + operator*(typeB coef) const; /// Multiply the matrix by the coefficient of the same type and store the result. /// @param coef Coefficient. /// @return Reference to the current matrix. @@ -561,7 +593,7 @@ namespace INMOST /// @param coef Coefficient. /// @return Matrix divided by the coefficient. template - Matrix::type> + Matrix::type, pool_array_t::type> > operator/(typeB coef) const; /// Divide the matrix by the coefficient of the same type and store the result. /// @param coef Coefficient. @@ -570,42 +602,49 @@ namespace INMOST AbstractMatrix & operator/=(typeB coef); /// Performs B^{-1}*A, multiplication by inverse matrix from left. + /// Throws exception if matrix is not invertable. See Mesh::PseudoSolve for + /// singular matrices. /// @param other Matrix to be inverted and multiplied from left. /// @return Multiplication of current matrix by inverse of another - /// matrix from left and boolean. - /// If boolean is true, then the matrix was inverted successfully. - /// \todo (test) Use Solve here. + /// @see Matrix::PseudoInvert. template - std::pair::type>,bool> + Matrix::type, pool_array_t::type> > operator/(const AbstractMatrix & other) const { - return other.Solve(*this); + Matrix::type, pool_array_t::type> > ret(other.Cols(),Cols()); + ret = other.Solve(*this); + return ret; } /// Kronecker product, latex symbol \otimes. /// @param other Matrix on the right of the Kronecker product. /// @return Result of the Kronecker product of current and another matrix. template - Matrix::type> + Matrix::type, pool_array_t::type> > Kronecker(const AbstractMatrix & other) const; /// Inverts matrix using Crout-LU decomposition with full pivoting for /// maximum element. Works well for non-singular matrices, for singular /// matrices look see Matrix::PseudoInvert. - /// @param print_fail Verbose output for singular matrices. - /// @return A pair of inverse matrix and boolean. If boolean is true, - /// then the matrix was inverted successfully. + /// @param ierr Returns error on fail. If ierr is NULL, then throws an exception. + /// If *ierr == -1 on input, then prints out information in case of failure. + /// In case of failure *ierr equal to a positive integer that represents the row + /// on which the failure occured (starting from 1), in case of no failure *ierr = 0. + /// @return Inverse matrix. /// @see Matrix::PseudoInvert. /// \todo (test) Activate and test implementation with Solve. - std::pair,bool> Invert(bool print_fail = false) const; + Matrix > Invert(int * ierr = NULL) const; /// Inverts symmetric positive-definite matrix using Cholesky decomposition. - std::pair,bool> CholeskyInvert(bool print_fail = false) const; + Matrix > CholeskyInvert(int * ierr = NULL) const; /// Finds X in A*X=B, where A and B are general matrices. /// Converts system into A^T*A*X=A^T*B. /// Inverts matrix A^T*A using Crout-LU decomposition with full pivoting for /// maximum element. Works well for non-singular matrices, for singular /// matrices see Matrix::PseudoInvert. - /// @param print_fail Verbose output for singular matrices. - /// @return A pair of inverse matrix and boolean. If boolean is true, - /// then the matrix was inverted successfully. + /// @param B Right hand side matrix. + /// @param ierr Returns error on fail. If ierr is NULL, then throws an exception. + /// If *ierr == -1 on input, then prints out information in case of failure. + /// In case of failure *ierr equal to a positive integer that represents the row + /// on which the failure occured (starting from 1), in case of no failure *ierr = 0. + /// @return Inverse matrix, /// @see Matrix::PseudoInvert. /// \warning Number of rows in matrices A and B should match. /// \todo @@ -613,13 +652,20 @@ namespace INMOST /// 2. Maximum product transversal + block pivoting instead of pivoting /// by maximum element. template - std::pair::type>,bool> - Solve(const AbstractMatrix & B, bool print_fail = false) const; + Matrix::type, pool_array_t::type> > + Solve(const AbstractMatrix & B, int * ierr = NULL) const; /// Finds X in A*X=B, where A is a square symmetric positive definite matrix. /// Uses Cholesky decomposition algorithm. + /// @param B Right hand side matrix. + /// @param ierr Returns error on fail. If ierr is NULL, then throws an exception. + /// If *ierr == -1 on input, then prints out information in case of failure. + /// In case of failure *ierr equal to a positive integer that represents the row + /// on which the failure occured (starting from 1), in case of no failure *ierr = 0. + /// @see Matrix::PseudoInvert. + /// @return Inverse matrix, template - std::pair::type>,bool> - CholeskySolve(const AbstractMatrix & B, bool print_fail = false) const; + Matrix::type, pool_array_t::type> > + CholeskySolve(const AbstractMatrix & B, int * ierr = NULL) const; /// Calculate sum of the diagonal elements of the matrix. /// @return Trace of the matrix. Var Trace() const @@ -638,7 +684,7 @@ namespace INMOST { for(enumerator l = 0; l < Cols(); ++l) { - if( std::isinf(get_value((*this)(k,l))) ) + if( __isinf__(get_value((*this)(k,l))) ) std::cout << std::setw(12) << "inf"; else if( std::isnan(get_value((*this)(k,l))) ) std::cout << std::setw(12) << "nan"; @@ -701,24 +747,26 @@ namespace INMOST } /// Calculates Moore-Penrose pseudo-inverse of the matrix. /// @param tol Thershold for singular values. Singular values smaller - /// then threshold are considered to be zero. - /// @param print_fail Verbose output if singular value decomposition - /// algorithm has failed. + /// then threshold are considered to be zero. + /// @param ierr Returns error on fail. If ierr is NULL, then throws an exception. + /// If *ierr == -1 on input, then prints out information in case of failure. + /// In case of failure *ierr = 1, in case of no failure *ierr = 0. /// @return A pair of pseudo-inverse matrix and boolean. If boolean is true, /// then the matrix was inverted successfully. - std::pair,bool> PseudoInvert(INMOST_DATA_REAL_TYPE tol = 0, bool print_fail = false) const; + Matrix > PseudoInvert(INMOST_DATA_REAL_TYPE tol = 0, int * ierr = NULL) const; /// Solves the system of equations of the form A*X=B, with A and B matrices. /// Uses Moore-Penrose pseudo-inverse of the matrix A and calculates X = A^+*B. /// @param B Matrix at the right hand side. /// @param tol Thershold for singular values. Singular values smaller - /// then threshold are considered to be zero. - /// @param print_fail Verbose output if singular value decomposition - /// algorithm has failed. + /// then threshold are considered to be zero. + /// @param ierr Returns error on fail. If ierr is NULL, then throws an exception. + /// If *ierr == -1 on input, then prints out information in case of failure. + /// In case of failure *ierr = 1, in case of no failure *ierr = 0. /// @return A pair of the solution matrix X and boolean. If boolean is true, /// then the matrix was inverted successfully. template - std::pair::type>,bool> - PseudoSolve(const AbstractMatrix & B, INMOST_DATA_REAL_TYPE tol = 0, bool print_fail = false) const; + Matrix::type, pool_array_t::type> > + PseudoSolve(const AbstractMatrix & B, INMOST_DATA_REAL_TYPE tol = 0, int * ierr = NULL) const; /// Extract submatrix of a matrix. /// Let A = {a_ij}, i in [0,n), j in [0,m) be original matrix. /// Then the method returns B = {a_ij}, i in [ibeg,iend), @@ -728,12 +776,12 @@ namespace INMOST /// @param jbeg Starting column in the original matrix. /// @param jend Last column (excluded) in the original matrix. /// @return Submatrix of the original matrix. - Matrix ExtractSubMatrix(enumerator ibeg, enumerator iend, enumerator jbeg, enumerator jend) const; + Matrix > ExtractSubMatrix(enumerator ibeg, enumerator iend, enumerator jbeg, enumerator jend) const; /// Change representation of the matrix into matrix of another size. /// Useful to change representation from matrix into vector and back. /// Replaces original number of columns and rows with a new one. /// @return Matrix with same entries and provided number of rows and columns. - Matrix Repack(enumerator rows, enumerator cols) const; + Matrix > Repack(enumerator rows, enumerator cols) const; }; @@ -919,9 +967,9 @@ namespace INMOST /// @param size Size of the input array. /// @param matsize Size of the final tensor. /// @return Matrix of the tensor of size matsize by matsize. - static SymmetricMatrix FromTensor(const Var * K, enumerator size, enumerator matsize = 3) + static SymmetricMatrix > FromTensor(const Var * K, enumerator size, enumerator matsize = 3) { - Matrix Kc(matsize,matsize); + SymmetricMatrix > Kc(matsize,matsize); if( matsize == 1 ) { assert(size == 1); @@ -1015,9 +1063,9 @@ namespace INMOST /// @param r Array of diagonal elements. /// @param size Size of the matrix. /// @return Matrix with diagonal defined by array, other elements are zero. - static SymmetricMatrix FromDiagonal(const Var * r, enumerator size) + static SymmetricMatrix > FromDiagonal(const Var * r, enumerator size) { - SymmetricMatrix ret(size); + SymmetricMatrix > ret(size); ret.Zero(); for(enumerator k = 0; k < size; ++k) ret(k,k) = r[k]; return ret; @@ -1026,9 +1074,9 @@ namespace INMOST /// @param r Array of diagonal elements. /// @param size Size of the matrix. /// @return Matrix with diagonal defined by inverse of array elements. - static SymmetricMatrix FromDiagonalInverse(const Var * r, enumerator size) + static SymmetricMatrix > FromDiagonalInverse(const Var * r, enumerator size) { - SymmetricMatrix ret(size); + SymmetricMatrix > ret(size); ret.Zero(); for(enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k]; return ret; @@ -1038,9 +1086,9 @@ namespace INMOST /// @param pn Number of rows and columns in the matrix. /// @param c Value to put onto diagonal. /// @return Returns a unit matrix. - static SymmetricMatrix Unit(enumerator pn, const Var & c = 1.0) + static SymmetricMatrix > Unit(enumerator pn, const Var & c = 1.0) { - SymmetricMatrix ret(pn,0.0); + SymmetricMatrix > ret(pn,0.0); for(enumerator i = 0; i < pn; ++i) ret(i,i) = c; return ret; } @@ -1359,9 +1407,9 @@ namespace INMOST /// @param Perm Array with new positions for rows. /// @param size Size of the array and the resulting matrix. /// @return Permutation matrix. - static Matrix Permutation(const INMOST_DATA_ENUM_TYPE * Perm, enumerator size) + static Matrix > Permutation(const INMOST_DATA_ENUM_TYPE * Perm, enumerator size) { - Matrix Ret(size,size); + Matrix > Ret(size,size); Ret.Zero(); for(enumerator k = 0; k < size; ++k) Ret(k,Perm[k]) = 1; @@ -1388,9 +1436,9 @@ namespace INMOST /// @param size Size of the input array. /// @param matsize Size of the final tensor. /// @return Matrix of the tensor of size matsize by matsize. - static Matrix FromTensor(const Var * K, enumerator size, enumerator matsize = 3) + static Matrix > FromTensor(const Var * K, enumerator size, enumerator matsize = 3) { - Matrix Kc(matsize,matsize); + Matrix > Kc(matsize,matsize); if( matsize == 1 ) { assert(size == 1); @@ -1506,17 +1554,17 @@ namespace INMOST /// @param r Array of elements of the vector. /// @param size Size of the vector. /// @return Vector with contents of the array. - static Matrix FromVector(const Var * r, enumerator size) + static Matrix > FromVector(const Var * r, enumerator size) { - return Matrix(r,size,1); + return Matrix >(r,size,1); } /// Create diagonal matrix from array /// @param r Array of diagonal elements. /// @param size Size of the matrix. /// @return Matrix with diagonal defined by array, other elements are zero. - static Matrix FromDiagonal(const Var * r, enumerator size) + static Matrix > FromDiagonal(const Var * r, enumerator size) { - Matrix ret(size,size); + Matrix > ret(size,size); ret.Zero(); for(enumerator k = 0; k < size; ++k) ret(k,k) = r[k]; return ret; @@ -1525,9 +1573,9 @@ namespace INMOST /// @param r Array of diagonal elements. /// @param size Size of the matrix. /// @return Matrix with diagonal defined by inverse of array elements. - static Matrix FromDiagonalInverse(const Var * r, enumerator size) + static Matrix > FromDiagonalInverse(const Var * r, enumerator size) { - Matrix ret(size,size); + Matrix > ret(size,size); ret.Zero(); for(enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k]; return ret; @@ -1538,12 +1586,12 @@ namespace INMOST /// and vector. For a x b equivalent is CrossProduct(a)*b. /// @param vec Array of elements representing a vector. /// @return A matrix representing cross product. - static Matrix CrossProduct(const Var vec[3]) + static Matrix > CrossProduct(const Var vec[3]) { // | 0 -z y | // | z 0 -x | // | -y x 0 | - Matrix ret(3,3); + Matrix > ret(3,3); ret(0,0) = 0.0; ret(0,1) = -vec[2]; //-z ret(0,2) = vec[1]; //y @@ -1560,9 +1608,9 @@ namespace INMOST /// @param pn Number of rows and columns in the matrix. /// @param c Value to put onto diagonal. /// @return Returns a unit matrix. - static Matrix Unit(enumerator pn, const Var & c = 1.0) + static Matrix > Unit(enumerator pn, const Var & c = 1.0) { - Matrix ret(pn,pn); + Matrix > ret(pn,pn); ret.Zero(); for(enumerator i = 0; i < pn; ++i) ret(i,i) = c; return ret; @@ -1572,18 +1620,18 @@ namespace INMOST /// @param pn Number of columns. /// @param c Value to fill the matrix. /// @return Returns a matrix with 1 row. - static Matrix Row(enumerator pn, const Var & c = 1.0) + static Matrix > Row(enumerator pn, const Var & c = 1.0) { - return Matrix(1,pn,c); + return Matrix >(1,pn,c); } /// Matix with 1 column, Create a matrix of size pn by 1 and /// fills it with c. /// @param pn Number of rows. /// @param c Value to fill the matrix. /// @return Returns a matrix with 1 column. - static Matrix Col(enumerator pn, const Var & c = 1.0) + static Matrix > Col(enumerator pn, const Var & c = 1.0) { - return Matrix(pn, 1, c); + return Matrix >(pn, 1, c); } /// Concatenate B matrix as columns of current matrix. /// Assumes that number of rows of current matrix is @@ -1591,10 +1639,10 @@ namespace INMOST /// @param B Matrix to be concatenated to current matrix. /// @return Result of concatenation of current matrix and parameter. /// @see Matrix::ConcatRows - Matrix ConcatCols(const Matrix & B) const + Matrix > ConcatCols(const Matrix & B) const { assert(Rows() == B.Rows()); - Matrix ret(Rows(),Cols()+B.Cols()); + Matrix > ret(Rows(),Cols()+B.Cols()); const Matrix & A = *this; for(enumerator i = 0; i < Rows(); ++i) { @@ -1611,10 +1659,10 @@ namespace INMOST /// @param B Matrix to be concatenated to current matrix. /// @return Result of concatenation of current matrix and parameter. /// @see Matrix::ConcatCols - Matrix ConcatRows(const Matrix & B) const + Matrix > ConcatRows(const Matrix & B) const { assert(Cols() == B.Cols()); - Matrix ret(Rows()+B.Rows(),Cols()); + Matrix > ret(Rows()+B.Rows(),Cols()); const Matrix & A = *this; for(enumerator i = 0; i < Rows(); ++i) { @@ -1640,13 +1688,13 @@ namespace INMOST /// @return A unitary n by n matrix V used to diagonalize array of /// initial matrices. Current matrix is replaced by concatination of /// V^T*A_i*V, a collection of diagonalized matrices. - Matrix JointDiagonalization(INMOST_DATA_REAL_TYPE threshold = 1.0e-7) + Matrix > JointDiagonalization(INMOST_DATA_REAL_TYPE threshold = 1.0e-7) { enumerator N = Rows(); enumerator M = Cols() / Rows(); - Matrix V = Matrix::Unit(m); - Matrix R(2,M); - Matrix G(2,2); + Matrix > V(Matrix >::Unit(m)); + Matrix > R(2,M); + Matrix > G(2,2); Matrix & A = *this; Var ton, toff, theta, c, s, Ap, Aq, Vp, Vq; bool repeat; @@ -1812,9 +1860,9 @@ namespace INMOST /// not affect elements of the submatrix or original matrix /// used to create submatrix. /// @return Matrix with same entries as submatrix. - ::INMOST::Matrix MakeMatrix() + ::INMOST::Matrix > MakeMatrix() { - ::INMOST::Matrix ret(Rows(),Cols()); + ::INMOST::Matrix > ret(Rows(),Cols()); for(enumerator i = 0; i < Rows(); ++i) for(enumerator j = 0; j < Cols(); ++j) ret(i,j) = (*this)(i,j); @@ -1888,9 +1936,9 @@ namespace INMOST /// not affect elements of the submatrix or original matrix /// used to create submatrix. /// @return Matrix with same entries as submatrix. - ::INMOST::Matrix MakeMatrix() + ::INMOST::Matrix > MakeMatrix() { - ::INMOST::Matrix ret(Rows(),Cols()); + ::INMOST::Matrix > ret(Rows(),Cols()); for(enumerator i = 0; i < Rows(); ++i) for(enumerator j = 0; j < Cols(); ++j) ret(i,j) = (*this)(i,j); @@ -1908,10 +1956,10 @@ namespace INMOST }; template - Matrix + Matrix > AbstractMatrix::Transpose() const { - Matrix ret(Cols(),Rows()); + Matrix > ret(Cols(),Rows()); for(enumerator i = 0; i < Rows(); ++i) for(enumerator j = 0; j < Cols(); ++j) ret(j,i) = (*this)(i,j); @@ -1922,16 +1970,16 @@ namespace INMOST void AbstractMatrix::Swap(AbstractMatrix & b) { - Matrix tmp = b; + Matrix > tmp = b; b = (*this); (*this) = tmp; } template - Matrix + Matrix > AbstractMatrix::operator-() const { - Matrix ret(Rows(),Cols()); + Matrix > ret(Rows(),Cols()); for(enumerator i = 0; i < Rows(); ++i) for(enumerator j = 0; j < Cols(); ++j) ret(i,j) = -(*this)(i,j); @@ -1940,12 +1988,68 @@ namespace INMOST template template - Matrix::type> + Matrix::type, pool_array_t::type> > + AbstractMatrix::CrossProduct(const AbstractMatrix & B) const + { + const AbstractMatrix & A = *this; + assert(A.Rows() == 3); + assert(A.Cols() == 1); + assert(B.Rows() == 3); + Matrix::type, pool_array_t::type> > ret(3,B.Cols()); //check RVO + for(unsigned k = 0; k < B.Cols(); ++k) + { + ret(0,k) = (A(1,0)*B(2,k) - A(2,0)*B(1,k)); + ret(1,k) = (A(2,0)*B(0,k) - A(0,0)*B(2,k)); + ret(2,k) = (A(0,0)*B(1,k) - A(1,0)*B(0,k)); + } + return ret; + } + + template + template + Matrix::type, pool_array_t::type> > + AbstractMatrix::Transform(const AbstractMatrix & B) const + { + const AbstractMatrix & A = *this; + assert(A.Rows() == 3); + assert(A.Cols() == 1); + assert(B.Rows() == 3); + assert(B.Cols() == 1); + typedef typename Promote::type var_t; + Matrix > Q(3,3); + var_t x,y,z,w,n,l; + + x = -(B(1,0)*A(2,0) - B(2,0)*A(1,0)); + y = -(B(2,0)*A(0,0) - B(0,0)*A(2,0)); + z = -(B(0,0)*A(1,0) - B(1,0)*A(0,0)); + w = A.FrobeniusNorm()*B.FrobeniusNorm() + A.DotProduct(B); + + l = B.FrobeniusNorm()/A.FrobeniusNorm(); + n = 2*l/(x*x+y*y+z*z+w*w); + + Q(0,0) = l - (y*y + z*z)*n; + Q(0,1) = (x*y - z*w)*n; + Q(0,2) = (x*z + y*w)*n; + + Q(1,0) = (x*y + z*w)*n; + Q(1,1) = l - (x*x + z*z)*n; + Q(1,2) = (y*z - x*w)*n; + + Q(2,0) = (x*z - y*w)*n; + Q(2,1) = (y*z + x*w)*n; + Q(2,2) = l - (x*x + y*y)*n; + + return Q; + } + + template + template + Matrix::type, pool_array_t::type> > AbstractMatrix::operator-(const AbstractMatrix & other) const { assert(Rows() == other.Rows()); assert(Cols() == other.Cols()); - Matrix::type> ret(Rows(),Cols()); //check RVO + Matrix::type, pool_array_t::type> > ret(Rows(),Cols()); //check RVO for(enumerator i = 0; i < Rows(); ++i) for(enumerator j = 0; j < Cols(); ++j) ret(i,j) = (*this)(i,j)-other(i,j); @@ -1967,12 +2071,12 @@ namespace INMOST template template - Matrix::type> + Matrix::type, pool_array_t::type> > AbstractMatrix::operator+(const AbstractMatrix & other) const { assert(Rows() == other.Rows()); assert(Cols() == other.Cols()); - Matrix::type> ret(Rows(),Cols()); //check RVO + Matrix::type, pool_array_t::type> > ret(Rows(),Cols()); //check RVO for(enumerator i = 0; i < Rows(); ++i) for(enumerator j = 0; j < Cols(); ++j) ret(i,j) = (*this)(i,j)+other(i,j); @@ -1995,11 +2099,11 @@ namespace INMOST template template - Matrix::type> + Matrix::type, pool_array_t::type> > AbstractMatrix::operator*(const AbstractMatrix & other) const { assert(Cols() == other.Rows()); - Matrix::type> ret(Rows(),other.Cols()); //check RVO + Matrix::type, pool_array_t::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 @@ -2061,10 +2165,10 @@ namespace INMOST template template - Matrix::type> + Matrix::type, pool_array_t::type> > AbstractMatrix::operator*(typeB coef) const { - Matrix::type> ret(Rows(),Cols()); //check RVO + Matrix::type, pool_array_t::type> > ret(Rows(),Cols()); //check RVO for(enumerator i = 0; i < Rows(); ++i) for(enumerator j = 0; j < Cols(); ++j) ret(i,j) = (*this)(i,j)*coef; @@ -2084,10 +2188,10 @@ namespace INMOST template template - Matrix::type> + Matrix::type, pool_array_t::type> > AbstractMatrix::operator/(typeB coef) const { - Matrix::type> ret(Rows(),Cols()); + Matrix::type, pool_array_t::type> > ret(Rows(),Cols()); for(enumerator i = 0; i < Rows(); ++i) for(enumerator j = 0; j < Cols(); ++j) ret(i,j) = (*this)(i,j)/coef; @@ -2107,10 +2211,10 @@ namespace INMOST template template - Matrix::type> + Matrix::type, pool_array_t::type> > AbstractMatrix::Kronecker(const AbstractMatrix & other) const { - Matrix::type> ret(Rows()*other.Rows(),Cols()*other.Cols()); + Matrix::type, pool_array_t::type> > ret(Rows()*other.Rows(),Cols()*other.Cols()); for(enumerator i = 0; i < Rows(); ++i) //loop rows { for(enumerator j = 0; j < other.Rows(); ++j) //loop columns @@ -2128,35 +2232,36 @@ namespace INMOST } template - std::pair,bool> - AbstractMatrix::Invert(bool print_fail) const + Matrix > + AbstractMatrix::Invert(int * ierr) const { - return Solve(Matrix::Unit(Rows()),print_fail); + Matrix > ret(Cols(),Rows()); + ret = Solve(Matrix >::Unit(Rows()),ierr); + return ret; } template - std::pair,bool> - AbstractMatrix::CholeskyInvert(bool print_fail) const + Matrix > + AbstractMatrix::CholeskyInvert(int * ierr) const { - return CholeskySolve(Matrix::Unit(Rows()),print_fail); + Matrix > ret(Rows(),Rows()); + ret = CholeskySolve(Matrix >::Unit(Rows()),ierr); + return ret; } template template - std::pair::type>,bool> - AbstractMatrix::CholeskySolve(const AbstractMatrix & B, bool print_fail) const + Matrix::type, pool_array_t::type> > + AbstractMatrix::CholeskySolve(const AbstractMatrix & B, int * ierr) const { - - const AbstractMatrix & A = *this; assert(A.Rows() == A.Cols()); assert(A.Rows() == B.Rows()); enumerator n = A.Rows(); enumerator l = B.Cols(); - std::pair::type>,bool> - ret = std::make_pair(Matrix::type>(B),true); - SymmetricMatrix L = A; + Matrix::type, pool_array_t::type> > ret(B); + SymmetricMatrix > L(A); //SAXPY /* @@ -2190,8 +2295,12 @@ namespace INMOST { if( L(k,k) < 0.0 ) { - ret.second = false; - if( print_fail ) std::cout << "Negative diagonal pivot " << get_value(L(k,k)) << std::endl; + if( ierr ) + { + if( *ierr == -1 ) std::cout << "Negative diagonal pivot " << get_value(L(k,k)) << " row " << k << std::endl; + *ierr = k+1; + } + else throw MatrixCholeskySolveFail; return ret; } @@ -2199,8 +2308,12 @@ namespace INMOST if( fabs(L(k,k)) < 1.0e-24 ) { - ret.second = false; - if( print_fail ) std::cout << "Diagonal pivot is too small " << get_value(L(k,k)) << std::endl; + if( ierr ) + { + if( *ierr == -1 ) std::cout << "Diagonal pivot is too small " << get_value(L(k,k)) << " row " << k << std::endl; + *ierr = k+1; + } + else throw MatrixCholeskySolveFail; return ret; } @@ -2214,7 +2327,7 @@ namespace INMOST } } // LY=B - Matrix::type> & Y = ret.first; + Matrix::type, pool_array_t::type> > & Y = ret; for(enumerator i = 0; i < n; ++i) { for(enumerator k = 0; k < l; ++k) @@ -2225,7 +2338,7 @@ namespace INMOST } } // L^TX = Y - Matrix::type> & X = ret.first; + Matrix::type, pool_array_t::type> > & X = ret; for(enumerator it = n; it > 0; --it) { enumerator i = it-1; @@ -2239,13 +2352,14 @@ namespace INMOST X(i,k) /= L(i,i); } } + if( ierr ) *ierr = 0; return ret; } template template - std::pair::type>,bool> - AbstractMatrix::Solve(const AbstractMatrix & B, bool print_fail) const + Matrix::type, pool_array_t::type> > + AbstractMatrix::Solve(const AbstractMatrix & B, int * ierr) const { // for A^T B assert(Rows() == B.Rows()); @@ -2258,15 +2372,18 @@ namespace INMOST Matrix AtB = B; //m by l matrix Matrix AtA = (*this); //m by m matrix */ - Matrix At = this->Transpose(); //m by n matrix - Matrix::type> AtB = At*B; //m by l matrix - Matrix AtA = At*(*this); //m by m matrix - enumerator l = AtB.Cols(); - //enumerator n = Rows(); + enumerator l = B.Cols(); enumerator m = Cols(); - enumerator * order = new enumerator [m]; - std::pair::type>,bool> - ret = std::make_pair(Matrix::type>(m,l),true); + Matrix::type, pool_array_t::type> > ret(m,l); + Matrix > At = this->Transpose(); //m by n matrix + Matrix::type, pool_array_t::type> > AtB = At*B; //m by l matrix + Matrix > AtA = At*(*this); //m by m matrix + assert(l == AtB.Cols()); + //enumerator l = AtB.Cols(); + //enumerator n = Rows(); + + dynarray order(m); + //Var temp; INMOST_DATA_REAL_TYPE max,v; //typeB tempb; @@ -2344,9 +2461,12 @@ namespace INMOST if( ok ) AtA(i,i) = AtA(i,i) < 0.0 ? - 1.0e-12 : 1.0e-12; else { - if( print_fail ) std::cout << "Failed to invert matrix" << std::endl; - ret.second = false; - delete [] order; + if( ierr ) + { + if( *ierr == -1 ) std::cout << "Failed to invert matrix" << std::endl; + *ierr = i+1; + } + else throw MatrixSolveFail; return ret; } } @@ -2381,21 +2501,22 @@ namespace INMOST AtB(i,k) -= AtB(j,k) * AtA(i,j); } for(enumerator i = 0; i < m; i++) - ret.first(order[i],k) = AtB(i,k); + ret(order[i],k) = AtB(i,k); } - delete [] order; + //delete [] order; + if( ierr ) *ierr = 0; return ret; } template - Matrix + Matrix > AbstractMatrix::ExtractSubMatrix(enumerator ibeg, enumerator iend, enumerator jbeg, enumerator jend) const { assert(ibeg < Rows()); assert(iend < Rows()); assert(jbeg < Cols()); assert(jend < Cols()); - Matrix ret(iend-ibeg,jend-jbeg); + Matrix > ret(iend-ibeg,jend-jbeg); for(enumerator i = ibeg; i < iend; ++i) { for(enumerator j = jbeg; j < jend; ++j) @@ -2405,25 +2526,33 @@ namespace INMOST } template - Matrix + Matrix > AbstractMatrix::Repack(enumerator rows, enumerator cols) const { assert(Cols()*Rows()==rows*cols); - Matrix ret(*this); + Matrix > ret(*this); ret.Rows() = rows; ret.Cols() = cols; return ret; } template - std::pair,bool> - AbstractMatrix::PseudoInvert(INMOST_DATA_REAL_TYPE tol, bool print_fail) const + Matrix > + AbstractMatrix::PseudoInvert(INMOST_DATA_REAL_TYPE tol, int * ierr) const { - std::pair,bool> ret = std::make_pair(Matrix(Cols(),Rows()),true); - Matrix U,S,V; - ret.second = SVD(U,S,V); - if( print_fail && !ret.second ) - std::cout << "Failed to compute Moore-Penrose inverse of the matrix" << std::endl; + Matrix > ret(Cols(),Rows()); + Matrix > U,S,V; + bool success = SVD(U,S,V); + if( !success ) + { + if( ierr ) + { + if( *ierr == -1 ) std::cout << "Failed to compute Moore-Penrose inverse of the matrix" << std::endl; + *ierr = 1; + return ret; + } + else throw MatrixPseudoSolveFail; + } for(INMOST_DATA_ENUM_TYPE k = 0; k < S.Cols(); ++k) { if( S(k,k) > tol ) @@ -2431,21 +2560,28 @@ namespace INMOST else S(k,k) = 0.0; } - ret.first = V*S*U.Transpose(); + ret = V*S*U.Transpose(); + if( ierr ) *ierr = 0; return ret; } template template - std::pair::type>,bool> - AbstractMatrix::PseudoSolve(const AbstractMatrix & B, INMOST_DATA_REAL_TYPE tol, bool print_fail) const + Matrix::type, pool_array_t::type> > + AbstractMatrix::PseudoSolve(const AbstractMatrix & B, INMOST_DATA_REAL_TYPE tol, int * ierr) const { - std::pair::type>,bool> - ret = std::make_pair(Matrix::type>(Cols(),B.Cols()),true); - Matrix U,S,V; - ret.second = SVD(U,S,V); - if( print_fail && !ret.second ) - std::cout << "Failed to compute Moore-Penrose inverse of the matrix" << std::endl; + Matrix::type, pool_array_t::type> > ret(Cols(),B.Cols()); + Matrix > U,S,V; + bool success = SVD(U,S,V); + if( !success ) + { + if( ierr ) + { + if( *ierr == -1 ) std::cout << "Failed to compute Moore-Penrose inverse of the matrix" << std::endl; + *ierr = 1; + } + else throw MatrixPseudoSolveFail; + } for(int k = 0; k < S.Cols(); ++k) { if( S(k,k) > tol ) @@ -2453,7 +2589,8 @@ namespace INMOST else S(k,k) = 0.0; } - ret.first = V*S*U.Transpose()*B; + ret = V*S*U.Transpose()*B; + if( ierr ) *ierr = 0; return ret; } @@ -2805,14 +2942,30 @@ namespace INMOST typedef Matrix iMatrix; /// shortcut for matrix of real values. typedef Matrix rMatrix; - /// shortcut for matrix of integer values. + /// shortcut for matrix of integer values in stack-preallocated and dynamically reallocated array. typedef Matrix > idMatrix; - /// shortcut for matrix of real values. + /// shortcut for matrix of real values in stack-preallocated and dynamically reallocated array. typedef Matrix > rdMatrix; + /// shortcut for matrix of integer values in pool-allocated array (beaware of deallocation order issue). + typedef Matrix > ipMatrix; + /// shortcut for matrix of real values in pool-allocated array (beaware of deallocation order issue). + typedef Matrix > rpMatrix; /// shortcut for matrix of integer values in existing array. typedef Matrix > iaMatrix; /// shortcut for matrix of real values in existing array. typedef Matrix > raMatrix; + /// shortcut for symmetric matrix of integer values. + typedef SymmetricMatrix iSymmetricMatrix; + /// shortcut for symmetric matrix of real values. + typedef SymmetricMatrix rSymmetricMatrix; + /// shortcut for symmetric matrix of integer values in stack-preallocated and dynamically reallocated array. + typedef SymmetricMatrix > idSymmetricMatrix; + /// shortcut for symmetric matrix of real values in stack-preallocated and dynamically reallocated array. + typedef SymmetricMatrix > rdSymmetricMatrix; + /// shortcut for symmetric matrix of integer values in pool-allocated array (beaware of deallocation order issue). + typedef SymmetricMatrix > ipSymmetricMatrix; + /// shortcut for symmetric matrix of real values in pool-allocated array (beaware of deallocation order issue). + typedef SymmetricMatrix > rpSymmetricMatrix; /// shortcut for matrix of integer values in existing array. typedef SymmetricMatrix > iaSymmetricMatrix; /// shortcut for matrix of real values in existing array. @@ -2829,12 +2982,18 @@ namespace INMOST typedef Matrix vMatrix; //< shortcut for matrix of variables with first and second order derivatives. typedef Matrix hMatrix; - /// shortcut for matrix of variables with single unit entry of first order derivative. + /// shortcut for matrix of variables with single unit entry of first order derivative in stack-preallocated and dynamically reallocated array. typedef Matrix > udMatrix; - /// shortcut for matrix of variables with first order derivatives. + /// shortcut for matrix of variables with first order derivatives in stack-preallocated and dynamically reallocated array. typedef Matrix > vdMatrix; - //< shortcut for matrix of variables with first and second order derivatives. + //< shortcut for matrix of variables with first and second order derivatives in stack-preallocated and dynamically reallocated array. typedef Matrix > hdMatrix; + /// shortcut for matrix of variables with single unit entry of first order derivative in pool-allocated array (beaware of deallocation order issue). + typedef Matrix > upMatrix; + /// shortcut for matrix of variables with first order derivatives in pool-allocated array (beaware of deallocation order issue). + typedef Matrix > vpMatrix; + //< shortcut for matrix of variables with first and second order derivatives in pool-allocated array (beaware of deallocation order issue). + typedef Matrix > hpMatrix; /// shortcut for matrix of unknowns in existing array. typedef Matrix > uaMatrix; /// shortcut for matrix of variables in existing array. @@ -2847,6 +3006,19 @@ namespace INMOST typedef SymmetricMatrix > vaSymmetricMatrix; /// shortcut for matrix of variables in existing array. typedef SymmetricMatrix > haSymmetricMatrix; + /// shortcut for symmetric matrix of variables with single unit entry of first order derivative in stack-preallocated and dynamically reallocated array. + typedef SymmetricMatrix > udSymmetricMatrix; + /// shortcut for symmetric matrix of variables with first order derivatives in stack-preallocated and dynamically reallocated array. + typedef SymmetricMatrix > vdSymmetricMatrix; + //< shortcut for symmetric matrix of variables with first and second order derivatives in stack-preallocated and dynamically reallocated array. + typedef SymmetricMatrix > hdSymmetricMatrix; + /// shortcut for symmetric matrix of variables with single unit entry of first order derivative in pool-allocated array (beaware of deallocation order issue). + typedef SymmetricMatrix > upSymmetricMatrix; + /// shortcut for symmetric matrix of variables with first order derivatives in pool-allocated array (beaware of deallocation order issue). + typedef SymmetricMatrix > vpSymmetricMatrix; + //< shortcut for symmetric matrix of variables with first and second order derivatives in pool-allocated array (beaware of deallocation order issue). + typedef SymmetricMatrix > hpSymmetricMatrix; + static uaMatrix uaMatrixMake(unknown * p, uaMatrix::enumerator n, uaMatrix::enumerator m) {return uaMatrix(shell(p,n*m),n,m);} @@ -2862,7 +3034,8 @@ namespace INMOST /// @param other Matrix to be multiplied. /// @return Matrix, each entry multiplied by a constant. template -INMOST::Matrix::type> operator *(INMOST_DATA_REAL_TYPE coef, const INMOST::AbstractMatrix & other) +INMOST::Matrix::type, INMOST::pool_array_t::type> > +operator *(INMOST_DATA_REAL_TYPE coef, const INMOST::AbstractMatrix & other) {return other*coef;} #if defined(USE_AUTODIFF) /// Multiplication of matrix by a unknown from left. @@ -2871,7 +3044,8 @@ INMOST::Matrix::type> oper /// @param other Matrix to be multiplied. /// @return Matrix, each entry multiplied by a variable. template -INMOST::Matrix::type> operator *(const INMOST::unknown & coef, const INMOST::AbstractMatrix & other) +INMOST::Matrix::type, INMOST::pool_array_t::type> > +operator *(const INMOST::unknown & coef, const INMOST::AbstractMatrix & other) {return other*coef;} /// Multiplication of matrix by a variable from left. /// Takes account for derivatives of variable. @@ -2879,7 +3053,8 @@ INMOST::Matrix::type> operator * /// @param other Matrix to be multiplied. /// @return Matrix, each entry multiplied by a variable. template -INMOST::Matrix::type> operator *(const INMOST::variable & coef, const INMOST::AbstractMatrix & other) +INMOST::Matrix::type, INMOST::pool_array_t::type> > +operator *(const INMOST::variable & coef, const INMOST::AbstractMatrix & other) {return other*coef;} /// Multiplication of matrix by a variable with first and /// second order derivatives from left. @@ -2888,7 +3063,8 @@ INMOST::Matrix::type> operator /// @param other Matrix to be multiplied. /// @return Matrix, each entry multiplied by a variable. template -INMOST::Matrix::type> operator *(const INMOST::hessian_variable & coef, const INMOST::AbstractMatrix & other) +INMOST::Matrix::type, INMOST::pool_array_t::type> > +operator *(const INMOST::hessian_variable & coef, const INMOST::AbstractMatrix & other) {return other*coef;} #endif diff --git a/Source/Headers/inmost_expression.h b/Source/Headers/inmost_expression.h index 7a51041..3dacdc2 100644 --- a/Source/Headers/inmost_expression.h +++ b/Source/Headers/inmost_expression.h @@ -21,6 +21,7 @@ //10.1 user should be able to provide RowMerger when Automatizator is not compiled //10.2 Automatizator may provide internal structure for RowMerger + #ifdef _MSC_VER #pragma warning(disable : 4503) #endif @@ -118,7 +119,7 @@ namespace INMOST } bool check_infs() const { - return std::isinf(value); + return __isinf__(value); } }; @@ -292,9 +293,9 @@ namespace INMOST } bool check_infs() const { - if( std::isinf(value) ) return true; + if( __isinf__(value) ) return true; for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it) - if( std::isinf(it->second) ) return true; + if( __isinf__(it->second) ) return true; return false; } /// Write variable into array of entries. @@ -571,11 +572,11 @@ namespace INMOST } bool check_infs() const { - if( std::isinf(value)) return true; + if( __isinf__(value)) return true; for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it) - if( std::isinf(it->second) ) return true; + if( __isinf__(it->second) ) return true; for(Sparse::HessianRow::const_iterator it = hessian_entries.Begin(); it != hessian_entries.End(); ++it) - if( std::isinf(it->second) ) return true; + if( __isinf__(it->second) ) return true; return false; } friend class hessian_multivar_expression_reference; @@ -766,9 +767,9 @@ namespace INMOST } bool check_infs() const { - if( std::isinf(value) ) return true; + if( __isinf__(value) ) return true; for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it) - if( std::isinf(it->second) ) return true; + if( __isinf__(it->second) ) return true; return false; } }; @@ -979,11 +980,11 @@ namespace INMOST } bool check_infs() const { - if( std::isinf(value) ) return true; + if( __isinf__(value) ) return true; for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it) - if( std::isinf(it->second) ) return true; + if( __isinf__(it->second) ) return true; for(Sparse::HessianRow::iterator it = hentries->Begin(); it != hentries->End(); ++it) - if( std::isinf(it->second) ) return true; + if( __isinf__(it->second) ) return true; return false; } }; @@ -2146,7 +2147,7 @@ __INLINE bool check_nans(INMOST_DATA_REAL_TYPE val) {return val != val;} __INLINE bool check_nans(INMOST::var_expression const & e) {return e.check_nans();} __INLINE bool check_nans(INMOST::multivar_expression const & e) {return e.check_nans();} __INLINE bool check_nans(INMOST::multivar_expression_reference const & e) {return e.check_nans();} -__INLINE bool check_infs(INMOST_DATA_REAL_TYPE val) {return std::isinf(val);} +__INLINE bool check_infs(INMOST_DATA_REAL_TYPE val) {return __isinf__(val);} __INLINE bool check_infs(INMOST::var_expression const & e) {return e.check_infs();} __INLINE bool check_infs(INMOST::multivar_expression const & e) {return e.check_infs();} __INLINE bool check_infs(INMOST::multivar_expression_reference const & e) {return e.check_infs();} @@ -2284,7 +2285,7 @@ __INLINE INMOST_DATA_REAL_TYPE get_table(INMOST_DATA_RE #else //USE_AUTODIFF __INLINE bool check_nans(INMOST_DATA_REAL_TYPE val) {return val != val;} -__INLINE bool check_infs(INMOST_DATA_REAL_TYPE val) {return std::isinf(val);} +__INLINE bool check_infs(INMOST_DATA_REAL_TYPE val) {return __isinf__(val);} __INLINE bool check_nans_infs(INMOST_DATA_REAL_TYPE val) {return check_nans(val) || check_infs(val);} __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = Val;} __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;} diff --git a/Source/IO/mesh_vtk_file.cpp b/Source/IO/mesh_vtk_file.cpp index 4b4f8aa..8437e04 100644 --- a/Source/IO/mesh_vtk_file.cpp +++ b/Source/IO/mesh_vtk_file.cpp @@ -27,10 +27,10 @@ #define R_QUIT 100 -static int __isnan__(double x) { return x != x; } +//static int __isnan__(double x) { return x != x; } //static int isinf(double x) { return !isnan(x) && isnan(x - x); } -static int __isinf__(double x) { return fabs(x) > DBL_MAX; } -static int __isbad(double x) { return __isnan__(x) || __isinf__(x); } +//static int __isinf__(double x) { return fabs(x) > DBL_MAX; } +//static int __isbad(double x) { return __isnan__(x) || __isinf__(x); } template void ReadCoords(FILE * f,INMOST_DATA_REAL_TYPE c[3]) diff --git a/Source/IO/mesh_vtu_file.cpp b/Source/IO/mesh_vtu_file.cpp index e7492ba..fec62c0 100644 --- a/Source/IO/mesh_vtu_file.cpp +++ b/Source/IO/mesh_vtu_file.cpp @@ -8,10 +8,10 @@ #if defined(USE_MESH) -static int __isnan__(double x) { return x != x; } +//static int __isnan__(double x) { return x != x; } //static int isinf(double x) { return !isnan(x) && isnan(x - x); } -static int __isinf__(double x) { return fabs(x) > DBL_MAX; } -static int __isbad(double x) { return __isnan__(x) || __isinf__(x); } +//static int __isinf__(double x) { return fabs(x) > DBL_MAX; } +//static int __isbad(double x) { return __isnan__(x) || __isinf__(x); } diff --git a/Source/Misc/CMakeLists.txt b/Source/Misc/CMakeLists.txt index 00f9beb..3d913e7 100644 --- a/Source/Misc/CMakeLists.txt +++ b/Source/Misc/CMakeLists.txt @@ -4,6 +4,7 @@ set(SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/xml.cpp${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp ${CMAKE_CURRENT_SOURCE_DIR}/base64.cpp +${CMAKE_CURRENT_SOURCE_DIR}/pool.cpp PARENT_SCOPE ) diff --git a/Source/Misc/pool.cpp b/Source/Misc/pool.cpp new file mode 100644 index 0000000..09ef7d7 --- /dev/null +++ b/Source/Misc/pool.cpp @@ -0,0 +1,15 @@ +#ifdef _MSC_VER //kill some warnings +#define _CRT_SECURE_NO_WARNINGS +#endif +#include "inmost.h" + +namespace INMOST +{ + thread_private _pool; + + + memory_pool & get_pool() + { + return *_pool; + } +} diff --git a/Source/Solver/solver_inner/SolverInner.h b/Source/Solver/solver_inner/SolverInner.h index 3d36821..ccd5580 100644 --- a/Source/Solver/solver_inner/SolverInner.h +++ b/Source/Solver/solver_inner/SolverInner.h @@ -5,7 +5,7 @@ #include "../Misc/utils.h" #include "solver_prototypes.hpp" #include "solver_bcgsl.hpp" -#define KSOLVER BCGS_solver +#define KSOLVER BCGSL_solver namespace INMOST { diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 08899c8..f065e88 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -1,27 +1,27 @@ -if(USE_MESH) -add_subdirectory(geom_test000) -endif(USE_MESH) - -if(USE_AUTODIFF) -add_subdirectory(autodiff_test000) -add_subdirectory(autodiff_test001) -add_subdirectory(autodiff_test002) -endif(USE_AUTODIFF) - -if(USE_SOLVER) -add_subdirectory(solver_test000) -#add_subdirectory(solver_test001) -add_subdirectory(solver_test002) -add_subdirectory(solver_test003) -endif(USE_SOLVER) - -if(USE_MESH AND USE_MPI) -add_subdirectory(pmesh_test000) -if(USE_PARTITIONER) -add_subdirectory(pmesh_test001) -endif() -endif() - - - -add_subdirectory(xml_reader_test000) +if(USE_MESH) +add_subdirectory(geom_test000) +endif(USE_MESH) + +if(USE_AUTODIFF) +add_subdirectory(autodiff_test000) +add_subdirectory(autodiff_test001) +add_subdirectory(autodiff_test002) +endif(USE_AUTODIFF) + +if(USE_SOLVER) +add_subdirectory(solver_test000) +#add_subdirectory(solver_test001) +add_subdirectory(solver_test002) +add_subdirectory(solver_test003) +endif(USE_SOLVER) + +if(USE_MESH AND USE_MPI) +add_subdirectory(pmesh_test000) +if(USE_PARTITIONER) +add_subdirectory(pmesh_test001) +endif() +endif() + + +add_subdirectory(linalg_test000) +add_subdirectory(xml_reader_test000) diff --git a/Tests/linalg_test000/CMakeLists.txt b/Tests/linalg_test000/CMakeLists.txt new file mode 100644 index 0000000..9e940b8 --- /dev/null +++ b/Tests/linalg_test000/CMakeLists.txt @@ -0,0 +1,37 @@ +project(linalg_test000) +set(SOURCE main.cpp) + +add_executable(linalg_test000 ${SOURCE}) +target_link_libraries(linalg_test000 inmost) + +if(USE_MPI) + message("linking linalg_test000 with MPI") + target_link_libraries(linalg_test000${MPI_LIBRARIES}) + if(MPI_LINK_FLAGS) + set_target_properties(linalg_test000 PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") + endif() +endif(USE_MPI) + + +add_test(NAME linalg_test000_matmul COMMAND$ 0) +add_test(NAME linalg_test000_matdiv COMMAND $1) +add_test(NAME linalg_test000_matadd COMMAND$ 2) +add_test(NAME linalg_test000_matsub COMMAND $3) +add_test(NAME linalg_test000_svd_square COMMAND$ 4) +add_test(NAME linalg_test000_svd_ngtm COMMAND $5) +add_test(NAME linalg_test000_svd_nltm COMMAND$ 6) +add_test(NAME linalg_test000_inv_square COMMAND $7) +add_test(NAME linalg_test000_inv_ngtm COMMAND$ 8) +add_test(NAME linalg_test000_inv_nltm COMMAND $9) +add_test(NAME linalg_test000_cholinv_square COMMAND$ 10) +add_test(NAME linalg_test000_pinv_square COMMAND $11) +add_test(NAME linalg_test000_pinv_ngtm COMMAND$ 12) +add_test(NAME linalg_test000_pinv_nltm COMMAND $13) +add_test(NAME linalg_test000_sol_square COMMAND$ 14) +add_test(NAME linalg_test000_sol_ngtm COMMAND $15) +add_test(NAME linalg_test000_sol_nltm COMMAND$ 16) +add_test(NAME linalg_test000_cholsol_square COMMAND $17) +add_test(NAME linalg_test000_psol_square COMMAND$ 18) +add_test(NAME linalg_test000_psol_ngtm COMMAND $19) +add_test(NAME linalg_test000_psol_nltm COMMAND$ 20) +add_test(NAME linalg_test000_transform COMMAND \$ 21) diff --git a/Tests/linalg_test000/main.cpp b/Tests/linalg_test000/main.cpp new file mode 100644 index 0000000..3d60170 --- /dev/null +++ b/Tests/linalg_test000/main.cpp @@ -0,0 +1,419 @@ +#include "inmost.h" +#include +#include +using namespace INMOST; + +int main(int argc,char ** argv) +{ + int test = 0; + if (argc > 1) test = atoi(argv[1]); + + double err = 1; + + if( test == 0 ) //A*B + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0 + }; + double B[] = + { + 4, 2, + 3, 5, + 9, 8, + 1, 7 + }; + double C[] = + { + 36, 68, + 55, 69, + 103, 90 + }; + err = (raMatrixMake(A,3,4)*raMatrixMake(B,4,2)-raMatrixMake(C,3,2)).FrobeniusNorm(); + } + else if( test == 1 ) // A/B == B^{-1}*A + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0 + }; + double B[] = + { + 4, 2, + 9, 8, + 2, 1 + }; + double C[] = + { + 0.942857142857144, 0.685714285714285, 2.685714285714290, 2.142857142857140, + -0.685714285714285, 0.228571428571429, -2.771428571428570, -2.285714285714290 + }; + err = (raMatrixMake(A,3,4)/raMatrixMake(B,3,2)-raMatrixMake(C,2,4)).FrobeniusNorm(); + } + else if( test == 2 ) // A+B + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0 + }; + double B[] = + { + 2, 6, 9, 1, + 9, 1, 5, 3, + 8, 21, 91, 9 + }; + double C[] = + { + 3, 9, 11, 6, + 12, 9, 7, 4, + 12, 23, 100, 9 + }; + err = (raMatrixMake(A,3,4)+raMatrixMake(B,3,4)-raMatrixMake(C,3,4)).FrobeniusNorm(); + } + else if( test == 3 ) // A-B + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0 + }; + double B[] = + { + 2, 6, 9, 1, + 9, 1, 5, 3, + 8, 21, 91, 9 + }; + double C[] = + { + -1, -3, -7, 4, + -6, 7, -3, -2, + -4, -19, -82, -9 + }; + err = (raMatrixMake(A,3,4)-raMatrixMake(B,3,4)-raMatrixMake(C,3,4)).FrobeniusNorm(); + } + else if( test == 4 ) // svd square + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0, + 2, 1, 10, 9 + }; + rMatrix U, S, V; + raMatrix mA = raMatrixMake(A,4,4); + mA.SVD(U,S,V); + err = (U*S*V.Transpose() - mA).FrobeniusNorm(); + } + else if( test == 5 ) // svd n > m + { + double A[] = + { + 1, 3, 2, + 3, 8, 2, + 4, 2, 9, + 2, 1, 10 + }; + rMatrix U, S, V; + raMatrix mA = raMatrixMake(A,4,3); + mA.SVD(U,S,V); + err = (U*S*V.Transpose() - mA).FrobeniusNorm(); + } + else if( test == 6 ) // svd n < m + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0 + }; + rMatrix U, S, V; + raMatrix mA = raMatrixMake(A,3,4); + mA.SVD(U,S,V); + err = (U*S*V.Transpose() - mA).FrobeniusNorm(); + } + else if( test == 7 ) // inv square + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0, + 2, 1, 10, 9 + }; + raMatrix mA = raMatrixMake(A,4,4); + err = (mA.Invert()*mA-rMatrix::Unit(4)).FrobeniusNorm(); + err+= (mA*mA.Invert()-rMatrix::Unit(4)).FrobeniusNorm(); + } + else if( test == 8 ) // inv n > m + { + double A[] = + { + 1, 3, 2, + 3, 8, 2, + 4, 2, 9, + 2, 1, 10 + }; + raMatrix mA = raMatrixMake(A,4,3); + err = (mA.Invert()*mA-rMatrix::Unit(3)).FrobeniusNorm(); + } + else if( test == 9 ) // inv n < m + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0 + }; + raMatrix mA = raMatrixMake(A,3,4); + err = (mA*mA.Invert()-rMatrix::Unit(3)).FrobeniusNorm(); + } + else if( test == 10 ) // cholinv square + { + double A[] = + { + 10, 3, 2, 5, + 3, 8, 2, 1, + 2, 2, 9, 0, + 5, 1, 0, 9 + }; + raMatrix mA = raMatrixMake(A,4,4); + err = (mA.CholeskyInvert()*mA-rMatrix::Unit(4)).FrobeniusNorm(); + err+= (mA*mA.CholeskyInvert()-rMatrix::Unit(4)).FrobeniusNorm(); + } + else if( test == 11 ) // pinv square + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0, + 2, 1, 10, 9 + }; + raMatrix mA = raMatrixMake(A,4,4); + err = (mA.PseudoInvert()*mA-rMatrix::Unit(4)).FrobeniusNorm(); + err+= (mA*mA.PseudoInvert()-rMatrix::Unit(4)).FrobeniusNorm(); + } + else if( test == 12 ) // pinv n > m + { + double A[] = + { + 1, 3, 2, + 3, 8, 2, + 4, 2, 9, + 2, 1, 10 + }; + raMatrix mA = raMatrixMake(A,4,3); + err = (mA.PseudoInvert()*mA-rMatrix::Unit(3)).FrobeniusNorm(); + } + else if( test == 13 ) // pinv n < m + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0 + }; + raMatrix mA = raMatrixMake(A,3,4); + err = (mA*mA.PseudoInvert()-rMatrix::Unit(3)).FrobeniusNorm(); + } + else if( test == 14 ) // sol square + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0, + 2, 1, 10, 9 + }; + double B[] = + { + 1, 2, + 3, 9, + 4, 1, + 2, 5 + }; + double C[] = + { + 1.000000000000000, -8.333333333333330, + 0.000000000000000, 3.666666666666670, + -0.000000000000001, 3.000000000000000, + 0.000000000000001, -1.333333333333330 + }; + raMatrix mA = raMatrixMake(A,4,4); + raMatrix mB = raMatrixMake(B,4,2); + raMatrix mC = raMatrixMake(C,4,2); + err = (mA.Solve(mB)-mC).FrobeniusNorm(); + } + else if( test == 15 ) // sol n > m + { + double A[] = + { + 1, 3, 2, + 3, 8, 2, + 4, 2, 9, + 2, 1, 10 + }; + double B[] = + { + 1, 2, + 3, 9, + 4, 1, + 2, 5 + }; + raMatrix mA = raMatrixMake(A,4,3); + raMatrix mB = raMatrixMake(B,4,2); + //raMatrix mC = raMatrixMake(C,4,2); + //err = (mA.PseudoSolve(mB)-mC).FrobeniusNorm(); + err = (mA.Transpose()*(mA*mA.Solve(mB)-mB)).FrobeniusNorm(); + } + else if( test == 16 ) // sol n < m + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0 + }; + double B[] = + { + 1, 2, + 3, 9, + 4, 1, + }; + raMatrix mA = raMatrixMake(A,3,4); + raMatrix mB = raMatrixMake(B,3,2); + //raMatrix mC = raMatrixMake(C,4,2); + //err = (mA.PseudoSolve(mB)-mC).FrobeniusNorm(); + err = (mA*mA.Solve(mB)-mB).FrobeniusNorm(); + } + else if( test == 17 ) // cholsol square + { + double A[] = + { + 10, 3, 2, 5, + 3, 8, 2, 1, + 2, 2, 9, 0, + 5, 1, 0, 9 + }; + double B[] = + { + 1, 2, + 3, 9, + 4, 1, + 2, 5 + }; + double C[] = + { + -0.241562583045445, -0.514217379750200, + 0.318628753654002, 1.242625564709010, + 0.427318628753653, -0.050757374435302, + 0.321020462397020, 0.703162370449109 + }; + raMatrix mA = raMatrixMake(A,4,4); + raMatrix mB = raMatrixMake(B,4,2); + raMatrix mC = raMatrixMake(C,4,2); + err = (mA.CholeskySolve(mB)-mC).FrobeniusNorm(); + } + else if( test == 18 ) // psol square + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0, + 2, 1, 10, 9 + }; + double B[] = + { + 1, 2, + 3, 9, + 4, 1, + 2, 5 + }; + double C[] = + { + 1.000000000000000, -8.333333333333330, + 0.000000000000000, 3.666666666666670, + -0.000000000000001, 3.000000000000000, + 0.000000000000001, -1.333333333333330 + }; + raMatrix mA = raMatrixMake(A,4,4); + raMatrix mB = raMatrixMake(B,4,2); + raMatrix mC = raMatrixMake(C,4,2); + err = (mA.PseudoSolve(mB)-mC).FrobeniusNorm(); + //err = (mA*mA.PseudoSolve(mB)-mB).FrobeniusNorm(); + } + else if( test == 19 ) // psol n > m + { + double A[] = + { + 1, 3, 2, + 3, 8, 2, + 4, 2, 9, + 2, 1, 10 + }; + double B[] = + { + 1, 2, + 3, 9, + 4, 1, + 2, 5 + }; + raMatrix mA = raMatrixMake(A,4,3); + raMatrix mB = raMatrixMake(B,4,2); + //raMatrix mC = raMatrixMake(C,4,2); + //err = (mA.PseudoSolve(mB)-mC).FrobeniusNorm(); + err = (mA.Transpose()*(mA*mA.PseudoSolve(mB)-mB)).FrobeniusNorm(); + } + else if( test == 20 ) // psol n < m + { + double A[] = + { + 1, 3, 2, 5, + 3, 8, 2, 1, + 4, 2, 9, 0 + }; + double B[] = + { + 1, 2, + 3, 9, + 4, 1, + }; + raMatrix mA = raMatrixMake(A,3,4); + raMatrix mB = raMatrixMake(B,3,2); + //raMatrix mC = raMatrixMake(C,4,2); + //err = (mA.PseudoSolve(mB)-mC).FrobeniusNorm(); + err = (mA*mA.PseudoSolve(mB)-mB).FrobeniusNorm(); + } + else if( test == 21 ) // transform + { + double A[] = {2,4,9}; + double B[] = {-1,5,-21}; + raMatrix mA = raMatrixMake(A,3,1); + raMatrix mB = raMatrixMake(B,3,1); + rdMatrix Q = mA.Transform(mB); + std::cout << "Q:" << std::endl; + Q.Print(); + std::cout << "Q*mA:" << std::endl; + (Q*mA).Print(); + std::cout << "Q*mB:" << std::endl; + (Q*mB).Print(); + err = (mA.Transform(mB)*mA-mB).FrobeniusNorm(); + } + + if( fabs(err) > 1.0e-10 ) + { + std::cout << "error is " << err << std::endl; + return -1; + } + else std::cout << "no error" << std::endl; + return 0; +}