Features

BCGStab(L) solver: perturbation of r_tilde to avoid breakdown.
Possibility to annotate matrix.

inmost_solver.h

@@ -6,6 +6,7 @@
 #include "inmost_common.h"
 //#include "solver_prototypes.hpp"
 
+#define MTX_ALLOW_ANNOTATION
 #define DEFAULT_ADDITIVE_SCHWARTZ_OVERLAP             1
 #define DEFAULT_ABSOLUTE_TOLERANCE                    1.0e-5
 #define DEFAULT_RELATIVE_TOLERANCE                    1.0e-12
@@ -210,7 +211,9 @@ namespace INMOST
 		class Row 
 		{
 		public:
-			
+#if defined(MTX_ALLOW_ANNOTATION)
+			std::string annotation;
+#endif
 			/// Entry of the sparse matrix row.
 			typedef struct entry_s 
 			{
@@ -251,35 +254,25 @@ namespace INMOST
 			bool marker;
 			Entries data;
 		public:
+
+      std::string GetAnnotation();
+      void SetAnnotation(std::string input);
 			void Report() {data.report_addr();}
 			void SetMarker() { marker = true; }
 			void RemMarker() { marker = false; }
 			bool GetMarker() { return marker; }
-			Row() :data() 
-			{
-#if defined(USE_OMP)
-				omp_init_lock(&lock);
-#endif
-				modified_pattern = marker = false;
-			}
-			Row(const Row & other) :marker(other.marker),data(other.data) 
-			{ 
-				//std::cout << __FUNCTION__ << " ";
-				//for(iterator it = Begin(); it != End(); ++it) std::cout << it->first << "," << it->second << " ";
-				//std::cout << std::endl;
-#if defined(USE_OMP)
-				omp_init_lock(&lock);
-#endif
-				modified_pattern = other.modified_pattern; 
-			}
-			Row(entry * pbegin, entry * pend) :data(pbegin, pend) 
-			{ 
+			Row();
+			Row(const Row & other);
+			Row(entry * pbegin, entry * pend);
+      bool HaveLock()
+      {
 #if defined(USE_OMP)
-				omp_init_lock(&lock);
+				return true;
+#else
+        return false;
 #endif
-				modified_pattern = true; marker = false; 
-			}
-			void Lock() 
+      }
+			void Lock()
 			{
 #if defined(USE_OMP)
 				omp_set_lock(&lock); 
@@ -291,8 +284,8 @@ namespace INMOST
 				omp_unset_lock(&lock); 
 #endif
 			}
-			~Row() {}
-			Row & operator = (Row const & other) { data = other.data; marker = other.marker; return *this; }
+			~Row();
+			Row & operator = (Row const & other);
 			/// The operator [] used to fill the sparse matrix row, but not to access individual elements of the row.
 			INMOST_DATA_REAL_TYPE & operator [](INMOST_DATA_ENUM_TYPE i) // use to fill matrix, not to access individual elements
 			{
@@ -325,7 +318,7 @@ namespace INMOST
 			//void           Reserve(INMOST_DATA_ENUM_TYPE num) { data.reserve(num);}
 			/// Clear all data of the current row.
 			void                    Clear() { data.clear(); }
-			void                    Swap(Solver::Row & other) { data.swap(other.data); bool tmp = marker; marker = other.marker; other.marker = tmp; }
+			void                    Swap(Solver::Row & other);
 			/// The size of the sparse row, i.e. the total number of nonzero elements.
 			INMOST_DATA_ENUM_TYPE   Size() const { return static_cast<INMOST_DATA_ENUM_TYPE>(data.size()); }
 			INMOST_DATA_ENUM_TYPE & GetIndex(INMOST_DATA_ENUM_TYPE k) {assert(k < data.size()); return (data.begin()+k)->first;}

solver.cpp

@@ -1295,8 +1295,13 @@ namespace INMOST
 		mtx.precision(15);
 		for(iterator it = Begin(); it != End(); ++it)
 		{
+#if defined(MTX_ALLOW_ANNOTATION)
+      if( it->GetAnnotation() != "" ) mtx << "% " << it->GetAnnotation() << "\n";
+#endif
 			for(Row::iterator jt = it->Begin(); jt != it->End(); ++jt)
-				mtx << row << " " << jt->first+1 << " " << jt->second << std::endl;
+      {
+				mtx << row << " " << jt->first+1 << " " << jt->second << "\n";
+      }
 			++row;
 		}
 #if defined(USE_MPI) && defined(USE_MPI_FILE) // USE_MPI2?
@@ -2874,6 +2879,76 @@ namespace INMOST
     }
     else return 0.0;
   }
+
+  std::string Solver::Row::GetAnnotation() 
+  {
+#if defined(MTX_ALLOW_ANNOTATION)
+    return annotation;
+#else
+    return "";
+#endif
+  }
+  void Solver::Row::SetAnnotation(std::string input) 
+  {
+#if defined(MTX_ALLOW_ANNOTATION)
+    annotation = input;
+#endif
+  }
+  Solver::Row::Row() :data() 
+	{
+#if defined(MTX_ALLOW_ANNOTATION)
+    annotation = "";
+#endif
+#if defined(USE_OMP)
+		omp_init_lock(&lock);
+#endif
+		modified_pattern = marker = false;
+	}
+  Solver::Row::Row(const Row & other) :marker(other.marker),data(other.data) 
+	{ 
+#if defined(MTX_ALLOW_ANNOTATION)
+    annotation = other.annotation;
+#endif
+#if defined(USE_OMP)
+		omp_init_lock(&lock);
+#endif
+		modified_pattern = other.modified_pattern; 
+	}
+  Solver::Row::Row(entry * pbegin, entry * pend) :data(pbegin, pend) 
+	{ 
+#if defined(MTX_ALLOW_ANNOTATION)
+    annotation = "";
+#endif
+#if defined(USE_OMP)
+		omp_init_lock(&lock);
+#endif
+		modified_pattern = true; marker = false; 
+	}
+  Solver::Row & Solver::Row::operator = (Row const & other) 
+  { 
+    data = other.data; 
+    marker = other.marker; 
+#if defined(MTX_ALLOW_ANNOTATION)
+    annotation = other.annotation;
+#endif
+    return *this; 
+  }
+  Solver::Row::~Row()
+  {
+  }
+  void  Solver::Row::Swap(Solver::Row & other) 
+  { 
+    data.swap(other.data); 
+    bool tmp = marker; 
+    marker = other.marker; 
+    other.marker = tmp; 
+#if defined(MTX_ALLOW_ANNOTATION)
+    annotation.swap(other.annotation);
+#endif
+#if defined(USE_OMP)
+    //swap locks?
+#endif
+  }
 }
 #endif
 

solver_bcgsl.hpp

@@ -433,6 +433,7 @@ namespace INMOST
 #endif //PSEUDOINVERSE
 	class BCGSL_solver : public IterativeMethod
 	{
+    INMOST_DATA_REAL_TYPE length, r_tilde_norm;
 		INMOST_DATA_REAL_TYPE rtol, atol, divtol, last_resid;
 		INMOST_DATA_ENUM_TYPE iters, maxits, l, last_it;
 		INMOST_DATA_REAL_TYPE resid;
@@ -608,12 +609,12 @@ namespace INMOST
 				resid += r[0][k]*r[0][k];
 			info->Integrate(&resid,1);
 			//info->ScalarProd(r[0],r[0],vlocbeg,vlocend,resid); //resid = dot(r[0],r[0])
+      last_resid = resid = resid0 = sqrt(resid/rhs_norm); //resid = sqrt(dot(r[0],r[0])
 #if defined(USE_OMP)
 #pragma omp parallel for
 #endif
 			for(INMOST_DATA_INTEGER_TYPE k = ivbeg; k < ivend; ++k) // r_tilde = r[0] / dot(r[0],r[0])
 				r_tilde[k] /= resid;
-			last_resid = resid = resid0 = sqrt(resid/rhs_norm); //resid = sqrt(dot(r[0],r[0])
 			last_it = 0;
 #if defined(REPORT_RESIDUAL)
 			if( info->GetRank() == 0 ) 
@@ -652,6 +653,56 @@ namespace INMOST
 					tt = Timer();
 					for(INMOST_DATA_ENUM_TYPE j = 0; j < l; j++)
 					{
+            if( false )
+            {
+perturbate:
+              //std::cout << "Rescue solver by perturbing orthogonal direction!" << std::endl;
+              //Compute length of the vector
+#if defined(USE_OMP)
+#pragma omp single
+#endif
+              {
+                length = static_cast<INMOST_DATA_REAL_TYPE>(ivlocend-ivlocbeg);
+                r_tilde_norm = 0.0;
+              }
+              info->Integrate(&length,1);
+              //perform perturbation (note that rand() with openmp may give identical sequences of random values
+#if defined(USE_OMP)
+#pragma omp for
+#endif
+              for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
+              {
+                INMOST_DATA_REAL_TYPE unit = 2*static_cast<INMOST_DATA_REAL_TYPE>(rand())/static_cast<INMOST_DATA_REAL_TYPE>(RAND_MAX)-1.0;
+                r_tilde[k] += unit/length;
+              }
+              //compute norm for orthogonal vector
+#if defined(USE_OMP)
+#pragma omp for reduction(+:r_tilde_norm)
+#endif
+              for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
+                r_tilde_norm += r_tilde[k]*r_tilde[k];
+              r_tilde_norm = sqrt(r_tilde_norm);
+              info->Integrate(&r_tilde_norm,1);
+              //normalize orthogonal vector to unity
+#if defined(USE_OMP)
+#pragma omp for
+#endif
+              for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
+                r_tilde[k] /= r_tilde_norm;
+              //Recompute rho1
+              /*
+#if defined(USE_OMP)
+#pragma omp single
+#endif
+	  					rho1 = 0;
+#if defined(USE_OMP)
+#pragma omp for reduction(+:rho1)
+#endif
+						  for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
+							  rho1+= r[j][k]*r_tilde[k];
+						  info->Integrate(&rho1,1);
+              */
+            }
 #if defined(USE_OMP)
 #pragma omp single
 #endif
@@ -662,6 +713,11 @@ namespace INMOST
 						for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k) 
 							rho1+= r[j][k]*r_tilde[k];
 						info->Integrate(&rho1,1);
+            if( fabs(rho1) < 1.0e-50 )
+            {
+              //std::cout << "Asking to perturbate r_tilde since rho1 is too small " << rho1 << std::endl;
+              goto perturbate;
+            }
 						//info->ScalarProd(r[j],r_tilde,vlocbeg,vlocend,rho1); // rho1 = dot(r[j],r_tilde)
 #if defined(USE_OMP)
 #pragma omp single
@@ -683,6 +739,7 @@ namespace INMOST
 
 						if( beta != beta )
 						{
+              //std::cout << "alpha " << alpha << " rho1 " << rho1 << " rho0 " << rho0 << " beta " << beta << std::endl;
 #if defined(USE_OMP)
 #pragma omp single
 #endif
@@ -719,6 +776,12 @@ namespace INMOST
 						info->Integrate(&eta,1);
 						//info->ScalarProd(u[j+1],r_tilde,vlocbeg,vlocend,eta); //eta = dot(u[j+1],r_tilde)
 
+            if( fabs(eta) < 1.0e-50 )
+            {
+              //std::cout << "Asking to perturbate r_tilde since eta is too small " << eta << std::endl;
+              goto perturbate;
+            }
+
 #if defined(USE_OMP)
 #pragma omp single
 #endif