Bug fix and features

Fixed element removal algorithm during modification.
Added features to RowMerger algorithm.
Changed reported number of iterations for BiCGStab(L)

element.cpp

@@ -973,11 +973,7 @@ namespace INMOST
 		GetMeshLink()->SetGeometricType(GetHandle(),t); 
 	}
 
-	Node Element::getAsNode() const {assert(GetElementType() == NODE); return Node(GetMeshLink(),GetHandle());}
-	Edge Element::getAsEdge() const {assert(GetElementType() == EDGE); return Edge(GetMeshLink(),GetHandle());}
-	Face Element::getAsFace() const {assert(GetElementType() == FACE); return Face(GetMeshLink(),GetHandle());} 
-	Cell Element::getAsCell() const {assert(GetElementType() == CELL); return Cell(GetMeshLink(),GetHandle());}
-	ElementSet Element::getAsSet() const {assert(GetElementType() == ESET); return ElementSet(GetMeshLink(),GetHandle());}
+	
 
 }
 #endif

inmost_mesh.h

@@ -449,6 +449,12 @@ namespace INMOST
 		bool                     isValid         () const;
 		__INLINE Mesh *          GetMeshLink     () const {return m_link;}
 		__INLINE HandleType      GetHandle       () const {return handle;}
+    Element                   getAsElement() const;
+    Node                      getAsNode   () const;
+		Edge                      getAsEdge   () const;
+		Face                      getAsFace   () const;
+		Cell                      getAsCell   () const;
+		ElementSet                getAsSet    () const;
 		friend class Mesh;
 	};
 	
@@ -729,11 +735,6 @@ namespace INMOST
 		virtual ElementArray<Edge>  getEdges                (MarkerType mask,bool invert_mask = false) const; //unordered
 		virtual ElementArray<Face>  getFaces                (MarkerType mask,bool invert_mask = false) const; //unordered
 		virtual ElementArray<Cell>  getCells                (MarkerType mask,bool invert_mask = false) const; //unordered
-		Node                        getAsNode               () const;
-		Edge                        getAsEdge               () const;
-		Face                        getAsFace               () const;
-		Cell                        getAsCell               () const;
-		ElementSet                  getAsSet                () const;
 		GeometricType               GetGeometricType        () const;
 		unsigned int                GetElementDimension     () const {return GetGeometricDimension(GetGeometricType());}
 		Status                      GetStatus               () const;

inmost_solver.h

@@ -452,7 +452,7 @@ namespace INMOST
 		/// advised to set total size of the matrix as interval of the
 		/// RowMerger. In future this may change, see todo 2 below.
 		/// \todo
-		/// 1. Add iterators over entries.
+		/// 1. (testing!) Add iterators over entries.
 		/// 2. Implement multiple intervals for distributed computation,
 		///    then in parallel the user may specify additional range of indexes
 		///    for elements that lay on the borders between each pair of processors.
@@ -461,6 +461,30 @@ namespace INMOST
 		public:
 			static const INMOST_DATA_ENUM_TYPE EOL = ENUMUNDEF-1; ///< End of linked list.
 			static const INMOST_DATA_ENUM_TYPE UNDEF = ENUMUNDEF; ///< Value not defined in linked list.
+      class iterator
+      {
+      private:
+        INMOST_DATA_ENUM_TYPE pos;
+        interval< INMOST_DATA_ENUM_TYPE, Row::entry > * LinkedList; ///< Link to associated storage for linked list.
+        iterator(interval< INMOST_DATA_ENUM_TYPE, Row::entry > * pLinkedList) : pos(pLinkedList->get_interval_beg()), LinkedList(pLinkedList) {}
+      public:
+        iterator(const iterator & other) : pos(other.pos), LinkedList(other.LinkedList) {}
+        ~iterator() {}
+        INMOST_DATA_REAL_TYPE & operator *() {return (*LinkedList)[pos].second;}
+        INMOST_DATA_REAL_TYPE operator *() const {return (*LinkedList)[pos].second;}
+        INMOST_DATA_REAL_TYPE * operator ->() {return &(*LinkedList)[pos].second;}
+        const INMOST_DATA_REAL_TYPE * operator ->() const {return &(*LinkedList)[pos].second;}
+        iterator & operator ++(){ pos = (*LinkedList)[pos].first; return *this;}
+        iterator operator ++(int){ iterator ret(LinkedList); ret.pos = (*LinkedList)[pos].first; return ret; }
+        iterator & operator = (const iterator & other) {LinkedList = other.LinkedList; pos = other.pos;}
+        bool operator ==(const iterator & other) const {return LinkedList == other.LinkedList && pos == other.pos;}
+        bool operator !=(const iterator & other) const {return LinkedList != other.LinkedList || pos != other.pos;}
+        bool operator < (const iterator & other) const {return LinkedList == other.LinkedList && pos < other.pos;}
+        bool operator <=(const iterator & other) const {return LinkedList == other.LinkedList && pos <= other.pos;}
+        bool operator > (const iterator & other) const {return LinkedList == other.LinkedList && pos > other.pos;}
+        bool operator >=(const iterator & other) const {return LinkedList == other.LinkedList && pos >= other.pos;}
+        friend class RowMerger;
+      };
 		private:
 			bool Sorted; ///< Contents of linked list should be sorted.
 			INMOST_DATA_ENUM_TYPE Nonzeros; ///< Number of nonzero in linked list.
@@ -497,7 +521,7 @@ namespace INMOST
 			/// Clear linked list.
 			void Clear();
 			/// Add a row with a coefficient into empty linked list.
-			/// This routing should be a bit faster then Solver::RowMerger::AddRow
+			/// This routine should be a bit faster then Solver::RowMerger::AddRow
 			/// for empty linked list. It may result in an unexpected behavior
 			/// for non-empty linked list, asserts will fire in debug mode.
 			/// @param coef Coefficient to multiply row values.
@@ -523,6 +547,17 @@ namespace INMOST
 			//INMOST_DATA_REAL_TYPE ScalarProd(RowMerger & other);
 			/// Get current number of nonzeros from linked list.
 			INMOST_DATA_ENUM_TYPE Size() {return Nonzeros;}
+      /// Retrive/add an entry from/to linked list.
+      /// @param pos Position in the list.
+      INMOST_DATA_REAL_TYPE & operator [] (INMOST_DATA_ENUM_TYPE pos);
+      /// Retrive an entry from linked list.
+      /// \warning
+      /// Will fire an exception if there is no entry.
+      /// @param pos Position in the list.
+      INMOST_DATA_REAL_TYPE operator [] (INMOST_DATA_ENUM_TYPE pos) const;
+
+      iterator Begin() {return iterator(&LinkedList);}
+      iterator End() {iterator ret(&LinkedList); ret.pos = EOL; return ret;}
 		};
 		
 	private:
@@ -565,6 +600,12 @@ namespace INMOST
 		Solver(const Solver & other);// prohibit copy
 		Solver & operator =(Solver const & other); //prohibit assignment
 	public:
+    /// Retrive approximate condition number produced by INNER_MPTILUC.
+    /// The number is cond(L^-1).
+    INMOST_DATA_REAL_TYPE GetConditionNumberL();
+    /// Retrive approximate condition number produced by INNER_MPTILUC.
+    /// The number is cond(U^-1).
+    INMOST_DATA_REAL_TYPE GetConditionNumberU();
 		/// Set the solver parameter of the integer type.
 		/// You can find defaults for parameters in the top of the file inmost_solver.h.
 		///

mesh.cpp

@@ -1661,7 +1661,7 @@ namespace INMOST
 				return ret;
 			}
 		}
-		assert(false); //if you reached here then you either don't release markers (it's a bug) or you should increase MarkerFields const in inmost_mesh.h
+		assert(false); //if you reached here then you either don't release markers (it's your bug) or you should increase MarkerFields const in inmost_mesh.h
 		return InvalidMarker();
 	}
 	void Mesh::ReleaseMarker(MarkerType n)

modify.cpp

@@ -1726,6 +1726,8 @@ public:
 	incident_matrix(Mesh * mesh, InputIterator beg, InputIterator end, typename ElementArray<T>::size_type num_inner)
 	: mesh(mesh), head_column(beg,end), min_loop()
 	{
+    min_loop.SetMeshLink(mesh);
+    temp_loop.SetMeshLink(mesh);
 		//isInputForwardIterators<T,InputIterator>();
 		if( !head_column.empty() )
 		{
@@ -1852,7 +1854,7 @@ public:
 		for(adj_type::size_type it = 0; it < hc.size(); ++it) if( !m->GetMarker(hc[it],hm) )
 			cells.push_back(hc[it]);
 
-		assert(cells.size() == 2);
+		//assert(cells.size() == 2);
 
 		temp.insert(temp.end(),edges.begin(),edges.end());
 
@@ -2130,7 +2132,7 @@ public:
 	{
 		if(!New(h) && Hide(h))
 		{
-			if( GetElementType() != CELL ) //mark all elements that rely on this that they should be deleted
+			if( GetHandleElementType(h) != CELL ) //mark all elements that rely on this that they should be deleted
 			{
 				Element::adj_type & hc = HighConn(h);
 				for(Element::adj_type::size_type it = 0; it < hc.size(); ++it) Delete(hc[it]);

solver.cpp

@@ -115,6 +115,44 @@ namespace INMOST
 
 	Solver::RowMerger::RowMerger() : Sorted(true), Nonzeros(0) {}
 
+  INMOST_DATA_REAL_TYPE & Solver::RowMerger::operator[] (INMOST_DATA_ENUM_TYPE pos)
+  {
+    if( LinkedList[pos+1].first != UNDEF ) return LinkedList[pos+1].second;
+    else 
+    {
+      INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg(), next;
+      if( Sorted )
+      {
+        next = index;
+        while(next < pos+1)
+        {
+          index = next;
+          next = LinkedList[index].first;
+        }
+        assert(index < pos+1);
+        assert(pos+1 < next);
+        ++Nonzeros;
+        LinkedList[index].first = pos+1;
+        LinkedList[pos+1].first = next;
+        return LinkedList[pos+1].second;
+      }
+      else
+      {
+        INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg();
+        ++Nonzeros;
+        LinkedList[pos+1].first = LinkedList[index].first;
+        LinkedList[index].first = pos+1;
+        return LinkedList[pos+1].second;
+      }
+    }
+  }
+
+  INMOST_DATA_REAL_TYPE Solver::RowMerger::operator[] (INMOST_DATA_ENUM_TYPE pos) const
+  {
+    if( LinkedList[pos+1].first != UNDEF ) return LinkedList[pos+1].second;
+    else throw -1;
+  }
+
 	Solver::RowMerger::RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted) 
 				: Sorted(Sorted), Nonzeros(0), LinkedList(interval_begin,interval_end+1,Row::make_entry(UNDEF,0.0)) 
 	{
@@ -158,6 +196,7 @@ namespace INMOST
 		{
 			j = LinkedList[i].first;
 			LinkedList[i].first = UNDEF;
+      LinkedList[i].second = 0.0;
 			i = j;
 		}
 		Nonzeros = 0;
@@ -2818,6 +2857,23 @@ namespace INMOST
 	{
 		return last_resid;
 	}
+
+  INMOST_DATA_REAL_TYPE Solver::GetConditionNumberL()
+  {
+    if(_pack == INNER_MPTILUC )
+    {
+      return static_cast<IterativeMethod *>(solver_data)->RealParameter(":condition_number_L");
+    }
+    else return 0.0;
+  }
+  INMOST_DATA_REAL_TYPE Solver::GetConditionNumberU()
+  {
+    if(_pack == INNER_MPTILUC )
+    {
+      return static_cast<IterativeMethod *>(solver_data)->RealParameter(":condition_number_U");
+    }
+    else return 0.0;
+  }
 }
 #endif
 

solver_bcgsl.hpp

@@ -223,7 +223,7 @@ namespace INMOST
 							h = 1.0 / h;
 							c = g * h;
 							s = (- f * h);
-							for (j = 0; j < n; j++) 
+							for (j = nm < 0 ? 1 : 0; j < n; j++) 
 							{
 								y = u[j*n+nm];
 								z = u[j*n+i];
@@ -775,6 +775,7 @@ namespace INMOST
 #endif
 						resid = sqrt(resid/rhs_norm); // resid = sqrt(dot(r[j],r[j]))
 
+            last_it++;
 					
 						if( resid < atol || resid < rtol*resid0 ) 
 						{
@@ -787,6 +788,8 @@ namespace INMOST
 							break;
 						}
 						ApplyOperator(r[j],r[j+1]); // r[j+1] = A*R*r[j]		
+
+            
 					}
 
 					if( halt ) break;
@@ -808,11 +811,30 @@ namespace INMOST
 #endif
 							for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k)
 								tau_sum += r[j][k]*r[m][k];
+
+              if( fabs(tau_sum) > 1.0e+100 )
+              {
 #if defined(USE_OMP)
 #pragma omp single
 #endif
+							  reason = "multiplier(tau) is too large";
+                halt = true;
+              }
+              if( tau_sum != tau_sum )
+              {
+#if defined(USE_OMP)
+#pragma omp single
+#endif
+							  reason = "multiplier(tau) is NaN";
+                halt = true;
+              }
+#if defined(USE_OMP)
+#pragma omp single
+#endif
+              {
 							  tau[(j-1) + (m-1)*size] = tau[(m-1) + (j-1)*size] = tau_sum;
               }
+						}
 #if defined(USE_OMP)
 #pragma omp single
 #endif
@@ -822,11 +844,33 @@ namespace INMOST
 #endif
 						for(INMOST_DATA_INTEGER_TYPE k = ivlocbeg; k < ivlocend; ++k)
 							sigma_sum += r[0][k]*r[j][k];
+
+            if( fabs(sigma_sum) > 1.0e+100 )
+            {
+#if defined(USE_OMP)
+#pragma omp single
+#endif
+							reason = "multiplier(sigma) is too large";
+              halt = true;
+            }
+            if( sigma_sum != sigma_sum )
+            {
+#if defined(USE_OMP)
+#pragma omp single
+#endif
+							reason = "multiplier(sigma) is NaN";
+              halt = true;
+            }
 #if defined(USE_OMP)
 #pragma omp single
 #endif
+            {
 						  sigma[j-1] = sigma_sum;
             }
+					}
+
+          if( halt ) break;
+
 #if defined(CONVEX_COMBINATION)
 #if defined(USE_OMP)
 #pragma omp single
@@ -1067,7 +1111,7 @@ namespace INMOST
 						}
 					}
 					*/
-					last_it = i+1;
+					//last_it = l+1;
 					{
 #if defined(USE_OMP)
 #pragma omp single

storage.cpp

@@ -231,6 +231,12 @@ namespace INMOST
 	{
 		return GetMeshLink()->DataLocalID(GetHandle());
 	}
+  Element Storage::getAsElement() const {assert(GetElementType() & (NODE | EDGE | FACE | CELL | ESET) ); return Element(GetMeshLink(), GetHandle());}
+  Node Storage::getAsNode() const {assert(GetElementType() == NODE); return Node(GetMeshLink(),GetHandle());}
+	Edge Storage::getAsEdge() const {assert(GetElementType() == EDGE); return Edge(GetMeshLink(),GetHandle());}
+	Face Storage::getAsFace() const {assert(GetElementType() == FACE); return Face(GetMeshLink(),GetHandle());} 
+	Cell Storage::getAsCell() const {assert(GetElementType() == CELL); return Cell(GetMeshLink(),GetHandle());}
+	ElementSet Storage::getAsSet() const {assert(GetElementType() == ESET); return ElementSet(GetMeshLink(),GetHandle());}
 }
 
 #endif