diff --git a/OndselSolver/FullColumn.cpp b/OndselSolver/FullColumn.cpp
index 75af3f7..9dab0db 100644
--- a/OndselSolver/FullColumn.cpp
+++ b/OndselSolver/FullColumn.cpp
@@ -164,16 +164,4 @@ namespace MbD {
         assert(false);
         return FColsptr<T>();
     }
-    template<typename T>
-    std::ostream& FullColumn<T>::printOn(std::ostream& s) const
-    {
-        s << "FullCol{";
-        s << this->at(0);
-        for (int i = 1; i < this->size(); i++)
-        {
-            s << ", " << this->at(i);
-        }
-        s << "}";
-        return s;
-    }
 }
\ No newline at end of file
diff --git a/OndselSolver/FullColumn.h b/OndselSolver/FullColumn.h
index 37e3d49..9bd13c4 100644
--- a/OndselSolver/FullColumn.h
+++ b/OndselSolver/FullColumn.h
@@ -15,6 +15,7 @@
 #include "FullColumn.ref.h"
 #include "FullRow.ref.h"
 #include "FullMatrix.ref.h"
+//#include "Symbolic.h"
 
 namespace MbD {
 	class Symbolic;
@@ -48,10 +49,21 @@ namespace MbD {
 
 		std::ostream& printOn(std::ostream& s) const override;
 	};
-    // instantiate on purpose:
+    // the following "printOn" needs to be in the header for unknown reasons linker
+    template<typename T>
+    std::ostream& FullColumn<T>::printOn(std::ostream& s) const
+    {
+        s << "FullCol{";
+        s << this->at(0);
+        for (int i = 1; i < this->size(); i++)
+        {
+            s << ", " << this->at(i);
+        }
+        s << "}";
+        return s;
+    }
+    // instantiate on purpose to make visible in library api:
     template class FullColumn<double>;
     template class FullColumn<int>;
-//    template class FullColumn<std::shared_ptr<MbD::FullMatrixDouble>>;
-//    template class FullColumn<std::shared_ptr<MbD::Symbolic>>;
 }
 
diff --git a/OndselSolver/FullMatrix.cpp b/OndselSolver/FullMatrix.cpp
index 6087540..1a12d03 100644
--- a/OndselSolver/FullMatrix.cpp
+++ b/OndselSolver/FullMatrix.cpp
@@ -313,15 +313,15 @@ namespace MbD {
         }
         return answer;
     }
-    std::shared_ptr<FullMatrixFullMatrixDouble> FullMatrixFullMatrixDouble::timesTransposeFullMatrix(std::shared_ptr<FullMatrixFullMatrixDouble> fullMat)
-    {
-        int nrow = this->nrow();
-        auto answer = std::make_shared<FullMatrixFullMatrixDouble>(nrow);
-        for (int i = 0; i < nrow; i++) {
-            answer->at(i) = this->at(i)->timesTransposeFullMatrixForFMFMDsptr(fullMat);
-        }
-        return answer;
-    }
+//    std::shared_ptr<FullMatrixFullMatrixDouble> FullMatrixFullMatrixDouble::timesTransposeFullMatrix(std::shared_ptr<FullMatrixFullMatrixDouble> fullMat)
+//    {
+//        int nrow = this->nrow();
+//        auto answer = std::make_shared<FullMatrixFullMatrixDouble>(nrow);
+//        for (int i = 0; i < nrow; i++) {
+//            answer->at(i) = this->at(i)->timesTransposeFullMatrixForFMFMDsptr(fullMat);
+//        }
+//        return answer;
+//    }
     std::shared_ptr<FullMatrixDouble> FullMatrixDouble::times(double a)
     {
         int m = this->nrow();
diff --git a/OndselSolver/FullMatrix.h b/OndselSolver/FullMatrix.h
index d0d53d1..2648732 100644
--- a/OndselSolver/FullMatrix.h
+++ b/OndselSolver/FullMatrix.h
@@ -119,7 +119,7 @@ namespace MbD {
         double maxMagnitude() override;
         void zeroSelf() override;
         std::shared_ptr<FullMatrixFullMatrixDouble> times(double a);
-        std::shared_ptr<FullMatrixFullMatrixDouble> timesTransposeFullMatrix(std::shared_ptr<FullMatrixFullMatrixDouble> fullMat);
+        // std::shared_ptr<FullMatrixFullMatrixDouble> timesTransposeFullMatrix(std::shared_ptr<FullMatrixFullMatrixDouble> fullMat);
         double sumOfSquares() override;
         void identity();
         static std::shared_ptr<MbD::FullMatrixFullMatrixDouble> identitysptr(int n);
diff --git a/OndselSolver/FullRow.cpp b/OndselSolver/FullRow.cpp
index 2c3d2a9..2716599 100644
--- a/OndselSolver/FullRow.cpp
+++ b/OndselSolver/FullRow.cpp
@@ -7,6 +7,7 @@
  ***************************************************************************/
  
 #include "FullRow.h"
+#include "FullMatrix.h"
 
 namespace MbD {
     template<typename T>
@@ -22,45 +23,37 @@ namespace MbD {
         return answer;
     }
 
-    template<typename T>
-    FRowsptr<T> FullRow<T>::timesTransposeFullMatrix(std::shared_ptr<FullMatrixDouble> fullMat)
+    template<>
+    FRowsptr<double> FullRow<double>::timesTransposeFullMatrix(std::shared_ptr<FullMatrixDouble> fullMat)
     {
         //"a*bT = a(1,j)b(k,j)"
         int ncol = fullMat->nrow();
-        auto answer = std::make_shared<FullRow<T>>(ncol);
+        auto answer = std::make_shared<FullRow<double>>(ncol);
         for (int k = 0; k < ncol; k++) {
             answer->at(k) = this->dot(fullMat->at(k));
         }
         return answer;
     }
 
-    template<typename T>
-    FRowsptr<T> FullRow<T>::timesTransposeFullMatrixForFMFMDsptr(std::shared_ptr<FullMatrixFullMatrixDouble> fullMat)
-    {
-        //"a*bT = a(1,j)b(k,j)"
-        int ncol = fullMat->nrow();
-        auto answer = std::make_shared<FullRow<T>>(ncol);
-        for (int k = 0; k < ncol; k++) {
-            answer->at(k) = this->dot(fullMat->at(k));
-        }
-        return answer;
-    }
+// TODO: can't get the following to work, but CLion says the routine that calls it in FullMatrixFullMatrixDouble is also
+//       never called.
+//    template<>
+//    FRowsptr<std::shared_ptr<FullMatrixDouble>> FullRow<std::shared_ptr<FullMatrixDouble>>::timesTransposeFullMatrixForFMFMDsptr(
+//            std::shared_ptr<FullMatrixFullMatrixDouble> fullMat)
+//    {
+//        //"a*bT = a(1,j)b(k,j)"
+//        int ncol = fullMat->nrow();
+//        auto answer = std::make_shared<FullRow<std::shared_ptr<FullMatrixDouble>>>(ncol);
+//        for (int k = 0; k < ncol; k++) {
+//            answer->at(k) = this->dot(fullMat->at(k));
+//        }
+//        return answer;
+//    }
 
-    template<typename T>
-    FRowsptr<T> FullRow<T>::timesFullMatrix(std::shared_ptr<FullMatrixDouble> fullMat)
+    template<>
+    FRowsptr<double> FullRow<double>::timesFullMatrix(std::shared_ptr<FullMatrixDouble> fullMat)
     {
-        FRowsptr<T> answer = fullMat->at(0)->times(this->at(0));
-        for (int j = 1; j < (int) this->size(); j++)
-        {
-            answer->equalSelfPlusFullRowTimes(fullMat->at(j), this->at(j));
-        }
-        return answer;
-    }
-
-    template<typename T>
-    FRowsptr<T> FullRow<T>::timesFullMatrixForFMFMDsptr(std::shared_ptr<FullMatrixFullMatrixDouble> fullMat)
-    {
-        FRowsptr<T> answer = fullMat->at(0)->times(this->at(0));
+        FRowsptr<double> answer = fullMat->at(0)->times(this->at(0));
         for (int j = 1; j < (int) this->size(); j++)
         {
             answer->equalSelfPlusFullRowTimes(fullMat->at(j), this->at(j));
diff --git a/OndselSolver/FullRow.h b/OndselSolver/FullRow.h
index 0ca02f4..8d9ecca 100644
--- a/OndselSolver/FullRow.h
+++ b/OndselSolver/FullRow.h
@@ -37,13 +37,10 @@ namespace MbD {
 		void atiplusFullRow(int j, FRowsptr<T> fullRow);
 		std::ostream& printOn(std::ostream& s) const override;
 
-        FRowsptr<T> timesFullMatrix(std::shared_ptr<FullMatrixDouble> fullMat);
-        FRowsptr<T> timesTransposeFullMatrix(std::shared_ptr<FullMatrixDouble> fullMat);
         std::shared_ptr<FullMatrixDouble> transposeTimesFullRow(FRowsptr<T> fullRow);
-
-        FRowsptr<T> timesFullMatrixForFMFMDsptr(std::shared_ptr<FullMatrixFullMatrixDouble> fullMat);
-        FRowsptr<T> timesTransposeFullMatrixForFMFMDsptr(std::shared_ptr<FullMatrixFullMatrixDouble> fullMat);
-        // std::shared_ptr<FullMatrixFullMatrixDouble> transposeTimesFullRow(FRowsptr<T> fullRow);
+        FRowsptr<double> timesTransposeFullMatrix(std::shared_ptr<FullMatrixDouble> fullMat);
+        // FRowsptr<std::shared_ptr<FullMatrixDouble>> timesTransposeFullMatrixForFMFMDsptr(std::shared_ptr<FullMatrixFullMatrixDouble> fullMat);
+        FRowsptr<double> timesFullMatrix(std::shared_ptr<FullMatrixDouble> fullMat);
     };
 
 	template<>
@@ -148,6 +145,6 @@ namespace MbD {
 		s << "}";
 		return s;
 	};
-    // TODO: template class FullRow<double>;
+
 }
 
diff --git a/OndselSolver/SparseMatrix.h b/OndselSolver/SparseMatrix.h
index eea8d69..f07232f 100644
--- a/OndselSolver/SparseMatrix.h
+++ b/OndselSolver/SparseMatrix.h
@@ -13,6 +13,7 @@
 #include "RowTypeMatrix.h"
 #include "SparseRow.h"
 #include "DiagonalMatrix.h"
+#include "FullMatrix.h"
 
 namespace MbD {
 	template<typename T>