Memory pool for dense matrix operations and more

Fix compilation issue with std::isinf

Finish memory_pool implementation in container.h

Use of matrices allocated in memory pool for intermediate results in
dense linear algebra operations in inmost_dense.h

Use of memory pool for matrices in inmost_autodiff.h

Change error reporting and return type in Invert, Solve,
CholeskyInvert, CholeskySolve, PseudoInvert, PseudoSolve

Add unit tests for linear algebra on dense matrices.

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<real>(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<real>(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);

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));

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<rMatrix, bool> 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<double>(d);
+	Q = (A*(D - p*p.Transpose())*A.Transpose()).Invert() / static_cast<double>(d);
 	c = A*p;
 	//check
 	if (d == 2)

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

Source/Autodiff/autodiff.cpp

@@ -18,11 +18,21 @@ namespace INMOST
 	template<> Demote<variable>::type                 AbstractEntry::Access<variable>                (const Storage& e, INMOST_DATA_ENUM_TYPE pos) const {return Unknown(e,pos);}
 	template<> Demote<hessian_variable>::type         AbstractEntry::Access<hessian_variable>        (const Storage& e, INMOST_DATA_ENUM_TYPE pos) const {return Unknown(e,pos);}
 	
-	template<> Matrix<Demote<INMOST_DATA_REAL_TYPE>::type>    AbstractEntry::Access<INMOST_DATA_REAL_TYPE>   (const Storage& e) const {return Value(e);}
-	template<> Matrix<Demote<INMOST_DATA_INTEGER_TYPE>::type> AbstractEntry::Access<INMOST_DATA_INTEGER_TYPE>(const Storage& e) const {return Index(e);}
-	template<> Matrix<Demote<unknown>::type>                  AbstractEntry::Access<unknown>                 (const Storage& e) const {return Unknown(e);}
-	template<> Matrix<Demote<variable>::type>                 AbstractEntry::Access<variable>                (const Storage& e) const {return Unknown(e);}
-	template<> Matrix<Demote<hessian_variable>::type>         AbstractEntry::Access<hessian_variable>        (const Storage& e) const {return Unknown(e);}
+	template<>
+	Matrix<Demote<INMOST_DATA_REAL_TYPE>::type, pool_array_t<Demote<INMOST_DATA_REAL_TYPE>::type> >
+	AbstractEntry::Access<INMOST_DATA_REAL_TYPE>   (const Storage& e) const {return Value(e);}
+	template<>
+	Matrix<Demote<INMOST_DATA_INTEGER_TYPE>::type, pool_array_t<Demote<INMOST_DATA_INTEGER_TYPE>::type> >
+	AbstractEntry::Access<INMOST_DATA_INTEGER_TYPE>(const Storage& e) const {return Index(e);}
+	template<>
+	Matrix<Demote<unknown>::type, pool_array_t<Demote<unknown>::type> >
+	AbstractEntry::Access<unknown>(const Storage& e) const {return Unknown(e);}
+	template<>
+	Matrix<Demote<variable>::type, pool_array_t<Demote<variable>::type> >
+	AbstractEntry::Access<variable>(const Storage& e) const {return Unknown(e);}
+	template<>
+	Matrix<Demote<hessian_variable>::type,pool_array_t<Demote<hessian_variable>::type> >
+	AbstractEntry::Access<hessian_variable>(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) )
 		{

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<typename T> Matrix<typename Demote<T>::type> Access(const Storage &e) const;
+		template<typename T>
+		Matrix<typename Demote<T>::type, pool_array_t<typename Demote<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 T> typename Demote<T>::type Access(const Storage &e, INMOST_DATA_ENUM_TYPE pos) const;
+		template<typename T>
+		typename Demote<T>::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<unknown_tag[e].size()); 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(int k = 0; k < (int)unknown_tag[e].size(); ++k) ret(k,0) = Value(e,k); return ret; }
+		rpMatrix Value(const Storage & e) const { rpMatrix ret(MatrixSize(e),1); for(int k = 0; k < (int)unknown_tag[e].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(int k = 0; k < (int)unknown_tag[e].size(); ++k) ret(k,0) = Index(e,k); return ret; }
+		ipMatrix Index(const Storage & e) const { ipMatrix ret(MatrixSize(e),1); for(int k = 0; k < (int)unknown_tag[e].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 VectorEntry::operator [](e);}
+		upMatrix Unknown(const Storage & e) const {return VectorEntry::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(int k = 0; k < (int)unknown_tag[e].size(); ++k) ret(0,0) = Unknown(e,k); return ret; }
+		upMatrix operator [](const Storage & e) const { upMatrix ret(MatrixSize(e),1); for(int k = 0; k < (int)unknown_tag[e].size(); ++k) ret(0,0) = Unknown(e,k); return ret; }
 		/// The intended size of the matrix for this entry.
 		INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const {return (INMOST_DATA_ENUM_TYPE)unknown_tag[e].size();}
 		/// Number of tags in block.
@@ -272,13 +276,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(INMOST_DATA_ENUM_TYPE 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(INMOST_DATA_ENUM_TYPE 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(INMOST_DATA_ENUM_TYPE 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(INMOST_DATA_ENUM_TYPE 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 StatusBlockEntry::operator [](e);}
+		upMatrix Unknown(const Storage & e) const {return StatusBlockEntry::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(INMOST_DATA_ENUM_TYPE 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(INMOST_DATA_ENUM_TYPE 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.
@@ -325,13 +329,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 vector filled with values of unknowns of the block.
-		rMatrix Value(const Storage & e) const;
+		rpMatrix Value(const Storage & e) const;
 		/// Return vector filled with indices of unknowns of the block.
-		iMatrix Index(const Storage & e) const;
+		ipMatrix Index(const Storage & e) const;
 		/// Return vector filled with unknowns of the block with their derivatives.
-		uMatrix Unknown(const Storage & e) const {return MultiEntry::operator [](e);}
+		upMatrix Unknown(const Storage & e) const {return MultiEntry::operator [](e);}
 		/// Return vector filled with unknowns of the block with their derivatives.
-		uMatrix operator [](const Storage & e) const;
+		upMatrix operator [](const Storage & e) const;
 		/// The intended size of the matrix for this entry.
 		INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const;
 		/// Number of tags in block.

Source/Headers/inmost_block_variable.h

@@ -433,11 +433,11 @@ namespace INMOST
 		}
 		/// Get value of variable expression on provided element e.
 		rMatrix Value(const Storage & e) const
-		{return (ArgA.Value(e)/ArgB.Value(e)).first;}
+		{return (ArgA.Value(e)/ArgB.Value(e));}
 		/// 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)/ArgB.Variable(e)).first;}
+		{return (ArgA.Variable(e)/ArgB.Variable(e));}
 		/// 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<abstract_dynamic_block_variable *>(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<abstract_dynamic_block_variable *>(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<abstract_dynamic_block_variable *>(new pseudo_inverse_block_variable(*this));}
 	};

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 <cfloat>
+
+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

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;}

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<typename T>
 void ReadCoords(FILE * f,INMOST_DATA_REAL_TYPE c[3])

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); }
 
 
 

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
 )
 

Source/Misc/pool.cpp

+#ifdef _MSC_VER //kill some warnings
+#define _CRT_SECURE_NO_WARNINGS
+#endif
+#include "inmost.h"
+
+namespace INMOST
+{
+	thread_private<memory_pool> _pool;
+
+
+	memory_pool & get_pool()
+	{
+		return *_pool;
+	}
+}

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 {
 

Tests/CMakeLists.txt

@@ -23,5 +23,5 @@ endif()
 endif()
 
 
-
+add_subdirectory(linalg_test000)
 add_subdirectory(xml_reader_test000)

Tests/linalg_test000/CMakeLists.txt

+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 $<TARGET_FILE:linalg_test000> 0)
+add_test(NAME linalg_test000_matdiv         COMMAND $<TARGET_FILE:linalg_test000> 1)
+add_test(NAME linalg_test000_matadd         COMMAND $<TARGET_FILE:linalg_test000> 2)
+add_test(NAME linalg_test000_matsub         COMMAND $<TARGET_FILE:linalg_test000> 3)
+add_test(NAME linalg_test000_svd_square     COMMAND $<TARGET_FILE:linalg_test000> 4)
+add_test(NAME linalg_test000_svd_ngtm       COMMAND $<TARGET_FILE:linalg_test000> 5)
+add_test(NAME linalg_test000_svd_nltm       COMMAND $<TARGET_FILE:linalg_test000> 6)
+add_test(NAME linalg_test000_inv_square     COMMAND $<TARGET_FILE:linalg_test000> 7)
+add_test(NAME linalg_test000_inv_ngtm       COMMAND $<TARGET_FILE:linalg_test000> 8)
+add_test(NAME linalg_test000_inv_nltm       COMMAND $<TARGET_FILE:linalg_test000> 9)
+add_test(NAME linalg_test000_cholinv_square COMMAND $<TARGET_FILE:linalg_test000> 10)
+add_test(NAME linalg_test000_pinv_square    COMMAND $<TARGET_FILE:linalg_test000> 11)
+add_test(NAME linalg_test000_pinv_ngtm      COMMAND $<TARGET_FILE:linalg_test000> 12)
+add_test(NAME linalg_test000_pinv_nltm      COMMAND $<TARGET_FILE:linalg_test000> 13)
+add_test(NAME linalg_test000_sol_square     COMMAND $<TARGET_FILE:linalg_test000> 14)
+add_test(NAME linalg_test000_sol_ngtm       COMMAND $<TARGET_FILE:linalg_test000> 15)
+add_test(NAME linalg_test000_sol_nltm       COMMAND $<TARGET_FILE:linalg_test000> 16)
+add_test(NAME linalg_test000_cholsol_square COMMAND $<TARGET_FILE:linalg_test000> 17)
+add_test(NAME linalg_test000_psol_square    COMMAND $<TARGET_FILE:linalg_test000> 18)
+add_test(NAME linalg_test000_psol_ngtm      COMMAND $<TARGET_FILE:linalg_test000> 19)
+add_test(NAME linalg_test000_psol_nltm      COMMAND $<TARGET_FILE:linalg_test000> 20)
+add_test(NAME linalg_test000_transform      COMMAND $<TARGET_FILE:linalg_test000> 21)

Tests/linalg_test000/main.cpp

+#include "inmost.h"
+#include <cstdio>
+#include <cmath>
+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();
+	}