MatrixSolver::solvewithsaveOriginal

MbDCode/AbsConstraint.cpp

@@ -7,7 +7,7 @@ AbsConstraint::AbsConstraint() {}
 
 AbsConstraint::AbsConstraint(const char* str) : Constraint(str) {}
 
-AbsConstraint::AbsConstraint(int i)
+AbsConstraint::AbsConstraint(size_t i)
 {
     axis = i;
 }
@@ -27,3 +27,14 @@ void MbD::AbsConstraint::calcPostDynCorrectorIteration()
         aG = static_cast<PartFrame*>(owner)->qE->at(axis - 3);
     }
 }
+
+void MbD::AbsConstraint::useEquationNumbers()
+{
+    iqXminusOnePlusAxis = static_cast<PartFrame*>(owner)->iqX - 1 + axis;
+}
+
+void MbD::AbsConstraint::fillPosICJacob(SpMatDsptr mat)
+{
+    mat->atijplusNumber(iG, iqXminusOnePlusAxis, 1.0);
+    mat->atijplusNumber(iqXminusOnePlusAxis, iG, 1.0);
+}

MbDCode/AbsConstraint.h

@@ -7,12 +7,14 @@ namespace MbD {
     public:
         AbsConstraint();
         AbsConstraint(const char* str);
-        AbsConstraint(int axis);
+        AbsConstraint(size_t axis);
         void initialize();
         void calcPostDynCorrectorIteration() override;
+        void useEquationNumbers() override;
+        void fillPosICJacob(SpMatDsptr mat) override;
 
-        int axis;
-        int iqXminusOnePlusAxis;
+        size_t axis = -1;
+        size_t iqXminusOnePlusAxis = -1;
     };
 }
 

MbDCode/AnyPosICNewtonRaphson.cpp

@@ -1,7 +1,49 @@
 #include "AnyPosICNewtonRaphson.h"
+#include "SystemSolver.h"
+#include "Item.h"
+
+using namespace MbD;
 
 void MbD::AnyPosICNewtonRaphson::initialize()
 {
 	NewtonRaphson::initialize();
 	nSingularMatrixError = 0;
 }
+
+void MbD::AnyPosICNewtonRaphson::initializeGlobally()
+{
+	SystemNewtonRaphson::initializeGlobally();
+	system->partsJointsMotionsDo([&](std::shared_ptr<Item> item) {
+		item->fillqsu(qsuOld);
+		item->fillqsuWeights(qsuWeights);
+		item->fillqsulam(x);
+		});
+}
+
+void MbD::AnyPosICNewtonRaphson::createVectorsAndMatrices()
+{
+	qsuOld = std::make_shared<FullColumn<double>>(nqsu);
+	qsuWeights = std::make_shared<DiagonalMatrix<double>>(nqsu);
+	SystemNewtonRaphson::createVectorsAndMatrices();
+}
+
+void MbD::AnyPosICNewtonRaphson::fillY()
+{
+	auto newMinusOld = qsuOld->negated();
+	newMinusOld->equalSelfPlusFullColumnAt(x, 0);
+	y->zeroSelf();
+	y->atiminusFullColumn(0, (qsuWeights->timesFullColumn(newMinusOld)));
+	system->partsJointsMotionsDo([&](std::shared_ptr<Item> item) { item->fillPosICError(y); });
+}
+
+void MbD::AnyPosICNewtonRaphson::fillPyPx()
+{
+	pypx->zeroSelf();
+	pypx->atijminusDiagonalMatrix(0, 0, qsuWeights);
+	system->partsJointsMotionsDo([&](std::shared_ptr<Item> item) {item->fillPosICJacob(pypx); });
+}
+
+void MbD::AnyPosICNewtonRaphson::passRootToSystem()
+{
+	system->partsJointsMotionsDo([&](std::shared_ptr<Item> item) {item->setqsulam(x); });
+}

MbDCode/AnyPosICNewtonRaphson.h

@@ -1,17 +1,25 @@
 #pragma once
 
 #include "PosNewtonRaphson.h"
+#include "DiagonalMatrix.h"
 
 namespace MbD {
+
     class AnyPosICNewtonRaphson : public PosNewtonRaphson
     {
         //nqsu qsuOld qsuWeights nSingularMatrixError 
     public:
         void initialize() override;
+        void initializeGlobally() override;
+        void createVectorsAndMatrices() override;
+        void fillY() override;
+        void fillPyPx() override;
+        void passRootToSystem() override;
 
-        int nqsu;
-        std::shared_ptr<FullColumn<double>> qsuOld, qsuWeights;
-        int nSingularMatrixError;
+        size_t nqsu = -1;
+        std::shared_ptr<FullColumn<double>> qsuOld;
+        std::shared_ptr<DiagonalMatrix<double>> qsuWeights;
+        size_t nSingularMatrixError = -1;
     };
 }
 

MbDCode/Array.cpp

@@ -1 +1,3 @@
 #include "Array.h"
+
+using namespace MbD;

MbDCode/Array.h

@@ -1,7 +1,10 @@
 #pragma once
 #include <vector>
 #include <memory>
-#include <type_traits>
+//#include <type_traits>
+#include <cmath>
+#include <cassert>
+
 
 namespace MbD {
 	template <typename T>
@@ -13,9 +16,10 @@ namespace MbD {
 		Array(size_t count, const T& value) : std::vector<T>(count, value) {}
 		Array(std::initializer_list<T> list) : std::vector<T>{ list } {}
 		void copyFrom(std::shared_ptr<Array<T>> x);
-		void zeroSelf();
-		double sumOfSquares();
-		double sumOfSquaresOfVector();
+		virtual void zeroSelf() = 0;
+		virtual double sumOfSquares() = 0;
+		double rootMeanSquare();
+		virtual size_t numberOfElements() = 0;
 	};
 	template<typename T>
 	inline void Array<T>::copyFrom(std::shared_ptr<Array<T>> x)
@@ -24,37 +28,16 @@ namespace MbD {
 			this->at(i) = x->at(i);
 		}
 	}
-	template <>
-	inline void Array<double>::zeroSelf() {
-		for (size_t i = 0; i < this->size(); i++) {
-			this->at(i) = 0.0;;
-		}
-	}
+	//template <>
+	//inline void Array<double>::zeroSelf() {
+	//	for (size_t i = 0; i < this->size(); i++) {
+	//		this->at(i) = 0.0;;
+	//	}
+	//}
 	template<typename T>
-	inline double Array<T>::sumOfSquares()
+	inline double Array<T>::rootMeanSquare()
 	{
-		if (std::is_arithmetic<T>) {
-			return this->sumOfSquaresOfVector();
-		}
-		else {
-			double sum = 0.0;
-			for (size_t i = 0; i < this->size(); i++)
-			{
-				sum += this->at(i)->sumOfSquares();
-			}
-			return sum;
-		}
-	}
-	template<typename T>
-	inline double Array<T>::sumOfSquaresOfVector()
-	{
-		double sum = 0.0;
-		for (size_t i = 0; i < this->size(); i++)
-		{
-			double element = this->at(i);
-			sum += element * element;
-		}
-		return sum;
+		return std::sqrt(this->sumOfSquares() / this->numberOfElements());
 	}
 	using ListD = std::initializer_list<double>;
 	using ListListD = std::initializer_list<std::initializer_list<double>>;

MbDCode/AtPointConstraintIJ.cpp

@@ -4,7 +4,7 @@
 
 using namespace MbD;
 
-AtPointConstraintIJ::AtPointConstraintIJ(EndFrmcptr frmi, EndFrmcptr frmj, int axisi) :
+AtPointConstraintIJ::AtPointConstraintIJ(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi) :
 	ConstraintIJ(frmi, frmj), axis(axisi)
 {
 }
@@ -50,3 +50,9 @@ MbD::ConstraintType MbD::AtPointConstraintIJ::type()
 {
 	return MbD::displacement;
 }
+
+void MbD::AtPointConstraintIJ::postPosICIteration()
+{
+	riIeJeO->postPosICIteration();
+	Item::postPosICIteration();
+}

MbDCode/AtPointConstraintIJ.h

@@ -9,7 +9,7 @@ namespace MbD {
     {
         //axis riIeJeO 
     public:
-        AtPointConstraintIJ(EndFrmcptr frmi, EndFrmcptr frmj, int axisi);
+        AtPointConstraintIJ(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi);
         void initialize();
         void initializeLocally() override;
         void initializeGlobally() override;
@@ -18,8 +18,9 @@ namespace MbD {
         void calcPostDynCorrectorIteration() override;
         void prePosIC() override;
         MbD::ConstraintType type() override;
+        void postPosICIteration() override;
 
-        int axis;
+        size_t axis;
         std::shared_ptr<DispCompIecJecO> riIeJeO;
     };
 }

MbDCode/AtPointConstraintIqcJc.cpp

@@ -5,7 +5,7 @@
 
 using namespace MbD;
 
-AtPointConstraintIqcJc::AtPointConstraintIqcJc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi) :
+AtPointConstraintIqcJc::AtPointConstraintIqcJc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi) :
 	AtPointConstraintIJ(frmi, frmj, axisi)
 {
 }

MbDCode/AtPointConstraintIqcJc.h

@@ -7,7 +7,7 @@ namespace MbD {
 	{
 		//pGpEI ppGpEIpEI iqXIminusOnePlusAxis iqEI 
 	public:
-		AtPointConstraintIqcJc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi);
+		AtPointConstraintIqcJc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi);
 		void initializeGlobally() override;
 		void initriIeJeO() override;
 		void calcPostDynCorrectorIteration() override;
@@ -15,8 +15,8 @@ namespace MbD {
 
 		FRowDsptr pGpEI;
 		FMatDsptr ppGpEIpEI;
-		int iqXIminusOnePlusAxis = -1;
-		int iqEI = -1;
+		size_t iqXIminusOnePlusAxis = -1;
+		size_t iqEI = -1;
 	};
 }
 

MbDCode/AtPointConstraintIqcJqc.cpp

@@ -5,7 +5,7 @@
 
 using namespace MbD;
 
-AtPointConstraintIqcJqc::AtPointConstraintIqcJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi) :
+AtPointConstraintIqcJqc::AtPointConstraintIqcJqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi) :
 	AtPointConstraintIqcJc(frmi, frmj, axisi)
 {
 }

MbDCode/AtPointConstraintIqcJqc.h

@@ -7,7 +7,7 @@ namespace MbD {
     {
         //pGpEJ ppGpEJpEJ iqXJminusOnePlusAxis iqEJ 
     public:
-        AtPointConstraintIqcJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi);
+        AtPointConstraintIqcJqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi);
         void initializeGlobally() override;
         void initriIeJeO() override;
         void calcPostDynCorrectorIteration() override;
@@ -15,7 +15,7 @@ namespace MbD {
 
         FRowDsptr pGpEJ;
         FMatDsptr ppGpEJpEJ;
-        int iqXJminusOnePlusAxis, iqEJ;
+        size_t iqXJminusOnePlusAxis = -1, iqEJ = -1;
     };
 }
 

MbDCode/AtPointConstraintIqctJqc.cpp

@@ -4,7 +4,7 @@
 
 using namespace MbD;
 
-AtPointConstraintIqctJqc::AtPointConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi) :
+AtPointConstraintIqctJqc::AtPointConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi) :
 	AtPointConstraintIqcJqc(frmi, frmj, axisi)
 {
 }

MbDCode/AtPointConstraintIqctJqc.h

@@ -7,7 +7,7 @@ namespace MbD {
     {
         //pGpt ppGpEIpt ppGptpt 
     public:
-        AtPointConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi);
+        AtPointConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi);
         void initializeGlobally() override;
         void initriIeJeO() override;
         void calcPostDynCorrectorIteration() override;

MbDCode/BasicIntegrator.cpp

@@ -1 +1,3 @@
 #include "BasicIntegrator.h"
+
+using namespace MbD;

MbDCode/BasicQuasiIntegrator.cpp

@@ -1 +1,3 @@
 #include "BasicQuasiIntegrator.h"
+
+using namespace MbD;

MbDCode/BasicUserFunction.cpp

@@ -1 +1,3 @@
 #include "BasicUserFunction.h"
+
+using namespace MbD;

MbDCode/CREATE.cpp

@@ -1 +1,3 @@
 #include "CREATE.h"
+
+using namespace MbD;

MbDCode/CREATE.h

@@ -23,7 +23,7 @@ namespace MbD {
 			inst->initialize();
 			return inst;
 		}
-		static std::shared_ptr<T> With(int n) {
+		static std::shared_ptr<T> With(size_t n) {
 			auto inst = std::make_shared<T>(n);
 			inst->initialize();
 			return inst;
@@ -33,17 +33,17 @@ namespace MbD {
 			inst->initialize();
 			return inst;
 		}
-		static std::shared_ptr<T> With(EndFrmcptr frmi, EndFrmcptr frmj, int axis) {
+		static std::shared_ptr<T> With(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis) {
 			auto inst = std::make_shared<T>(frmi, frmj, axis);
 			inst->initialize();
 			return inst;
 		}
-		static std::shared_ptr<T> With(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk) {
+		static std::shared_ptr<T> With(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk) {
 			auto inst = std::make_shared<T>(frmi, frmj, frmk, axisk);
 			inst->initialize();
 			return inst;
 		}
-		static std::shared_ptr<Constraint> ConstraintWith(EndFrmcptr frmi, EndFrmcptr frmj, int axis) {
+		static std::shared_ptr<Constraint> ConstraintWith(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis) {
 			std::shared_ptr<Constraint> inst;
 			std::string str = typeid(T(frmi, frmj, axis)).name();
 			if (str == "class MbD::AtPointConstraintIJ") {
@@ -65,12 +65,12 @@ namespace MbD {
 			inst->initialize();
 			return inst;
 		}
-		static std::shared_ptr<T> With(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) {
+		static std::shared_ptr<T> With(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) {
 			auto inst = std::make_shared<T>(frmi, frmj, axisi, axisj);
 			inst->initialize();
 			return inst;
 		}
-		static std::shared_ptr<Constraint> ConstraintWith(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) {
+		static std::shared_ptr<Constraint> ConstraintWith(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) {
 			std::shared_ptr<Constraint> inst;
 			std::string str = typeid(T(frmi, frmj, axisi, axisj)).name();
 			if (str == "class MbD::DirectionCosineConstraintIJ") {

MbDCode/Constraint.cpp

@@ -55,7 +55,7 @@ void MbD::Constraint::fillDispConstraints(std::shared_ptr<Constraint>sptr, std::
 }
 void MbD::Constraint::fillPerpenConstraints(std::shared_ptr<Constraint>sptr, std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> perpenConstraints)
 {
-	if (this->type() == MbD::displacement) {
+	if (this->type() == MbD::perpendicular) {
 		perpenConstraints->push_back(sptr);
 	}
 }
@@ -64,3 +64,13 @@ MbD::ConstraintType MbD::Constraint::type()
 {
 	return MbD::essential;
 }
+
+void MbD::Constraint::fillqsulam(FColDsptr col)
+{
+	col->at(iG) = lam;
+}
+
+void MbD::Constraint::setqsulam(FColDsptr col)
+{
+	lam = col->at(iG);
+}

MbDCode/Constraint.h

@@ -20,12 +20,13 @@ namespace MbD {
 		virtual void fillDispConstraints(std::shared_ptr<Constraint>sptr, std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> dispConstraints);
 		virtual void fillPerpenConstraints(std::shared_ptr<Constraint>sptr, std::shared_ptr<std::vector<std::shared_ptr<Constraint>>> perpenConstraints);
 		virtual MbD::ConstraintType type();
-
+		void fillqsulam(FColDsptr col) override;
+		void setqsulam(FColDsptr col) override;
 		
-		int iG;
+		size_t iG;
 		double aG;		//Constraint function
 		double lam;		//Lambda is Lagrange Multiplier
-		Item* owner;	//A Joint or PartFrame owns the constraint
+		Item* owner;	//A Joint or PartFrame owns the constraint.  //Use raw pointer when pointing backwards.
 	};
 }
 

MbDCode/DiagonalMatrix.cpp

@@ -1 +1,3 @@
 #include "DiagonalMatrix.h"
+
+using namespace MbD;

MbDCode/DiagonalMatrix.h

@@ -1,13 +1,72 @@
 #pragma once
+
 #include "Array.h"
+#include "FullColumn.h"
+
 namespace MbD {
-    template <typename T>
-    class DiagonalMatrix : public Array<T>
-    {
-        //
-    public:
-        DiagonalMatrix(size_t count) : Array<T>(count) {}
-        DiagonalMatrix(std::initializer_list<T> list) : Array<T>{ list } {}
-    };
+	template <typename T>
+	class DiagonalMatrix : public Array<T>
+	{
+		//
+	public:
+		DiagonalMatrix(size_t count) : Array<T>(count) {}
+		DiagonalMatrix(std::initializer_list<T> list) : Array<T>{ list } {}
+		void atiputDiagonalMatrix(size_t i, std::shared_ptr < DiagonalMatrix<T>> diagMat);
+		std::shared_ptr<FullColumn<T>> timesFullColumn(std::shared_ptr<FullColumn<T>> fullCol);
+		size_t nrow() {
+			return this->size();
+		}
+		size_t ncol() {
+			return this->size();
+		}
+		double sumOfSquares() override;
+		size_t numberOfElements() override;
+		void zeroSelf() override;
+	};
+	template<typename T>
+	inline void DiagonalMatrix<T>::atiputDiagonalMatrix(size_t i, std::shared_ptr<DiagonalMatrix<T>> diagMat)
+	{
+		for (size_t ii = 0; ii < diagMat->size(); ii++)
+		{
+			this->at(i + ii) = diagMat->at(ii);
+		}
+	}
+	template<typename T>
+	inline std::shared_ptr<FullColumn<T>> DiagonalMatrix<T>::timesFullColumn(std::shared_ptr<FullColumn<T>> fullCol)
+	{
+		//"a*b = a(i,j)b(j) sum j."
+
+		auto nrow = this->size();
+		auto answer = std::make_shared<FullColumn<T>>(nrow);
+		for (size_t i = 0; i < nrow; i++)
+		{
+			answer->at(i) = this->at(i) * fullCol->at(i);
+		}
+		return answer;
+	}
+	template<>
+	inline double DiagonalMatrix<double>::sumOfSquares()
+	{
+		double sum = 0.0;
+		for (size_t i = 0; i < this->size(); i++)
+		{
+			double element = this->at(i);
+			sum += element * element;
+		}
+		return sum;
+	}
+	template<typename T>
+	inline size_t DiagonalMatrix<T>::numberOfElements()
+	{
+		auto n = this->size();
+		return n * n;
+	}
+	template<>
+	inline void DiagonalMatrix<double>::zeroSelf()
+	{
+		for (size_t i = 0; i < this->size(); i++) {
+			this->at(i) = 0.0;;
+		}
+	}
 }
 

MbDCode/DirectionCosineConstraintIJ.cpp

@@ -5,7 +5,7 @@
 
 using namespace MbD;
 
-DirectionCosineConstraintIJ::DirectionCosineConstraintIJ(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) :
+DirectionCosineConstraintIJ::DirectionCosineConstraintIJ(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) :
 	ConstraintIJ(frmi, frmj), axisI(axisi), axisJ(axisj)
 {
 }
@@ -47,6 +47,12 @@ void MbD::DirectionCosineConstraintIJ::prePosIC()
 	Constraint::prePosIC();
 }
 
+void MbD::DirectionCosineConstraintIJ::postPosICIteration()
+{
+	aAijIeJe->postPosICIteration();
+	Item::postPosICIteration();
+}
+
 MbD::ConstraintType MbD::DirectionCosineConstraintIJ::type()
 {
 	return MbD::perpendicular;

MbDCode/DirectionCosineConstraintIJ.h

@@ -9,7 +9,7 @@ namespace MbD {
     {
         //axisI axisJ aAijIeJe 
     public:
-        DirectionCosineConstraintIJ(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj);
+        DirectionCosineConstraintIJ(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj);
         void initialize();
         void initializeLocally() override;
         void initializeGlobally() override;
@@ -17,9 +17,10 @@ namespace MbD {
         void postInput() override;
         void calcPostDynCorrectorIteration() override;
         void prePosIC() override;
+        void postPosICIteration() override;
         MbD::ConstraintType type() override;
 
-        int axisI, axisJ;
+        size_t axisI, axisJ;
         std::shared_ptr<DirectionCosineIecJec> aAijIeJe;
     };
 }

MbDCode/DirectionCosineConstraintIqcJc.cpp

@@ -5,7 +5,7 @@
 
 using namespace MbD;
 
-DirectionCosineConstraintIqcJc::DirectionCosineConstraintIqcJc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) :
+DirectionCosineConstraintIqcJc::DirectionCosineConstraintIqcJc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) :
 	DirectionCosineConstraintIJ(frmi, frmj, axisi, axisj)
 {
 }

MbDCode/DirectionCosineConstraintIqcJc.h

@@ -7,14 +7,14 @@ namespace MbD {
     {
         //pGpEI ppGpEIpEI iqEI 
     public:
-        DirectionCosineConstraintIqcJc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj);
+        DirectionCosineConstraintIqcJc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj);
         void initaAijIeJe() override;
         void calcPostDynCorrectorIteration() override;
         void useEquationNumbers() override;
 
         FRowDsptr pGpEI;
         FMatDsptr ppGpEIpEI;
-        int iqEI;
+        size_t iqEI = -1;
     };
 }
 

MbDCode/DirectionCosineConstraintIqcJqc.cpp

@@ -5,7 +5,7 @@
 
 using namespace MbD;
 
-DirectionCosineConstraintIqcJqc::DirectionCosineConstraintIqcJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) :
+DirectionCosineConstraintIqcJqc::DirectionCosineConstraintIqcJqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) :
 	DirectionCosineConstraintIqcJc(frmi, frmj, axisi, axisj)
 {
 }

MbDCode/DirectionCosineConstraintIqcJqc.h

@@ -7,7 +7,7 @@ namespace MbD {
     {
         //pGpEJ ppGpEIpEJ ppGpEJpEJ iqEJ 
     public:
-        DirectionCosineConstraintIqcJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj);
+        DirectionCosineConstraintIqcJqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj);
         void initaAijIeJe() override;
         void calcPostDynCorrectorIteration() override;
         void useEquationNumbers() override;
@@ -15,7 +15,7 @@ namespace MbD {
         FRowDsptr pGpEJ;
         FMatDsptr ppGpEIpEJ;
         FMatDsptr ppGpEJpEJ;
-        int iqEJ;
+        size_t iqEJ = -1;
     };
 }
 

MbDCode/DirectionCosineConstraintIqctJqc.cpp

@@ -4,7 +4,7 @@
 
 using namespace MbD;
 
-DirectionCosineConstraintIqctJqc::DirectionCosineConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) :
+DirectionCosineConstraintIqctJqc::DirectionCosineConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) :
 	DirectionCosineConstraintIqcJqc(frmi, frmj, axisi, axisj)
 {
 }

MbDCode/DirectionCosineConstraintIqctJqc.h

@@ -7,7 +7,7 @@ namespace MbD {
 	{
 		//pGpt ppGpEIpt ppGpEJpt ppGptpt 
 	public:
-		DirectionCosineConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj);
+		DirectionCosineConstraintIqctJqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj);
 		void initaAijIeJe() override;
 		MbD::ConstraintType type() override;
 

MbDCode/DirectionCosineIecJec.cpp

@@ -9,7 +9,7 @@ DirectionCosineIecJec::DirectionCosineIecJec()
 {
 }
 
-DirectionCosineIecJec::DirectionCosineIecJec(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) :
+DirectionCosineIecJec::DirectionCosineIecJec(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) :
 	KinematicIeJe(frmi, frmj), axisI(axisi), axisJ(axisj)
 {
 

MbDCode/DirectionCosineIecJec.h

@@ -12,10 +12,10 @@ namespace MbD {
         //aAijIeJe axisI axisJ aAjOIe aAjOJe 
     public:
         DirectionCosineIecJec();
-        DirectionCosineIecJec(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj);
+        DirectionCosineIecJec(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj);
         void calcPostDynCorrectorIteration() override;
 
-        int axisI, axisJ;   //0, 1, 2 = x, y, z
+        size_t axisI, axisJ;   //0, 1, 2 = x, y, z
         double aAijIeJe;
         std::shared_ptr<FullColumn<double>> aAjOIe, aAjOJe;
     };

MbDCode/DirectionCosineIeqcJec.cpp

@@ -7,7 +7,7 @@ DirectionCosineIeqcJec::DirectionCosineIeqcJec()
 {
 }
 
-DirectionCosineIeqcJec::DirectionCosineIeqcJec(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) :
+DirectionCosineIeqcJec::DirectionCosineIeqcJec(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) :
 	DirectionCosineIecJec(frmi, frmj, axisi, axisj)
 {
 }

MbDCode/DirectionCosineIeqcJec.h

@@ -8,7 +8,7 @@ namespace MbD {
         //pAijIeJepEI ppAijIeJepEIpEI pAjOIepEIT ppAjOIepEIpEI 
     public:
         DirectionCosineIeqcJec();
-        DirectionCosineIeqcJec(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj);
+        DirectionCosineIeqcJec(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj);
         void initialize();
         void initializeGlobally();
         void calcPostDynCorrectorIteration() override;

MbDCode/DirectionCosineIeqcJeqc.cpp

@@ -7,7 +7,7 @@ DirectionCosineIeqcJeqc::DirectionCosineIeqcJeqc()
 {
 }
 
-DirectionCosineIeqcJeqc::DirectionCosineIeqcJeqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) :
+DirectionCosineIeqcJeqc::DirectionCosineIeqcJeqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) :
 	DirectionCosineIeqcJec(frmi, frmj, axisi, axisj)
 {
 }

MbDCode/DirectionCosineIeqcJeqc.h

@@ -8,7 +8,7 @@ namespace MbD {
         //pAijIeJepEJ ppAijIeJepEIpEJ ppAijIeJepEJpEJ pAjOJepEJT ppAjOJepEJpEJ 
     public:
         DirectionCosineIeqcJeqc();
-        DirectionCosineIeqcJeqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj);
+        DirectionCosineIeqcJeqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj);
         void initialize();
         void initializeGlobally();
         void calcPostDynCorrectorIteration() override;

MbDCode/DirectionCosineIeqctJeqc.cpp

@@ -7,7 +7,7 @@ DirectionCosineIeqctJeqc::DirectionCosineIeqctJeqc()
 {
 }
 
-DirectionCosineIeqctJeqc::DirectionCosineIeqctJeqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj) :
+DirectionCosineIeqctJeqc::DirectionCosineIeqctJeqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj) :
 	DirectionCosineIeqcJeqc(frmi, frmj, axisi, axisj)
 {
 }

MbDCode/DirectionCosineIeqctJeqc.h

@@ -8,7 +8,7 @@ namespace MbD {
         //pAijIeJept ppAijIeJepEIpt ppAijIeJepEJpt ppAijIeJeptpt 
     public:
         DirectionCosineIeqctJeqc();
-        DirectionCosineIeqctJeqc(EndFrmcptr frmi, EndFrmcptr frmj, int axisi, int axisj);
+        DirectionCosineIeqctJeqc(EndFrmcptr frmi, EndFrmcptr frmj, size_t axisi, size_t axisj);
         void initialize();
         void initializeGlobally();
         void calcPostDynCorrectorIteration() override;

MbDCode/DispCompIecJecKec.cpp

@@ -1,10 +1,12 @@
 #include "DispCompIecJecKec.h"
 
+using namespace MbD;
+
 MbD::DispCompIecJecKec::DispCompIecJecKec()
 {
 }
 
-MbD::DispCompIecJecKec::DispCompIecJecKec(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk): KinematicIeJe(frmi, frmj), efrmK(frmk), axisK(axisk)
+MbD::DispCompIecJecKec::DispCompIecJecKec(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk): KinematicIeJe(frmi, frmj), efrmK(frmk), axisK(axisk)
 {
 }
 

MbDCode/DispCompIecJecKec.h

@@ -8,11 +8,11 @@ namespace MbD {
         //efrmK axisK riIeJeKe aAjOKe rIeJeO 
     public:
         DispCompIecJecKec();
-        DispCompIecJecKec(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk);
+        DispCompIecJecKec(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk);
         
         virtual double value();
         EndFrmcptr efrmK;
-        int axisK;
+        size_t axisK;
         double riIeJeKe;
         FColDsptr aAjOKe;
         FColDsptr rIeJeO;

MbDCode/DispCompIecJecKeqc.cpp

@@ -7,7 +7,7 @@ MbD::DispCompIecJecKeqc::DispCompIecJecKeqc()
 {
 }
 
-MbD::DispCompIecJecKeqc::DispCompIecJecKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk) : DispCompIecJecKec(frmi, frmj, frmk, axisk)
+MbD::DispCompIecJecKeqc::DispCompIecJecKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk) : DispCompIecJecKec(frmi, frmj, frmk, axisk)
 {
 }
 

MbDCode/DispCompIecJecKeqc.h

@@ -8,7 +8,7 @@ namespace MbD {
         //priIeJeKepEK ppriIeJeKepEKpEK pAjOKepEKT ppAjOKepEKpEK 
     public:
         DispCompIecJecKeqc();
-        DispCompIecJecKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk);
+        DispCompIecJecKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk);
         void initialize() override;
         void initializeGlobally() override;
         void calcPostDynCorrectorIteration() override;

MbDCode/DispCompIecJecO.cpp

@@ -1,10 +1,12 @@
 #include "DispCompIecJecO.h"
 
+using namespace MbD;
+
 MbD::DispCompIecJecO::DispCompIecJecO()
 {
 }
 
-MbD::DispCompIecJecO::DispCompIecJecO(EndFrmcptr frmi, EndFrmcptr frmj, int axis) : KinematicIeJe(frmi, frmj), axis(axis)
+MbD::DispCompIecJecO::DispCompIecJecO(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis) : KinematicIeJe(frmi, frmj), axis(axis)
 {
 }
 

MbDCode/DispCompIecJecO.h

@@ -8,10 +8,10 @@ namespace MbD {
 		//axis riIeJeO 
 	public:
 		DispCompIecJecO();
-		DispCompIecJecO(EndFrmcptr frmi, EndFrmcptr frmj, int axis);
+		DispCompIecJecO(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis);
 		void calcPostDynCorrectorIteration() override;
 
-		int axis;
+		size_t axis = -1;
 		double riIeJeO;
 	};
 }

MbDCode/DispCompIeqcJecKeqc.cpp

@@ -7,7 +7,7 @@ MbD::DispCompIeqcJecKeqc::DispCompIeqcJecKeqc()
 {
 }
 
-MbD::DispCompIeqcJecKeqc::DispCompIeqcJecKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk) : DispCompIecJecKeqc(frmi, frmj, frmk, axisk)
+MbD::DispCompIeqcJecKeqc::DispCompIeqcJecKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk) : DispCompIecJecKeqc(frmi, frmj, frmk, axisk)
 {
 }
 

MbDCode/DispCompIeqcJecKeqc.h

@@ -8,7 +8,7 @@ namespace MbD {
         //priIeJeKepXI priIeJeKepEI ppriIeJeKepXIpEK ppriIeJeKepEIpEI ppriIeJeKepEIpEK 
     public:
         DispCompIeqcJecKeqc();
-        DispCompIeqcJecKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk);
+        DispCompIeqcJecKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk);
         void initialize() override;
         void calcPostDynCorrectorIteration() override;
 

MbDCode/DispCompIeqcJecO.cpp

@@ -1,11 +1,13 @@
 #include "DispCompIeqcJecO.h"
 #include "EndFrameqc.h"
 
+using namespace MbD;
+
 MbD::DispCompIeqcJecO::DispCompIeqcJecO()
 {
 }
 
-MbD::DispCompIeqcJecO::DispCompIeqcJecO(EndFrmcptr frmi, EndFrmcptr frmj, int axis) : DispCompIecJecO(frmi, frmj, axis)
+MbD::DispCompIeqcJecO::DispCompIeqcJecO(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis) : DispCompIecJecO(frmi, frmj, axis)
 {
 }
 

MbDCode/DispCompIeqcJecO.h

@@ -8,7 +8,7 @@ namespace MbD {
         //priIeJeOpXI priIeJeOpEI ppriIeJeOpEIpEI 
     public:
         DispCompIeqcJecO();
-        DispCompIeqcJecO(EndFrmcptr frmi, EndFrmcptr frmj, int axis);
+        DispCompIeqcJecO(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis);
         void initializeGlobally() override;
         void calcPostDynCorrectorIteration() override;
 

MbDCode/DispCompIeqcJeqcKeqc.cpp

@@ -1,11 +1,13 @@
 #include "DispCompIeqcJeqcKeqc.h"
 #include "EndFrameqc.h"
 
+using namespace MbD;
+
 MbD::DispCompIeqcJeqcKeqc::DispCompIeqcJeqcKeqc()
 {
 }
 
-MbD::DispCompIeqcJeqcKeqc::DispCompIeqcJeqcKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk) : DispCompIeqcJecKeqc(frmi, frmj, frmk, axisk)
+MbD::DispCompIeqcJeqcKeqc::DispCompIeqcJeqcKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk) : DispCompIeqcJecKeqc(frmi, frmj, frmk, axisk)
 {
 }
 

MbDCode/DispCompIeqcJeqcKeqc.h

@@ -8,7 +8,7 @@ namespace MbD {
         //priIeJeKepXJ priIeJeKepEJ ppriIeJeKepXJpEK ppriIeJeKepEJpEJ ppriIeJeKepEJpEK 
     public:
         DispCompIeqcJeqcKeqc();
-        DispCompIeqcJeqcKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk);
+        DispCompIeqcJeqcKeqc(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk);
         void initialize() override;
         void calcPostDynCorrectorIteration() override;
 

MbDCode/DispCompIeqcJeqcKeqct.cpp

@@ -1,11 +1,13 @@
 #include "DispCompIeqcJeqcKeqct.h"
 #include "EndFrameqc.h"
 
+using namespace MbD;
+
 MbD::DispCompIeqcJeqcKeqct::DispCompIeqcJeqcKeqct()
 {
 }
 
-MbD::DispCompIeqcJeqcKeqct::DispCompIeqcJeqcKeqct(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk) : DispCompIeqcJeqcKeqc(frmi, frmj, frmk, axisk)
+MbD::DispCompIeqcJeqcKeqct::DispCompIeqcJeqcKeqct(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk) : DispCompIeqcJeqcKeqc(frmi, frmj, frmk, axisk)
 {
 }
 

MbDCode/DispCompIeqcJeqcKeqct.h

@@ -8,7 +8,7 @@ namespace MbD {
         //priIeJeKept ppriIeJeKepXIpt ppriIeJeKepEIpt ppriIeJeKepXJpt ppriIeJeKepEJpt ppriIeJeKepEKpt ppriIeJeKeptpt 
     public:
         DispCompIeqcJeqcKeqct();
-        DispCompIeqcJeqcKeqct(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk);
+        DispCompIeqcJeqcKeqct(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk);
         void initialize() override;
         void calcPostDynCorrectorIteration() override;
 

MbDCode/DispCompIeqcJeqcO.cpp

@@ -1,11 +1,13 @@
 #include "DispCompIeqcJeqcO.h"
 #include "EndFrameqc.h"
 
+using namespace MbD;
+
 MbD::DispCompIeqcJeqcO::DispCompIeqcJeqcO()
 {
 }
 
-MbD::DispCompIeqcJeqcO::DispCompIeqcJeqcO(EndFrmcptr frmi, EndFrmcptr frmj, int axis) : DispCompIeqcJecO(frmi, frmj, axis)
+MbD::DispCompIeqcJeqcO::DispCompIeqcJeqcO(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis) : DispCompIeqcJecO(frmi, frmj, axis)
 {
 }
 

MbDCode/DispCompIeqcJeqcO.h

@@ -8,7 +8,7 @@ namespace MbD {
         //priIeJeOpXJ priIeJeOpEJ ppriIeJeOpEJpEJ 
     public:
         DispCompIeqcJeqcO();
-        DispCompIeqcJeqcO(EndFrmcptr frmi, EndFrmcptr frmj, int axis);
+        DispCompIeqcJeqcO(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis);
         void initializeGlobally() override;
         void calcPostDynCorrectorIteration() override;
 

MbDCode/DispCompIeqctJeqcKeqct.cpp

@@ -1,9 +1,11 @@
 #include "DispCompIeqctJeqcKeqct.h"
 
+using namespace MbD;
+
 MbD::DispCompIeqctJeqcKeqct::DispCompIeqctJeqcKeqct()
 {
 }
 
-MbD::DispCompIeqctJeqcKeqct::DispCompIeqctJeqcKeqct(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk) : DispCompIeqcJeqcKeqct(frmi, frmj, frmk, axisk)
+MbD::DispCompIeqctJeqcKeqct::DispCompIeqctJeqcKeqct(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk) : DispCompIeqcJeqcKeqct(frmi, frmj, frmk, axisk)
 {
 }

MbDCode/DispCompIeqctJeqcKeqct.h

@@ -8,7 +8,7 @@ namespace MbD {
         //
     public:
         DispCompIeqctJeqcKeqct();
-        DispCompIeqctJeqcKeqct(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, int axisk);
+        DispCompIeqctJeqcKeqct(EndFrmcptr frmi, EndFrmcptr frmj, EndFrmcptr frmk, size_t axisk);
 
     };
 }

MbDCode/DispCompIeqctJeqcO.cpp

@@ -1,11 +1,13 @@
 #include "DispCompIeqctJeqcO.h"
 #include "EndFrameqct.h"
 
+using namespace MbD;
+
 MbD::DispCompIeqctJeqcO::DispCompIeqctJeqcO()
 {
 }
 
-MbD::DispCompIeqctJeqcO::DispCompIeqctJeqcO(EndFrmcptr frmi, EndFrmcptr frmj, int axis) : DispCompIeqcJeqcO(frmi, frmj, axis)
+MbD::DispCompIeqctJeqcO::DispCompIeqctJeqcO(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis) : DispCompIeqcJeqcO(frmi, frmj, axis)
 {
 }
 

MbDCode/DispCompIeqctJeqcO.h

@@ -8,7 +8,7 @@ namespace MbD {
         //priIeJeOpt ppriIeJeOpEIpt ppriIeJeOptpt 
     public:
         DispCompIeqctJeqcO();
-        DispCompIeqctJeqcO(EndFrmcptr frmi, EndFrmcptr frmj, int axis);
+        DispCompIeqctJeqcO(EndFrmcptr frmi, EndFrmcptr frmj, size_t axis);
         void initializeGlobally() override;
         void calcPostDynCorrectorIteration() override;
 

MbDCode/EndFramec.cpp

@@ -45,12 +45,12 @@ void MbD::EndFramec::calcPostDynCorrectorIteration()
 	aAOe = markerFrame->aAOm;
 }
 
-FColDsptr MbD::EndFramec::aAjOe(int j)
+FColDsptr MbD::EndFramec::aAjOe(size_t j)
 {
 	return aAOe->column(j);
 }
 
-double MbD::EndFramec::riOeO(int i)
+double MbD::EndFramec::riOeO(size_t i)
 {
 	return rOeO->at(i);
 }

MbDCode/EndFramec.h

@@ -26,10 +26,10 @@ namespace MbD {
         void initializeGlobally() override;
         virtual void initEndFrameqct();
         void calcPostDynCorrectorIteration() override;
-        FColDsptr aAjOe(int j);
-        double riOeO(int i);
+        FColDsptr aAjOe(size_t j);
+        double riOeO(size_t i);
 
-        MarkerFrame* markerFrame;
+        MarkerFrame* markerFrame; //Use raw pointer when pointing backwards.
         FColDsptr rOeO = std::make_shared<FullColumn<double>>(3);
         FMatDsptr aAOe = std::make_shared<FullMatrix<double>>(3, 3);
     };

MbDCode/EndFrameqc.cpp

@@ -38,7 +38,7 @@ void EndFrameqc::initEndFrameqct()
 	endFrameqct->setMarkerFrame(markerFrame);
 }
 
-FMatFColDsptr MbD::EndFrameqc::ppAjOepEpE(int jj)
+FMatFColDsptr MbD::EndFrameqc::ppAjOepEpE(size_t jj)
 {
 	auto answer = std::make_shared<FullMatrix<std::shared_ptr<FullColumn<double>>>>(4, 4);
 	for (size_t i = 0; i < 4; i++) {
@@ -59,7 +59,7 @@ void MbD::EndFrameqc::calcPostDynCorrectorIteration()
 	pAOepE = markerFrame->pAOmpE;
 }
 
-FMatDsptr MbD::EndFrameqc::pAjOepET(int axis)
+FMatDsptr MbD::EndFrameqc::pAjOepET(size_t axis)
 {
 	auto answer = std::make_shared<FullMatrix<double>>(4, 3);
 	for (size_t i = 0; i < 4; i++) {
@@ -73,7 +73,7 @@ FMatDsptr MbD::EndFrameqc::pAjOepET(int axis)
 	return answer;
 }
 
-FMatDsptr MbD::EndFrameqc::ppriOeOpEpE(int ii)
+FMatDsptr MbD::EndFrameqc::ppriOeOpEpE(size_t ii)
 {
 	auto answer = std::make_shared<FullMatrix<double>>(4, 4);
 	for (size_t i = 0; i < 4; i++) {
@@ -87,17 +87,17 @@ FMatDsptr MbD::EndFrameqc::ppriOeOpEpE(int ii)
 	return answer;
 }
 
-int MbD::EndFrameqc::iqX()
+size_t MbD::EndFrameqc::iqX()
 {
 	return markerFrame->iqX();
 }
 
-int MbD::EndFrameqc::iqE()
+size_t MbD::EndFrameqc::iqE()
 {
 	return markerFrame->iqE();
 }
 
-FRowDsptr MbD::EndFrameqc::priOeOpE(int i)
+FRowDsptr MbD::EndFrameqc::priOeOpE(size_t i)
 {
 	return prOeOpE->at(i);
 }

MbDCode/EndFrameqc.h

@@ -15,13 +15,13 @@ namespace MbD {
         void initialize();
         void initializeGlobally() override;
         void initEndFrameqct() override;
-        FMatFColDsptr ppAjOepEpE(int j);
+        FMatFColDsptr ppAjOepEpE(size_t j);
         void calcPostDynCorrectorIteration() override;
-        FMatDsptr pAjOepET(int j);
-        FMatDsptr ppriOeOpEpE(int i);
-        int iqX();
-        int iqE();
-        FRowDsptr priOeOpE(int i);
+        FMatDsptr pAjOepET(size_t j);
+        FMatDsptr ppriOeOpEpE(size_t i);
+        size_t iqX();
+        size_t iqE();
+        FRowDsptr priOeOpE(size_t i);
 
         FMatDsptr prOeOpE;
         std::shared_ptr<FullMatrix<std::shared_ptr<FullColumn<double>>>> pprOeOpEpE;

MbDCode/EulerAngleszxz.cpp

@@ -1 +1,3 @@
 #include "EulerAngleszxz.h"
+
+using namespace MbD;

MbDCode/EulerArray.cpp

@@ -1 +1,3 @@
 #include "EulerArray.h"
+
+using namespace MbD;

MbDCode/EulerConstraint.cpp

@@ -26,4 +26,9 @@ void MbD::EulerConstraint::calcPostDynCorrectorIteration()
 	{
 		pGpE->at(i) = 2.0 * qE->at(i);
 	}
-}
+}
+
+void MbD::EulerConstraint::useEquationNumbers()
+{
+	iqE = static_cast<PartFrame*>(owner)->iqE;
+}

MbDCode/EulerConstraint.h

@@ -14,9 +14,10 @@ namespace MbD {
 		EulerConstraint(const char* str);
 		void initialize();
 		void calcPostDynCorrectorIteration() override;
+		void useEquationNumbers() override;
 
 		FRowDsptr pGpE;		//partial derivative of G wrt pE
-		int iqE = -1;
+		size_t iqE = -1;
 	};
 }
 

MbDCode/EulerParameters.h

@@ -14,7 +14,6 @@ namespace MbD {
 		EulerParameters(size_t count) : EulerArray<T>(count) {}
 		EulerParameters(size_t count, const T& value) : EulerArray<T>(count, value) {}
 		EulerParameters(std::initializer_list<T> list) : EulerArray<T>{ list } {}
-		double sumOfSquares();
 
 		static std::shared_ptr<FullMatrix<std::shared_ptr<FullColumn<T>>>> ppApEpEtimesColumn(FColDsptr col);
 		static std::shared_ptr<FullMatrix<std::shared_ptr<FullMatrix<T>>>> ppApEpEtimesMatrix(FMatDsptr mat);
@@ -29,12 +28,6 @@ namespace MbD {
 		FColFMatDsptr pApE;
 	};
 
-	template<typename T>
-	inline double EulerParameters<T>::sumOfSquares()
-	{
-		return 0.0;
-	}
-
 	template<>
 	inline std::shared_ptr<FullMatrix<std::shared_ptr<FullColumn<double>>>> EulerParameters<double>::ppApEpEtimesColumn(FColDsptr col)
 	{

MbDCode/EulerParametersDot.cpp

@@ -1 +1,3 @@
 #include "EulerParametersDot.h"
+
+using namespace MbD;

MbDCode/ExpressionX.cpp

@@ -1 +1,3 @@
 #include "ExpressionX.h"
+
+using namespace MbD;

MbDCode/ForceTorqueItem.cpp

@@ -1 +1,3 @@
 #include "ForceTorqueItem.h"
+
+using namespace MbD;

MbDCode/FullColumn.h

@@ -2,21 +2,28 @@
 #include <string>
 #include <sstream> 
 
-#include "Vector.h"
+#include "FullVector.h"
 
 namespace MbD {
 
 	template <typename T>
-	class FullColumn : public Vector<T>
+	class FullColumn : public FullVector<T>
 	{
 	public:
-		FullColumn(size_t count) : Vector<T>(count) {}
-		FullColumn(size_t count, const T& value) : Vector<T>(count, value) {}
-		FullColumn(std::initializer_list<T> list) : Vector<T>{ list } {}
+		FullColumn(size_t count) : FullVector<T>(count) {}
+		FullColumn(size_t count, const T& value) : FullVector<T>(count, value) {}
+		FullColumn(std::initializer_list<T> list) : FullVector<T>{ list } {}
 		std::shared_ptr<FullColumn<T>> plusFullColumn(std::shared_ptr<FullColumn<T>> fullCol);
 		std::shared_ptr<FullColumn<T>> minusFullColumn(std::shared_ptr<FullColumn<T>> fullCol);
 		std::shared_ptr<FullColumn<T>> times(double a);
 		std::shared_ptr<FullColumn<T>> negated();
+		void atiputFullColumn(size_t i, std::shared_ptr<FullColumn<T>> fullCol);
+		void equalSelfPlusFullColumnAt(std::shared_ptr<FullColumn<T>> fullCol, size_t i);
+		void atiminusFullColumn(size_t i, std::shared_ptr<FullColumn<T>> fullCol);
+		void equalFullColumnAt(std::shared_ptr<FullColumn<T>> fullCol, size_t i);
+		std::shared_ptr<FullColumn<T>> copy();
+		
+
 		std::string toString()
 		{
 			std::stringstream ss;
@@ -65,5 +72,50 @@ namespace MbD {
 	{
 		return this->times(-1.0);
 	}
+	template<typename T>
+	inline void FullColumn<T>::atiputFullColumn(size_t i, std::shared_ptr<FullColumn<T>> fullCol)
+	{
+		for (size_t ii = 0; ii < fullCol->size(); ii++)
+		{
+			this->at(i + ii) = fullCol->at(ii);
+		}
+	}
+	template<typename T>
+	inline void FullColumn<T>::equalSelfPlusFullColumnAt(std::shared_ptr<FullColumn<T>> fullCol, size_t ii)
+	{
+		//self is subcolumn of fullCol
+		for (size_t i = 0; i < this->size(); i++)
+		{
+			this->at(i) += fullCol->at(ii + i);
+		}
+	}
+	template<typename T>
+	inline void FullColumn<T>::atiminusFullColumn(size_t i1, std::shared_ptr<FullColumn<T>> fullCol)
+	{
+		for (size_t ii = 0; ii < fullCol->size(); ii++)
+		{
+			size_t i = i1 + ii;
+			this->at(i) -= fullCol->at(ii);
+		}
+	}
+	template<typename T>
+	inline void FullColumn<T>::equalFullColumnAt(std::shared_ptr<FullColumn<T>> fullCol, size_t ii)
+	{
+		for (size_t i = 0; i < this->size(); i++)
+		{
+			this->at(i) = fullCol->at(ii + i);
+		}
+	}
+	template<>
+	inline std::shared_ptr<FullColumn<double>> FullColumn<double>::copy()
+	{
+		auto n = this->size();
+		auto answer = std::make_shared<FullColumn<double>>(n);
+		for (size_t i = 0; i < n; i++)
+		{
+			answer->at(i) = this->at(i);
+		}
+		return answer;
+	}
 }
 

MbDCode/FullMatrix.h

@@ -14,13 +14,10 @@ namespace MbD {
 	{
 	public:
 		FullMatrix() {}
-		FullMatrix(int m) {
-			for (size_t i = 0; i < m; i++) {
-				auto row = std::make_shared<FullRow<T>>();
-				this->push_back(row);
-			}
+		FullMatrix(size_t m) : RowTypeMatrix<std::shared_ptr<FullRow<T>>>(m) 
+		{
 		}
-		FullMatrix(int m, int n) {
+		FullMatrix(size_t m, size_t n) {
 			for (size_t i = 0; i < m; i++) {
 				auto row = std::make_shared<FullRow<T>>(n);
 				this->push_back(row);
@@ -40,7 +37,7 @@ namespace MbD {
 			}
 		}
 		void identity();
-		std::shared_ptr<FullColumn<T>> column(int j);
+		std::shared_ptr<FullColumn<T>> column(size_t j);
 		std::shared_ptr<FullColumn<T>> timesFullColumn(std::shared_ptr<FullColumn<T>> fullCol);
 		std::shared_ptr<FullMatrix<T>> timesFullMatrix(std::shared_ptr<FullMatrix<T>> fullMat);
 		std::shared_ptr<FullMatrix<T>> timesTransposeFullMatrix(std::shared_ptr<FullMatrix<T>> fullMat);
@@ -50,6 +47,9 @@ namespace MbD {
 		std::shared_ptr<FullMatrix<T>> negated();
 		void symLowerWithUpper();
 		void atijputFullColumn(size_t i, size_t j, std::shared_ptr<FullColumn<T>> fullCol);
+		double sumOfSquares() override;
+		void atijplusNumber(size_t i, size_t j, double value);
+		void zeroSelf() override;
 	};
 	template <>
 	inline void FullMatrix<double>::identity() {
@@ -59,7 +59,7 @@ namespace MbD {
 		}
 	}
 	template<typename T>
-	inline std::shared_ptr<FullColumn<T>> FullMatrix<T>::column(int j) {
+	inline std::shared_ptr<FullColumn<T>> FullMatrix<T>::column(size_t j) {
 		size_t n = this->size();
 		auto answer = std::make_shared<FullColumn<T>>(n);
 		for (size_t i = 0; i < n; i++) {
@@ -82,7 +82,7 @@ namespace MbD {
 	template<typename T>
 	inline std::shared_ptr<FullMatrix<T>> FullMatrix<T>::timesFullMatrix(std::shared_ptr<FullMatrix<T>> fullMat)
 	{
-		size_t m = this->nRow();
+		size_t m = this->nrow();
 		auto answer = std::make_shared<FullMatrix<T>>(m);
 		for (size_t i = 0; i < m; i++) {
 			answer->at(i) = this->at(i)->timesFullMatrix(fullMat);
@@ -92,7 +92,7 @@ namespace MbD {
 	template<typename T>
 	inline std::shared_ptr<FullMatrix<T>> FullMatrix<T>::timesTransposeFullMatrix(std::shared_ptr<FullMatrix<T>> fullMat)
 	{
-		size_t nrow = this->nRow();
+		size_t nrow = this->nrow();
 		auto answer = std::make_shared<FullMatrix<T>>(nrow);
 		for (size_t i = 0; i < nrow; i++) {
 			answer->at(i) = this->at(i)->timesTransposeFullMatrix(fullMat);
@@ -102,7 +102,7 @@ namespace MbD {
 	template<typename T>
 	inline std::shared_ptr<FullMatrix<T>> FullMatrix<T>::times(double a)
 	{
-		size_t m = this->nRow();
+		size_t m = this->nrow();
 		auto answer = std::make_shared<FullMatrix<T>>(m);
 		for (size_t i = 0; i < m; i++) {
 			answer->at(i) = this->at(i)->times(a);
@@ -122,8 +122,8 @@ namespace MbD {
 	template<typename T>
 	inline std::shared_ptr<FullMatrix<T>> FullMatrix<T>::transpose()
 	{
-		size_t nrow = this->nRow();
-		auto ncol = this->nCol();
+		size_t nrow = this->nrow();
+		auto ncol = this->ncol();
 		auto answer = std::make_shared<FullMatrix<T>>(ncol, nrow);
 		for (size_t i = 0; i < nrow; i++) {
 			auto& row = this->at(i);
@@ -143,7 +143,7 @@ namespace MbD {
 	{
 		size_t n = this->size();
 		for (size_t i = 0; i < n; i++) {
-			for (size_t j = i+1; j < n; j++) {
+			for (size_t j = i + 1; j < n; j++) {
 				this->at(j)->at(i) = this->at(i)->at(j);
 			}
 		}
@@ -157,6 +157,44 @@ namespace MbD {
 			this->at(iOffset + ii)->at(j1) = fullCol->at(ii);
 		}
 	}
+	template<>
+	inline double FullMatrix<double>::sumOfSquares()
+	{
+		double sum = 0.0;
+		for (size_t i = 0; i < this->size(); i++)
+		{
+			sum += this->at(i)->sumOfSquares();
+		}
+		return sum;
+	}
+	template<typename T>
+	inline double FullMatrix<T>::sumOfSquares()
+	{
+		assert(false);
+		return 0.0;
+	}
+	template<>
+	inline void FullMatrix<double>::atijplusNumber(size_t i, size_t j, double value)
+	{
+		this->at(i)->atiplusNumber(j, value);
+	}
+	template<typename T>
+	inline void FullMatrix<T>::atijplusNumber(size_t i, size_t j, double value)
+	{
+		assert(false);
+	}
+	template<>
+	inline void FullMatrix<double>::zeroSelf()
+	{
+		for (size_t i = 0; i < this->size(); i++) {
+			this->at(i)->zeroSelf();
+		}
+	}
+	template<typename T>
+	inline void FullMatrix<T>::zeroSelf()
+	{
+		assert(false);
+	}
 	using FMatDsptr = std::shared_ptr<FullMatrix<double>>;
 	using FMatDsptr = std::shared_ptr<FullMatrix<double>>;
 	using FMatFColDsptr = std::shared_ptr<FullMatrix<std::shared_ptr<FullColumn<double>>>>;

MbDCode/FullRow.cpp

@@ -1 +1,3 @@
 #include "FullRow.h"
+
+using namespace MbD;

MbDCode/FullRow.h

@@ -1,5 +1,5 @@
 #pragma once
-#include "Vector.h"
+#include "FullVector.h"
 #include "FullColumn.h"
 
 namespace MbD {
@@ -7,13 +7,13 @@ namespace MbD {
 	class FullMatrix;
 
 	template <typename T>
-	class FullRow : public Vector<T>
+	class FullRow : public FullVector<T>
 	{
 	public:
 		FullRow() {}
-		FullRow(size_t count) : Vector<T>(count) {}
-		FullRow(size_t count, const T& value) : Vector<T>(count, value) {}
-		FullRow(std::initializer_list<T> list) : Vector<T>{ list } {}
+		FullRow(size_t count) : FullVector<T>(count) {}
+		FullRow(size_t count, const T& value) : FullVector<T>(count, value) {}
+		FullRow(std::initializer_list<T> list) : FullVector<T>{ list } {}
 		std::shared_ptr<FullRow<T>> times(double a);
 		std::shared_ptr<FullRow<T>> negated();
 		std::shared_ptr<FullRow<T>> plusFullRow(std::shared_ptr<FullRow<T>> fullRow);
@@ -73,7 +73,7 @@ namespace MbD {
 	inline std::shared_ptr<FullRow<T>> FullRow<T>::timesTransposeFullMatrix(std::shared_ptr<FullMatrix<T>> fullMat)
 	{
 		//"a*bT = a(1,j)b(k,j)"
-		size_t ncol = fullMat->nRow();
+		size_t ncol = fullMat->nrow();
 		auto answer = std::make_shared<FullRow<T>>(ncol);
 		for (size_t k = 0; k < ncol; k++) {
 			answer->at(k) = this->dot(fullMat->at(k));

MbDCode/FullVector.cpp

@@ -1,4 +1,4 @@
-#include "Vector.h"
+#include "FullVector.h"
 
 using namespace MbD;
 

MbDCode/FullVector.h

@@ -0,0 +1,74 @@
+#pragma once
+#include "Array.h"
+
+namespace MbD {
+    template <typename T>
+    class FullVector : public Array<T>
+    {
+    public:
+        FullVector() {}
+        FullVector(size_t count) : Array<T>(count) {}
+        FullVector(size_t count, const T& value) : Array<T>(count, value) {}
+        FullVector(std::initializer_list<T> list) : Array<T>{ list } {}
+        double dot(std::shared_ptr<FullVector<T>> vec);
+        void atiplusNumber(size_t i, T value);
+        double sumOfSquares() override;
+        size_t numberOfElements() override;
+        void zeroSelf() override;
+        void atitimes(size_t i, double factor);
+    };
+    template<typename T>
+    inline double FullVector<T>::dot(std::shared_ptr<FullVector<T>> vec)
+    {
+        size_t n = this->size();
+        double answer = 0.0;
+        for (size_t i = 0; i < n; i++) {
+            answer += this->at(i) * vec->at(i);
+        }
+        return answer;
+    }
+    template<typename T>
+    inline void FullVector<T>::atiplusNumber(size_t i, T value)
+    {
+        this->at(i) += value;
+    }
+    template<>
+    inline double FullVector<double>::sumOfSquares()
+    {
+        double sum = 0.0;
+        for (size_t i = 0; i < this->size(); i++)
+        {
+            double element = this->at(i);
+            sum += element * element;
+        }
+        return sum;
+    }
+    template<typename T>
+    inline double FullVector<T>::sumOfSquares()
+    {
+        assert(false);
+        return 0.0;
+    }
+    template<typename T>
+    inline size_t FullVector<T>::numberOfElements()
+    {
+        return this->size();
+    }
+    template<>
+    inline void FullVector<double>::zeroSelf()
+    {
+        for (size_t i = 0; i < this->size(); i++) {
+            this->at(i) = 0.0;;
+        }
+    }
+    template<typename T>
+    inline void FullVector<T>::zeroSelf()
+    {
+        assert(false);
+    }
+    template<typename T>
+    inline void FullVector<T>::atitimes(size_t i, double factor)
+    {
+        this->at(i) *= factor;
+    }
+}

MbDCode/Function.cpp

@@ -1 +1,3 @@
 #include "Function.h"
+
+using namespace MbD;

MbDCode/FunctionX.cpp

@@ -1 +1,3 @@
 #include "FunctionX.h"
+
+using namespace MbD;

MbDCode/GEFullMat.cpp

@@ -1 +1,25 @@
+#include <cassert>
+
 #include "GEFullMat.h"
+
+using namespace MbD;
+
+void MbD::GEFullMat::forwardEliminateWithPivot(size_t p)
+{
+	assert(false);
+}
+
+void MbD::GEFullMat::backSubstituteIntoDU()
+{
+	assert(false);
+}
+
+void MbD::GEFullMat::postSolve()
+{
+	assert(false);
+}
+
+void MbD::GEFullMat::preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal)
+{
+	assert(false);
+}

MbDCode/GEFullMat.h

@@ -7,7 +7,11 @@ namespace MbD {
     {
         //
     public:
-
+        std::shared_ptr<FullMatrix<double>> matrixA;
+        void forwardEliminateWithPivot(size_t p) override;
+        void backSubstituteIntoDU() override;
+        void postSolve() override;
+        void preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal) override;
     };
 }
 

MbDCode/GEFullMatFullPv.cpp

@@ -1 +1,15 @@
+#include <cassert>
+
 #include "GEFullMatFullPv.h"
+
+using namespace MbD;
+
+void MbD::GEFullMatFullPv::doPivoting(size_t p)
+{
+	assert(false);
+}
+
+void MbD::GEFullMatFullPv::postSolve()
+{
+	assert(false);
+}

MbDCode/GEFullMatFullPv.h

@@ -7,7 +7,8 @@ namespace MbD {
     {
         //
     public:
-
+        void doPivoting(size_t p) override;
+        void postSolve() override;
     };
 }
 

MbDCode/GEFullMatParPv.cpp

@@ -1 +1,20 @@
+#include <cassert>
+
 #include "GEFullMatParPv.h"
+
+using namespace MbD;
+
+void MbD::GEFullMatParPv::doPivoting(size_t p)
+{
+	assert(false);
+}
+
+void MbD::GEFullMatParPv::postSolve()
+{
+	assert(false);
+}
+
+void MbD::GEFullMatParPv::preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal)
+{
+	assert(false);
+}

MbDCode/GEFullMatParPv.h

@@ -7,7 +7,9 @@ namespace MbD {
     {
         //
     public:
-
+        void doPivoting(size_t p) override;
+        void postSolve() override;
+        void preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal) override;
     };
 }
 

MbDCode/GESpMat3.cpp

@@ -1 +1,11 @@
 #include "GESpMat3.h"
+#include "FullColumn.h"
+#include "SparseMatrix.h"
+
+using namespace MbD;
+
+FColDsptr MbD::GESpMat::solvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal)
+{
+	this->timedSolvewithsaveOriginal(spMat, fullCol, saveOriginal);
+	return answerX;
+}

MbDCode/GESpMat3.h

@@ -1,13 +1,18 @@
 #pragma once
 
 #include "MatrixGaussElimination.h"
+#include "SparseMatrix.h"
 
 namespace MbD {
     class GESpMat : public MatrixGaussElimination
     {
-        //
+        //markowitzPivotRowCount markowitzPivotColCount privateIndicesOfNonZerosInPivotRow rowPositionsOfNonZerosInPivotColumn 
     public:
+        FColDsptr solvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal);
 
+        std::shared_ptr<SparseMatrix<double>> matrixA;
+        size_t markowitzPivotRowCount, markowitzPivotColCount;
+        std::shared_ptr<std::vector<size_t>> privateIndicesOfNonZerosInPivotRow, rowPositionsOfNonZerosInPivotColumn;
     };
 }
 

MbDCode/GESpMatFullPv.cpp

@@ -1 +1,30 @@
+#include <cassert>
+
 #include "GESpMatFullPv.h"
+
+using namespace MbD;
+
+void MbD::GESpMatFullPv::doPivoting(size_t p)
+{
+	assert(false);
+}
+
+void MbD::GESpMatFullPv::forwardEliminateWithPivot(size_t p)
+{
+	assert(false);
+}
+
+void MbD::GESpMatFullPv::backSubstituteIntoDU()
+{
+	assert(false);
+}
+
+void MbD::GESpMatFullPv::postSolve()
+{
+	assert(false);
+}
+
+void MbD::GESpMatFullPv::preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal)
+{
+	assert(false);
+}

MbDCode/GESpMatFullPv.h

@@ -5,9 +5,16 @@
 namespace MbD {
     class GESpMatFullPv : public GESpMat
     {
-        //
+        //positionsOfOriginalCols rowPositionsOfNonZerosInColumns 
     public:
+        void doPivoting(size_t p) override;
+        void forwardEliminateWithPivot(size_t p) override;
+        void backSubstituteIntoDU() override;
+        void postSolve() override;
+        void preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal) override;
 
+        std::shared_ptr<std::vector<size_t>> positionsOfOriginalCols;
+        std::shared_ptr<std::vector<std::shared_ptr<std::vector<size_t>>>> rowPositionsOfNonZerosInColumns;
     };
 }
 

MbDCode/GESpMatFullPvPosIC.cpp

@@ -1 +1,46 @@
+#include <algorithm>
+#include <cassert>
+
 #include "GESpMatFullPvPosIC.h"
+//#include "FullColumn.h"
+//#include "SparseMatrix.h"
+#include "PosICNewtonRaphson.h"
+
+using namespace MbD;
+
+void MbD::GESpMatFullPvPosIC::preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal)
+{
+	GESpMatFullPv::preSolvewithsaveOriginal(spMat, fullCol, saveOriginal);
+	if (system == nullptr) {
+		pivotRowLimits = std::make_shared<std::vector<size_t>>();
+	}
+	else {
+		pivotRowLimits = system->pivotRowLimits;
+	}
+	pivotRowLimit = -1;
+}
+
+void MbD::GESpMatFullPvPosIC::doPivoting(size_t p)
+{
+	//"Used by Gauss Elimination only."
+	//"Swap rows but keep columns in place."
+	//"The elements below the diagonal are removed column by column."
+
+	assert(false);
+	auto max = 0.0;
+	auto pivotRow = p;
+	auto pivotCol = p;
+	for (size_t j = p; j < n; j++)
+	{
+		rowPositionsOfNonZerosInColumns->at(colOrder->at(j))->clear();
+	}
+	if (p > pivotRowLimit) {
+		pivotRowLimit = *std::find_if(
+			pivotRowLimits->begin(), pivotRowLimits->end(),
+			[&](size_t limit) { return limit >= p; });
+	}
+	for (size_t i = p; i < pivotRowLimit; i++)
+	{
+		auto& rowi = matrixA->at(i);
+	}
+}

MbDCode/GESpMatFullPvPosIC.h

@@ -3,11 +3,18 @@
 #include "GESpMatFullPv.h"
 
 namespace MbD {
+    class PosICNewtonRaphson;
+
     class GESpMatFullPvPosIC : public GESpMatFullPv
     {
-        //
+        //system pivotRowLimits pivotRowLimit 
     public:
+        void preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal) override;
+        void doPivoting(size_t p) override;
 
+        PosICNewtonRaphson* system; //Use raw pointer when pointing backwards.
+        std::shared_ptr<std::vector<size_t>> pivotRowLimits;
+        size_t pivotRowLimit;
     };
 }
 

MbDCode/GESpMatParPv.cpp

@@ -1 +1,20 @@
+#include <cassert>
+
 #include "GESpMatParPv.h"
+
+using namespace MbD;
+
+void MbD::GESpMatParPv::forwardEliminateWithPivot(size_t p)
+{
+	assert(false);
+}
+
+void MbD::GESpMatParPv::backSubstituteIntoDU()
+{
+	assert(false);
+}
+
+void MbD::GESpMatParPv::postSolve()
+{
+	assert(false);
+}

MbDCode/GESpMatParPv.h

@@ -7,6 +7,9 @@ namespace MbD {
     {
         //
     public:
+        void forwardEliminateWithPivot(size_t p) override;
+        void backSubstituteIntoDU() override;
+        void postSolve() override;
 
     };
 }

MbDCode/GESpMatParPvMarko.cpp

@@ -1 +1,15 @@
+#include <cassert>
+
 #include "GESpMatParPvMarko.h"
+
+using namespace MbD;
+
+void MbD::GESpMatParPvMarko::doPivoting(size_t p)
+{
+	assert(false);
+}
+
+void MbD::GESpMatParPvMarko::preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal)
+{
+	assert(false);
+}

MbDCode/GESpMatParPvMarko.h

@@ -7,6 +7,8 @@ namespace MbD {
     {
         //
     public:
+        void doPivoting(size_t p) override;
+        void preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal) override;
 
     };
 }

MbDCode/GESpMatParPvMarkoFast.cpp

@@ -1 +1,42 @@
+#include <cassert>
+#include <memory>
+
 #include "GESpMatParPvMarkoFast.h"
+#include "SingularMatrixError.h"
+
+using namespace MbD;
+
+void MbD::GESpMatParPvMarkoFast::preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal)
+{
+	//assert(false);
+	//"Optimized for speed."
+	if (m != spMat->nrow() || n != spMat->ncol()) {
+		m = spMat->nrow();
+		n = spMat->ncol();
+		matrixA = std::make_shared<SparseMatrix<double>>(m);
+		privateIndicesOfNonZerosInPivotRow = std::make_shared<std::vector<size_t>>();
+		rowPositionsOfNonZerosInPivotColumn = std::make_shared<std::vector<size_t>>();
+	}
+	if (saveOriginal) {
+		rightHandSideB = fullCol->copy();
+	}
+	else {
+		rightHandSideB = fullCol;
+	}
+	for (size_t i = 0; i < m; i++)
+	{
+		auto& spRowi = spMat->at(i);
+		auto maxRowElement = spRowi->maxElement();
+		if (maxRowElement == 0) {
+			throw SingularMatrixError("");
+		}
+		auto scaling = 1.0 / maxRowElement;
+		matrixA->at(i) = spRowi->timesconditionedWithTol(scaling, singularPivotTolerance);
+		rightHandSideB->atitimes(i, scaling);
+	}
+}
+
+void MbD::GESpMatParPvMarkoFast::doPivoting(size_t p)
+{
+	assert(false);
+}

MbDCode/GESpMatParPvMarkoFast.h

@@ -7,6 +7,8 @@ namespace MbD {
     {
         //
     public:
+        void preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal) override;
+        void doPivoting(size_t p) override;
 
     };
 }

MbDCode/GESpMatParPvPrecise.cpp

@@ -1 +1,15 @@
+#include <cassert>
+
 #include "GESpMatParPvPrecise.h"
+
+using namespace MbD;
+
+void MbD::GESpMatParPvPrecise::doPivoting(size_t p)
+{
+	assert(false);
+}
+
+void MbD::GESpMatParPvPrecise::preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal)
+{
+	assert(false);
+}

MbDCode/GESpMatParPvPrecise.h

@@ -7,6 +7,8 @@ namespace MbD {
     {
         //
     public:
+        void doPivoting(size_t p) override;
+        void preSolvewithsaveOriginal(SpMatDsptr spMat, FColDsptr fullCol, bool saveOriginal) override;
 
     };
 }

MbDCode/Integrator.cpp

@@ -1 +1,3 @@
 #include "Integrator.h"
+
+using namespace MbD;