Fixes and improvements.

Fixed svd for bcgstab(l);
Transposed matvec;
Interfaces for collective operations for Mesh class.

inmost_mesh.h

@@ -2828,16 +2828,32 @@ namespace INMOST
 		///
 		/// Collective operation.
 		///
-		/// @param input value on current processor
-		/// @return sum over all processors
+		/// @param input Value on current processor
+		/// @return Sum over all processors
 		real                              Integrate          (real input);
 		/// Integrate integer value over all processors.
 		///
 		/// Collective operation.
 		///
-		/// @param input on current processor
-		/// @return sum over all processors
+		/// @param input Value on current processor
+		/// @return Sum over all processors
 		integer                           Integrate          (integer input);
+		/// Integrate an array of real values over all processors.
+		///
+		/// Collective operation.
+		///
+		/// @param input An array of values on current processor.
+		/// @param size Size of the array.
+		/// @return sum over all processors.
+		void                              Integrate          (real * input, integer size);
+		/// Integrate an array of integer values over all processors.
+		///
+		/// Collective operation.
+		///
+		/// @param input An array of values on current processor.
+		/// @param size Size of the array.
+		/// @return Sum over all processors.
+		void                              Integrate          (integer * input, integer size);
 		/// Integrate data corresponding to tag between all processors.
 		/// Elements without the data defined on them or when entry not present will be skipped.
 		///
@@ -2857,6 +2873,13 @@ namespace INMOST
 		integer                           ExclusiveSum       (integer input);
 		real                              AggregateMax       (real input);
 		integer                           AggregateMax       (integer input);
+		void                              AggregateMax       (real * input, integer size);
+		void                              AggregateMax       (integer * input, integer size);
+
+		real                              AggregateMin       (real input);
+		integer                           AggregateMin       (integer input);
+		void                              AggregateMin       (real * input, integer size);
+		void                              AggregateMin       (integer * input, integer size);
 		/// Regather ghosted and shared element sets for data exchange.
 		/// This function will be quite useful if you change statuses of elements
 		/// or modify mesh on your own bypassing internal algorithms.

inmost_solver.h

@@ -415,8 +415,8 @@ namespace INMOST
 			
 			/// Matrix-vector product of the form: y = alpha*A*x + beta * y.
 			/// @param y Input/output vector.
-			/// @see Solver::Vector::Zero
 			void MatVec(INMOST_DATA_REAL_TYPE alpha, Solver::Vector & x, INMOST_DATA_REAL_TYPE beta, Solver::Vector & y) const; //y = alpha*A*x + beta * y
+			void MatVecTranspose(INMOST_DATA_REAL_TYPE alpha, Solver::Vector & x, INMOST_DATA_REAL_TYPE beta, Solver::Vector & y) const; //y = alpha*At*x + beta * y
 
 			/// Clear all data of the matrix.
 			void Clear() {for(Matrix::iterator it = Begin(); it != End(); ++it) it->Clear(); data.clear();}

mesh_parallel.cpp

@@ -350,6 +350,9 @@ namespace INMOST
 #endif//USE_MPI
 		return output;
 	}
+
+
+	
 	
 	Storage::integer Mesh::Integrate(Storage::integer input)
 	{
@@ -361,6 +364,28 @@ namespace INMOST
 #endif//USE_MPI
 		return output;
 	}
+
+	void Mesh::Integrate(Storage::real * input, Storage::integer size)
+	{
+#if defined(USE_MPI)
+		static dynarray<Storage::real,64> temp(size);
+		memcpy(temp.data(),input,sizeof(Storage::real)*size);
+		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_SUM,comm);
+#else//USE_MPI
+		(void) input;
+#endif//USE_MPI
+	}
+
+	void Mesh::Integrate(Storage::integer * input, Storage::integer size)
+	{
+#if defined(USE_MPI)
+		static dynarray<Storage::integer,64> temp(size);
+		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
+		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_SUM,comm);
+#else//USE_MPI
+		(void) input;
+#endif//USE_MPI
+	}
 	
 	Storage::integer Mesh::ExclusiveSum(Storage::integer input)
 	{
@@ -412,6 +437,72 @@ namespace INMOST
 		return output;
 	}
 
+	void Mesh::AggregateMax(Storage::real * input, Storage::integer size)
+	{
+#if defined(USE_MPI)
+		static dynarray<Storage::real,64> temp(size);
+		memcpy(temp.data(),input,sizeof(Storage::real)*size);
+		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_MAX,comm);
+#else//USE_MPI
+		(void) input;
+#endif//USE_MPI
+	}
+
+	void Mesh::AggregateMax(Storage::integer * input, Storage::integer size)
+	{
+#if defined(USE_MPI)
+		static dynarray<Storage::integer,64> temp(size);
+		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
+		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MAX,comm);
+#else//USE_MPI
+		(void) input;
+#endif//USE_MPI
+	}
+
+	Storage::real Mesh::AggregateMin(Storage::real input)
+	{
+		Storage::real output = input;
+#if defined(USE_MPI)
+		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_REAL_TYPE,MPI_MIN,comm);
+#else //USE_MPI
+		(void) input;
+#endif //USE_MPI
+		return output;
+	}
+	
+	Storage::integer Mesh::AggregateMin(Storage::integer input)
+	{
+		Storage::integer output = input;
+#if defined(USE_MPI)
+		MPI_Allreduce(&input,&output,1,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MIN,comm);
+#else //USE_MPI
+		(void) input;
+#endif //USE_MPI
+		return output;
+	}
+
+	void Mesh::AggregateMin(Storage::real * input, Storage::integer size)
+	{
+#if defined(USE_MPI)
+		static dynarray<Storage::real,64> temp(size);
+		memcpy(temp.data(),input,sizeof(Storage::real)*size);
+		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_REAL_TYPE,MPI_MIN,comm);
+#else//USE_MPI
+		(void) input;
+#endif//USE_MPI
+	}
+
+	void Mesh::AggregateMin(Storage::integer * input, Storage::integer size)
+	{
+#if defined(USE_MPI)
+		static dynarray<Storage::integer,64> temp(size);
+		memcpy(temp.data(),input,sizeof(Storage::integer)*size);
+		MPI_Allreduce(temp.data(),input,size,INMOST_MPI_DATA_INTEGER_TYPE,MPI_MIN,comm);
+#else//USE_MPI
+		(void) input;
+#endif//USE_MPI
+	}
+
 
 	INMOST_MPI_Comm Mesh::GetCommunicator()
 	{

solver.cpp

@@ -960,7 +960,35 @@ namespace INMOST
 #endif
 		for(ind = imbeg; ind < imend; ++ind) //iterate rows of matrix
 			out[ind] = beta * out[ind] + alpha * (*this)[ind].RowVec(x);
-		// outer procedure should update Sout vector, if needed
+		// outer procedure should update out vector, if needed
+	}
+
+	void Solver::Matrix::MatVecTranspose(INMOST_DATA_REAL_TYPE alpha, Solver::Vector & x, INMOST_DATA_REAL_TYPE beta, Solver::Vector & out) const //y = alpha*A*x + beta * y
+	{
+		INMOST_DATA_ENUM_TYPE mbeg, mend;
+		INMOST_DATA_INTEGER_TYPE ind, imbeg, imend;
+		if( out.Empty() )
+		{
+			INMOST_DATA_ENUM_TYPE vbeg,vend;
+			GetInterval(vbeg,vend);
+			out.SetInterval(vbeg,vend);
+		}
+		//CHECK SOMEHOW FOR DEBUG THAT PROVIDED VECTORS ARE OK
+		//~ assert(GetFirstIndex() == out.GetFirstIndex());
+		//~ assert(Size() == out.Size());
+		GetInterval(mbeg,mend);
+		imbeg = mbeg;
+		imend = mend;
+		if( beta ) for(Vector::iterator it = out.Begin(); it != out.End(); ++it) (*it) *= beta;
+#if defined(USE_OMP)
+#pragma omp for private(ind)
+#endif
+		for(ind = imbeg; ind < imend; ++ind)
+		{
+			for(Row::const_iterator it = (*this)[ind].Begin(); it != (*this)[ind].End(); ++it)
+				out[it->first] += alpha * x[ind] * it->second;
+		}
+		// outer procedure should update out vector, if needed
 	}
 	
 	Solver::Matrix::Matrix(std::string _name, INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end, INMOST_MPI_Comm _comm) 

solver_bcgsl.hpp

@@ -25,40 +25,49 @@ namespace INMOST
 	{
 		void dgesvd_(char*,char*,int*,int*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*);
 	}
-#else // SVD adopted from http://stackoverflow.com/questions/3856072/svd-implementation-c answer by Dhairya Malhotra
-	void GivensL(INMOST_DATA_REAL_TYPE * S, const int N, int m, INMOST_DATA_REAL_TYPE a, INMOST_DATA_REAL_TYPE b)
+#else 
+	// SVD adopted from http://stackoverflow.com/questions/3856072/svd-implementation-c answer by Dhairya Malhotra
+	// corrected later to prevent divisions by zero
+
+	static void GivensL(INMOST_DATA_REAL_TYPE * S, const int N, int m, INMOST_DATA_REAL_TYPE a, INMOST_DATA_REAL_TYPE b)
 	{
 		INMOST_DATA_REAL_TYPE r = sqrt(a*a+b*b);
-		INMOST_DATA_REAL_TYPE c = a/r;
-		INMOST_DATA_REAL_TYPE s = -b/r;
-		for(int i=0;i<N;i++)
+		if( r )
 		{
-			INMOST_DATA_REAL_TYPE S0 = S[(m+0)*N+i];
-			INMOST_DATA_REAL_TYPE S1 = S[(m+1)*N+i];
-			S[(m+0)*N + i] += S0*(c-1);
-			S[(m+0)*N + i] += S1*( -s);
-			S[(m+1)*N + i] += S0*(s  );
-			S[(m+1)*N + i] += S1*(c-1);
+			INMOST_DATA_REAL_TYPE c = a/r;
+			INMOST_DATA_REAL_TYPE s = -b/r;
+			for(int i=0;i<N;i++)
+			{
+				INMOST_DATA_REAL_TYPE S0 = S[(m+0)*N+i];
+				INMOST_DATA_REAL_TYPE S1 = S[(m+1)*N+i];
+				S[(m+0)*N + i] += S0*(c-1);
+				S[(m+0)*N + i] += S1*( -s);
+				S[(m+1)*N + i] += S0*(s  );
+				S[(m+1)*N + i] += S1*(c-1);
+			}
 		}
 	}
 
-	void GivensR(INMOST_DATA_REAL_TYPE * S, const int N, int m, INMOST_DATA_REAL_TYPE a, INMOST_DATA_REAL_TYPE b)
+	static void GivensR(INMOST_DATA_REAL_TYPE * S, const int N, int m, INMOST_DATA_REAL_TYPE a, INMOST_DATA_REAL_TYPE b)
 	{
 		INMOST_DATA_REAL_TYPE r = sqrt(a*a+b*b);
-		INMOST_DATA_REAL_TYPE c = a/r;
-		INMOST_DATA_REAL_TYPE s = -b/r;
-		for(int i=0;i<N;i++)
+		if( r )
 		{
-			INMOST_DATA_REAL_TYPE S0 = S[i*N+(m+0)];
-			INMOST_DATA_REAL_TYPE S1 = S[i*N+(m+1)];
-			S[i*N+(m+0)] += S0*(c-1);
-			S[i*N+(m+0)] += S1*( -s);
-			S[i*N+(m+1)] += S0*(s  );
-			S[i*N+(m+1)] += S1*(c-1);
+			INMOST_DATA_REAL_TYPE c = a/r;
+			INMOST_DATA_REAL_TYPE s = -b/r;
+			for(int i=0;i<N;i++)
+			{
+				INMOST_DATA_REAL_TYPE S0 = S[i*N+(m+0)];
+				INMOST_DATA_REAL_TYPE S1 = S[i*N+(m+1)];
+				S[i*N+(m+0)] += S0*(c-1);
+				S[i*N+(m+0)] += S1*( -s);
+				S[i*N+(m+1)] += S0*(s  );
+				S[i*N+(m+1)] += S1*(c-1);
+			}
 		}
 	}
 
-	void svdnxn(INMOST_DATA_REAL_TYPE * A, INMOST_DATA_REAL_TYPE * U, INMOST_DATA_REAL_TYPE * S,  INMOST_DATA_REAL_TYPE * V, const int N)
+	static void svdnxn(INMOST_DATA_REAL_TYPE * A, INMOST_DATA_REAL_TYPE * U, INMOST_DATA_REAL_TYPE * S,  INMOST_DATA_REAL_TYPE * V, const int N)
 	{
 		memset(S,0,sizeof(INMOST_DATA_REAL_TYPE)*N*N);
 		memset(U,0,sizeof(INMOST_DATA_REAL_TYPE)*N*N);
@@ -82,13 +91,21 @@ namespace INMOST
 
 					INMOST_DATA_REAL_TYPE x_inv_norm=0;
 					for(int j=i;j<N;j++) x_inv_norm+= S[j*N+i]*S[j*N+i];
-					x_inv_norm=1/sqrt(x_inv_norm);
+					if( x_inv_norm ) 
+					{
+						x_inv_norm=1/sqrt(x_inv_norm);
 
-					INMOST_DATA_REAL_TYPE alpha=sqrt(1+x1*x_inv_norm);
-					INMOST_DATA_REAL_TYPE beta=x_inv_norm/alpha;
+						INMOST_DATA_REAL_TYPE alpha=sqrt(1+x1*x_inv_norm);
+						INMOST_DATA_REAL_TYPE beta=x_inv_norm/alpha;
+
+						house_vec[i]=-alpha;
+						for(int j=i+1;j<N;j++) house_vec[j]=-beta*S[j*N+i];
+					}
+					else
+					{
+						for(int j=i;j<N;j++) house_vec[j]=0;
+					}
 
-					house_vec[i]=-alpha;
-					for(int j=i+1;j<N;j++) house_vec[j]=-beta*S[j*N+i];
 					if(S[i*N+i]<0) for(int j=i+1;j<N;j++) house_vec[j]=-house_vec[j];
 				}
 				
@@ -116,13 +133,18 @@ namespace INMOST
 
 					INMOST_DATA_REAL_TYPE x_inv_norm=0;
 					for(int j=i+1;j<N;j++) x_inv_norm+=S[i*N+j]*S[i*N+j];
-					x_inv_norm=1/sqrt(x_inv_norm);
+					if( x_inv_norm ) 
+					{
+						x_inv_norm=1/sqrt(x_inv_norm);
 
-					INMOST_DATA_REAL_TYPE alpha=sqrt(1+x1*x_inv_norm);
-					INMOST_DATA_REAL_TYPE beta=x_inv_norm/alpha;
+						INMOST_DATA_REAL_TYPE alpha=sqrt(1+x1*x_inv_norm);
+						INMOST_DATA_REAL_TYPE beta=x_inv_norm/alpha;
 
-					house_vec[i+1]=-alpha;
-					for(int j=i+2;j<N;j++) house_vec[j]=-beta*S[i*N+j];
+						house_vec[i+1]=-alpha;
+						for(int j=i+2;j<N;j++) house_vec[j]=-beta*S[i*N+j];
+					}
+					else for(int j=i+1;j<N;j++) house_vec[j] = 0;
+					
 					if(S[i*N+(i+1)]<0) for(int j=i+2;j<N;j++) house_vec[j]=-house_vec[j];
 				}
 				
@@ -161,13 +183,22 @@ namespace INMOST
 			INMOST_DATA_REAL_TYPE mu=0;
 			{ // Compute mu
 				INMOST_DATA_REAL_TYPE C[2][2];
-				C[0][0]=S[(n-2)*N+(n-2)]*S[(n-2)*N+(n-2)]+S[(n-3)*N+(n-2)]*S[(n-3)*N+(n-2)]; 
+				C[0][0]=S[(n-2)*N+(n-2)]*S[(n-2)*N+(n-2)];
+				if( n-k0 > 2 ) C[0][0] += S[(n-3)*N+(n-2)]*S[(n-3)*N+(n-2)]; 
 				C[0][1]=S[(n-2)*N+(n-2)]*S[(n-2)*N+(n-1)];
 				C[1][0]=S[(n-2)*N+(n-2)]*S[(n-2)*N+(n-1)]; 
 				C[1][1]=S[(n-1)*N+(n-1)]*S[(n-1)*N+(n-1)]+S[(n-2)*N+(n-1)]*S[(n-2)*N+(n-1)];
 				INMOST_DATA_REAL_TYPE b =-(C[0][0]+C[1][1])/2;
 				INMOST_DATA_REAL_TYPE c =  C[0][0]*C[1][1] - C[0][1]*C[1][0];
-				INMOST_DATA_REAL_TYPE d = sqrt(b*b-c);
+				INMOST_DATA_REAL_TYPE d = 0;
+				if( b*b-c > 0 )
+					d= sqrt(b*b-c); //may there be any roundoff problem?
+				else
+				{
+					b =(C[0][0]-C[1][1])/2;
+					c =-C[0][1]*C[1][0];
+					if( b*b-c > 0 ) d = sqrt(b*b-c);
+				}
 				INMOST_DATA_REAL_TYPE lambda1 = -b+d;
 				INMOST_DATA_REAL_TYPE lambda2 = -b-d;
 				INMOST_DATA_REAL_TYPE d1 = lambda1-C[1][1]; d1 = (d1<0?-d1:d1);
@@ -192,7 +223,6 @@ namespace INMOST
 				beta  = S[k*N+(k+2)];
 			}
 		}
-		
 		for(int i=0;i<N;i++)
 		{
 			INMOST_DATA_REAL_TYPE temp;
@@ -207,8 +237,6 @@ namespace INMOST
 				V[j+N*i] = temp;
 			}
 		}
-		
-
 		for(int i=0;i<N;i++) if( S[i*N+i] < 0.0 )
 		{
 			for(int j=0;j<N;j++)
@@ -217,11 +245,10 @@ namespace INMOST
 			}
 			S[i*N+i] *= -1;
 		}
-		
 	}
 #endif //USE_LAPACK_SVD
 #endif //PSEUDOINVERSE
-	int solvenxn(INMOST_DATA_REAL_TYPE * A, INMOST_DATA_REAL_TYPE * x, INMOST_DATA_REAL_TYPE * b, int n, int * order)
+	static int solvenxn(INMOST_DATA_REAL_TYPE * A, INMOST_DATA_REAL_TYPE * x, INMOST_DATA_REAL_TYPE * b, int n, int * order)
 	{
 		INMOST_DATA_REAL_TYPE temp, max;
 		int temp2;

solver_ddpqiluc2.cpp

@@ -3,6 +3,11 @@
 #include "inmost_solver.h"
 #include "solver_ddpqiluc2.hpp"
 
+
+/// \todo
+/// 1. Correct estimator (as in MPTILUC)
+/// 2. Remove premature dropping.
+
 //#define PRINT_DDPQ
 //#define PRINT_HISTOGRAM
 //#define REPORT_ILU