From ca0393dabc83020a617fe6b24305d7d7364b9402 Mon Sep 17 00:00:00 2001
From: Jesper Andersson <jesper@fmrib.ox.ac.uk>
Date: Mon, 7 Apr 2008 17:30:21 +0000
Subject: [PATCH] BFSparseMatrix now templated to allow double or float
representation
---
bfmatrix.cpp | 382 ++++++++++++++-------------------------------------
bfmatrix.h | 250 ++++++++++++++++++++++++++++++---
2 files changed, 334 insertions(+), 298 deletions(-)
diff --git a/bfmatrix.cpp b/bfmatrix.cpp
index 9f59d9e..9305b73 100644
--- a/bfmatrix.cpp
+++ b/bfmatrix.cpp
@@ -17,216 +17,6 @@
namespace MISCMATHS {
-//
-// Member functions for SparseBFMatrix
-//
-
-//
-// Transpose
-//
-
-boost::shared_ptr<BFMatrix> SparseBFMatrix::Transpose()
-const
-{
- boost::shared_ptr<SparseBFMatrix> tm(new SparseBFMatrix(mp->t()));
- return(tm);
-}
-
-//
-// Concatenation of two matrices returning a third
-//
-void SparseBFMatrix::HorConcat(const BFMatrix& B, BFMatrix& AB) const
-{
- if (B.Nrows() && Nrows() != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
-
- try {
- SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB); // Means that output is sparse
- lAB = *this;
- lAB.HorConcat2MyRight(B);
- }
- catch (std::bad_cast) {
- try {
- FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB); // Means that output is full
- lAB = FullBFMatrix(this->AsMatrix());
- lAB.HorConcat2MyRight(B);
- }
- catch (std::bad_cast) {
- throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error");
- }
- }
-}
-
-void SparseBFMatrix::HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
-{
- if (B.Nrows() && int(Nrows()) != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
-
- try {
- SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB); // Means that output is sparse
- lAB = *this;
- lAB.HorConcat2MyRight(B);
- }
- catch (std::bad_cast) {
- try {
- FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB); // Means that output is full
- lAB = FullBFMatrix(this->AsMatrix());
- lAB.HorConcat2MyRight(B);
- }
- catch (std::bad_cast) {
- throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error");
- }
- }
-}
-
-void SparseBFMatrix::VertConcat(const BFMatrix& B, BFMatrix& AB) const
-{
- if (B.Ncols() && Ncols() != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
-
- try {
- SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB); // Means that output is sparse
- lAB = *this;
- lAB.VertConcatBelowMe(B);
- }
- catch (std::bad_cast) {
- try {
- FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB); // Means that output is full
- lAB = FullBFMatrix(this->AsMatrix());
- lAB.VertConcatBelowMe(B);
- }
- catch (std::bad_cast) {
- throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error");
- }
- }
-}
-
-void SparseBFMatrix::VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
-{
- if (B.Ncols() && int(Ncols()) != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
-
- try {
- SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB); // Means that output is sparse
- lAB = *this;
- lAB.VertConcatBelowMe(B);
- }
- catch (std::bad_cast) {
- try {
- FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB); // Means that output is full
- lAB = FullBFMatrix(this->AsMatrix());
- lAB.VertConcatBelowMe(B);
- }
- catch (std::bad_cast) {
- throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error");
- }
- }
-}
-
-//
-// Concatenate another matrix to *this
-//
-void SparseBFMatrix::HorConcat2MyRight(const BFMatrix& B)
-{
- if (!B.Nrows()) return;
-
- if (Nrows() != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
-
- try {
- const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B); // Means that we want to concatenate a sparse matrix
- *mp |= *(lB.mp);
- }
- catch (std::bad_cast) {
- try {
- const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B); // Means that we want to concatenate a full
- this->HorConcat2MyRight(lB.AsMatrix());
- }
- catch (std::bad_cast) {
- throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: dynamic cast error");
- }
- }
-}
-
-void SparseBFMatrix::HorConcat2MyRight(const NEWMAT::Matrix& B)
-{
- if (!B.Nrows()) return;
-
- if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
- *mp |= B;
-}
-
-void SparseBFMatrix::VertConcatBelowMe(const BFMatrix& B)
-{
- if (!B.Ncols()) return;
-
- if (Ncols() != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
-
- try {
- const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B); // Means that we want to concatenate a sparse matrix
- *mp &= *(lB.mp);
- }
- catch (std::bad_cast) {
- try {
- const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B); // Means that we want to concatenate a full
- this->VertConcatBelowMe(lB.AsMatrix());
- }
- catch (std::bad_cast) {
- throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: dynamic cast error");
- }
- }
-}
-
-void SparseBFMatrix::VertConcatBelowMe(const NEWMAT::Matrix& B)
-{
- if (!B.Ncols()) return;
-
- if (int(Ncols()) != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
- *mp &= B;
-}
-
-// Multiply by vector
-NEWMAT::ReturnMatrix SparseBFMatrix::MulByVec(const NEWMAT::ColumnVector& invec) const
-{
- if (invec.Nrows() != int(Ncols())) {throw BFMatrixException("Matrix-vector size mismatch");}
- NEWMAT::ColumnVector outvec = *mp * invec;
- outvec.Release();
- return(outvec);
-}
-
-// Add another matrix to this one
-void SparseBFMatrix::AddToMe(const BFMatrix& M, double s)
-{
- if (Ncols() != M.Ncols() || Nrows() != M.Nrows()) {
- throw BFMatrixException("SparseBFMatrix::AddToMe: Matrix size mismatch");
- }
-
- try {
- const SparseBFMatrix& lM = dynamic_cast<const SparseBFMatrix&>(M); // Add sparse matrix to this sparse matrix
- if (s == 1.0) *mp += *(lM.mp);
- else *mp += s * *(lM.mp);
- }
- catch (std::bad_cast) {
- try {
- const FullBFMatrix& lM = dynamic_cast<const FullBFMatrix&>(M); // Add full matrix to this sparse matrix
- if (s == 1.0) *mp += SpMat<double>(lM.ReadAsMatrix());
- else *mp += s * SpMat<double>(lM.ReadAsMatrix());
- }
- catch (std::bad_cast) {
- throw BFMatrixException("SparseBFMatrix::AddToMe: dynamic cast error");
- }
- }
-}
-
-// Given A*x=b, solve for x
-NEWMAT::ReturnMatrix SparseBFMatrix::SolveForx(const NEWMAT::ColumnVector& b,
- MISCMATHS::MatrixType type,
- double tol,
- int miter) const
-{
- if (b.Nrows() != int(Nrows())) {
- throw BFMatrixException("SparseBFMatrix::SolveForx: Matrix-vector size mismatch");
- }
- NEWMAT::ColumnVector x = mp->SolveForx(b,type,tol,miter);
- x.Release();
- return(x);
-}
-
//
// Member functions for FullBFMatrix
//
@@ -252,19 +42,24 @@ void FullBFMatrix::HorConcat(const BFMatrix& B, BFMatrix& AB) const
{
if (B.Nrows() && Nrows() != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat: Matrices must have same # of rows");}
- try {
- FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB); // This means output is a full matrix
- lAB = *this;
- lAB.HorConcat2MyRight(B);
+ FullBFMatrix *pAB = dynamic_cast<FullBFMatrix *>(&AB);
+ if (pAB) { // This means output is a full matrix
+ *pAB = *this;
+ pAB->HorConcat2MyRight(B);
}
- catch (std::bad_cast) {
- try {
- SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
- lAB = SparseBFMatrix(this->AsMatrix());
- lAB.HorConcat2MyRight(B);
+ else {
+ SparseBFMatrix<double> *psdAB = dynamic_cast<SparseBFMatrix<double> *>(&AB);
+ if (psdAB) {
+ *psdAB = SparseBFMatrix<double>(this->AsMatrix());
+ psdAB->HorConcat2MyRight(B);
}
- catch (std::bad_cast) {
- throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error");
+ else {
+ SparseBFMatrix<float> *psfAB = dynamic_cast<SparseBFMatrix<float> *>(&AB);
+ if (psfAB) {
+ *psfAB = SparseBFMatrix<float>(this->AsMatrix());
+ psfAB->HorConcat2MyRight(B);
+ }
+ else throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error");
}
}
}
@@ -273,19 +68,24 @@ void FullBFMatrix::HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
{
if (B.Nrows() && int(Nrows()) != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat: Matrices must have same # of rows");}
- try {
- FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB); // This means output is a full matrix
- lAB = *this;
- lAB.HorConcat2MyRight(B);
+ FullBFMatrix *pAB = dynamic_cast<FullBFMatrix *>(&AB);
+ if (pAB) { // This means output is a full matrix
+ *pAB = *this;
+ pAB->HorConcat2MyRight(B);
}
- catch (std::bad_cast) {
- try {
- SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
- lAB = SparseBFMatrix(this->AsMatrix());
- lAB.HorConcat2MyRight(B);
+ else {
+ SparseBFMatrix<double> *psdAB = dynamic_cast<SparseBFMatrix<double> *>(&AB);
+ if (psdAB) {
+ *psdAB = SparseBFMatrix<double>(this->AsMatrix());
+ psdAB->HorConcat2MyRight(B);
}
- catch (std::bad_cast) {
- throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error");
+ else {
+ SparseBFMatrix<float> *psfAB = dynamic_cast<SparseBFMatrix<float> *>(&AB);
+ if (psfAB) {
+ *psfAB = SparseBFMatrix<float>(this->AsMatrix());
+ psfAB->HorConcat2MyRight(B);
+ }
+ else throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error");
}
}
}
@@ -294,19 +94,24 @@ void FullBFMatrix::VertConcat(const BFMatrix& B, BFMatrix& AB) const
{
if (B.Ncols() && Ncols() != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcat: Matrices must have same # of columns");}
- try {
- FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB); // This means output is a full matrix
- lAB = *this;
- lAB.VertConcatBelowMe(B);
+ FullBFMatrix *pAB = dynamic_cast<FullBFMatrix *>(&AB);
+ if (pAB) { // This means output is a full matrix
+ *pAB = *this;
+ pAB->VertConcatBelowMe(B);
}
- catch (std::bad_cast) {
- try {
- SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
- lAB = SparseBFMatrix(this->AsMatrix());
- lAB.VertConcatBelowMe(B);
+ else {
+ SparseBFMatrix<double> *psdAB = dynamic_cast<SparseBFMatrix<double> *>(&AB);
+ if (psdAB) {
+ *psdAB = SparseBFMatrix<double>(this->AsMatrix());
+ psdAB->VertConcatBelowMe(B);
}
- catch (std::bad_cast) {
- throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error");
+ else {
+ SparseBFMatrix<float> *psfAB = dynamic_cast<SparseBFMatrix<float> *>(&AB);
+ if (psfAB) {
+ *psfAB = SparseBFMatrix<float>(this->AsMatrix());
+ psfAB->VertConcatBelowMe(B);
+ }
+ else throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error");
}
}
}
@@ -315,19 +120,24 @@ void FullBFMatrix::VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
{
if (B.Ncols() && int(Ncols()) != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcat: Matrices must have same # of columns");}
- try {
- FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB); // This means output is a full matrix
- lAB = *this;
- lAB.VertConcatBelowMe(B);
+ FullBFMatrix *pAB = dynamic_cast<FullBFMatrix *>(&AB);
+ if (pAB) { // This means output is a full matrix
+ *pAB = *this;
+ pAB->VertConcatBelowMe(B);
}
- catch (std::bad_cast) {
- try {
- SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
- lAB = SparseBFMatrix(this->AsMatrix());
- lAB.VertConcatBelowMe(B);
+ else {
+ SparseBFMatrix<double> *psdAB = dynamic_cast<SparseBFMatrix<double> *>(&AB);
+ if (psdAB) {
+ *psdAB = SparseBFMatrix<double>(this->AsMatrix());
+ psdAB->VertConcatBelowMe(B);
}
- catch (std::bad_cast) {
- throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error");
+ else {
+ SparseBFMatrix<float> *psfAB = dynamic_cast<SparseBFMatrix<float> *>(&AB);
+ if (psfAB) {
+ *psfAB = SparseBFMatrix<float>(this->AsMatrix());
+ psfAB->VertConcatBelowMe(B);
+ }
+ else throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error");
}
}
}
@@ -341,17 +151,21 @@ void FullBFMatrix::HorConcat2MyRight(const BFMatrix& B)
if (Nrows() != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
- try {
- const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
- *mp |= *(lB.mp);
+ const FullBFMatrix *pB = dynamic_cast<const FullBFMatrix *>(&B);
+ if (pB) { // If B was full
+ *mp |= *(pB->mp);
}
- catch (std::bad_cast) {
- try {
- const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
- this->HorConcat2MyRight(lB.AsMatrix());
+ else {
+ const SparseBFMatrix<double> *psdB = dynamic_cast<const SparseBFMatrix<double> *>(&B);
+ if (psdB) {
+ this->HorConcat2MyRight(psdB->AsMatrix());
}
- catch (std::bad_cast) {
- throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: dynamic cast error");
+ else {
+ const SparseBFMatrix<float> *psfB = dynamic_cast<const SparseBFMatrix<float> *>(&B);
+ if (psfB) {
+ this->HorConcat2MyRight(psfB->AsMatrix());
+ }
+ else throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: dynamic cast error");
}
}
}
@@ -370,17 +184,21 @@ void FullBFMatrix::VertConcatBelowMe(const BFMatrix& B)
if (Ncols() != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
- try {
- const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
- *mp &= *(lB.mp);
+ const FullBFMatrix *pB = dynamic_cast<const FullBFMatrix *>(&B);
+ if (pB) { // Means B is full
+ *mp &= *(pB->mp);
}
- catch (std::bad_cast) {
- try {
- const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
- this->VertConcatBelowMe(lB.AsMatrix());
+ else {
+ const SparseBFMatrix<double> *psdB = dynamic_cast<const SparseBFMatrix<double> *>(&B);
+ if (psdB) {
+ this->VertConcatBelowMe(psdB->AsMatrix());
}
- catch (std::bad_cast) {
- throw BFMatrixException("FullBFMatrix::HorConcatBelowMe: dynamic cast error");
+ else {
+ const SparseBFMatrix<float> *psfB = dynamic_cast<const SparseBFMatrix<float> *>(&B);
+ if (psfB) {
+ this->VertConcatBelowMe(psfB->AsMatrix());
+ }
+ else throw BFMatrixException("FullBFMatrix::HorConcatBelowMe: dynamic cast error");
}
}
}
@@ -417,17 +235,17 @@ void FullBFMatrix::AddToMe(const BFMatrix& m, double s)
throw BFMatrixException("FullBFMatrix::AddToMe: Matrix size mismatch");
}
- try {
- const FullBFMatrix& lm = dynamic_cast<const FullBFMatrix&>(m);
- *mp += s*(*lm.mp);
+ const FullBFMatrix *pm = dynamic_cast<const FullBFMatrix *>(&m);
+ if (pm) { // If m is full matrix
+ *mp += s*(*(pm->mp));
}
- catch (std::bad_cast) {
- try {
- const SparseBFMatrix& lm = dynamic_cast<const SparseBFMatrix&>(m);
- *mp += s*lm.AsMatrix();
- }
- catch (std::bad_cast) {
- throw BFMatrixException("FullBFMatrix::AddToMe: dynamic cast error");
+ else {
+ const SparseBFMatrix<double> *psdm = dynamic_cast<const SparseBFMatrix<double> *>(&m);
+ if (psdm) *mp += s*psdm->AsMatrix();
+ else {
+ const SparseBFMatrix<float> *psfm = dynamic_cast<const SparseBFMatrix<float> *>(&m);
+ if (psfm) *mp += s*psfm->AsMatrix();
+ else throw BFMatrixException("FullBFMatrix::AddToMe: dynamic cast error");
}
}
}
diff --git a/bfmatrix.h b/bfmatrix.h
index a6032a7..604f92c 100644
--- a/bfmatrix.h
+++ b/bfmatrix.h
@@ -49,6 +49,8 @@ public:
~BFMatrixException() throw() {}
};
+enum BFMatrixPrecisionType {BFMatrixDoublePrecision, BFMatrixFloatPrecision};
+
class BFMatrix
{
protected:
@@ -70,7 +72,8 @@ public:
virtual void Print(const std::string fname=std::string("")) const = 0;
// Setting, deleting or resizing the whole sparse matrix.
- virtual void SetMatrix(const MISCMATHS::SpMat<double>& M) = 0;
+ // virtual void SetMatrix(const MISCMATHS::SpMat<double>& M) = 0;
+ // virtual void SetMatrix(const MISCMATHS::SpMat<float>& M) = 0;
virtual void SetMatrix(const NEWMAT::Matrix& M) = 0;
virtual void Clear() = 0;
virtual void Resize(unsigned int m, unsigned int n) = 0;
@@ -114,26 +117,27 @@ public:
int miter=200) const = 0;
};
+template<class T>
class SparseBFMatrix : public BFMatrix
{
private:
- boost::shared_ptr<MISCMATHS::SpMat<double> > mp;
-
+ boost::shared_ptr<MISCMATHS::SpMat<T> > mp;
+
public:
// Constructors, destructor and assignment
SparseBFMatrix()
- : mp(boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>())) {}
+ : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>())) {}
SparseBFMatrix(unsigned int m, unsigned int n)
- : mp(boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(m,n))) {}
+ : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(m,n))) {}
SparseBFMatrix(unsigned int m, unsigned int n, const unsigned int *irp, const unsigned int *jcp, const double *sp)
- : mp(boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(m,n,irp,jcp,sp))) {}
- SparseBFMatrix(const MISCMATHS::SpMat<double>& M)
- : mp(boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(M))) {}
+ : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(m,n,irp,jcp,sp))) {}
+ SparseBFMatrix(const MISCMATHS::SpMat<T>& M)
+ : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M))) {}
SparseBFMatrix(const NEWMAT::Matrix& M)
- : mp(boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(M))) {}
+ : mp(boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M))) {}
virtual ~SparseBFMatrix() {}
- virtual const SparseBFMatrix& operator=(const SparseBFMatrix& M) {
- mp = boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(*(M.mp))); return(*this);
+ virtual const SparseBFMatrix& operator=(const SparseBFMatrix<T>& M) {
+ mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(*(M.mp))); return(*this);
}
// Access as NEWMAT::Matrix
@@ -147,11 +151,12 @@ public:
virtual void Print(const std::string fname=std::string("")) const {mp->Print(fname);}
// Setting, deleting or resizing the whole sparse matrix.
- virtual void SetMatrix(const MISCMATHS::SpMat<double>& M) {mp = boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(M));}
- virtual void SetMatrix(const NEWMAT::Matrix& M) {mp = boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(M));}
- virtual void SetMatrixPtr(boost::shared_ptr<MISCMATHS::SpMat<double> >& mptr) {mp = mptr;}
- virtual void Clear() {mp = boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>());}
- virtual void Resize(unsigned int m, unsigned int n) {mp = boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(m,n));}
+ virtual void SetMatrix(const MISCMATHS::SpMat<T>& M) {mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M));}
+ // virtual void SetMatrix(const MISCMATHS::SpMat<float>& M) {mp = boost::shared_ptr<MISCMATHS::SpMat<float> >(new MISCMATHS::SpMat<float>(M));}
+ virtual void SetMatrix(const NEWMAT::Matrix& M) {mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(M));}
+ virtual void SetMatrixPtr(boost::shared_ptr<MISCMATHS::SpMat<T> >& mptr) {mp = mptr;}
+ virtual void Clear() {mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>());}
+ virtual void Resize(unsigned int m, unsigned int n) {mp = boost::shared_ptr<MISCMATHS::SpMat<T> >(new MISCMATHS::SpMat<T>(m,n));}
// Accessing values
virtual double Peek(unsigned int r, unsigned int c) const {return(mp->Peek(r,c));}
@@ -221,6 +226,7 @@ public:
// Setting, deleting or resizing the whole matrix.
virtual void SetMatrix(const MISCMATHS::SpMat<double>& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M.AsNEWMAT()));}
+ virtual void SetMatrix(const MISCMATHS::SpMat<float>& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M.AsNEWMAT()));}
virtual void SetMatrix(const NEWMAT::Matrix& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M));}
virtual void SetMatrixPtr(boost::shared_ptr<NEWMAT::Matrix>& mptr) {mp = mptr;}
virtual void Clear() {mp->ReSize(0,0);}
@@ -266,6 +272,218 @@ public:
};
+//
+// Here comes member functions for SparseBFMatrix. Since it is templated
+// these need to go here rather than in bfmatrix.cpp.
+//
+
+//
+// Member functions for SparseBFMatrix
+//
+
+//
+// Transpose
+//
+
+template<class T>
+boost::shared_ptr<BFMatrix> SparseBFMatrix<T>::Transpose() const
+{
+ boost::shared_ptr<SparseBFMatrix<T> > tm(new SparseBFMatrix<T>(mp->t()));
+ return(tm);
+}
+
+//
+// Concatenation of two matrices returning a third
+//
+template<class T>
+void SparseBFMatrix<T>::HorConcat(const BFMatrix& B, BFMatrix& AB) const
+{
+ if (B.Nrows() && Nrows() != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
+
+ SparseBFMatrix<T> *pAB = dynamic_cast<SparseBFMatrix<T> *>(&AB);
+ if (pAB) { // Means that output is sparse of type T
+ *pAB = *this;
+ pAB->HorConcat2MyRight(B);
+ }
+ else {
+ FullBFMatrix *fpAB = dynamic_cast<FullBFMatrix *>(&AB);
+ if (fpAB) { // Means that output is full
+ *fpAB = FullBFMatrix(this->AsMatrix());
+ fpAB->HorConcat2MyRight(B);
+ }
+ else throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error");
+ }
+}
+
+template<class T>
+void SparseBFMatrix<T>::HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
+{
+ if (B.Nrows() && int(Nrows()) != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
+
+ SparseBFMatrix<T> *pAB = dynamic_cast<SparseBFMatrix<T> *>(&AB);
+ if (pAB) { // Means that output is sparse
+ *pAB = *this;
+ pAB->HorConcat2MyRight(B);
+ }
+ else {
+ FullBFMatrix *fpAB = dynamic_cast<FullBFMatrix *>(&AB);
+ if (fpAB) {// Means that output is full
+ *fpAB = FullBFMatrix(this->AsMatrix());
+ fpAB->HorConcat2MyRight(B);
+ }
+ else throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error");
+ }
+}
+
+template<class T>
+void SparseBFMatrix<T>::VertConcat(const BFMatrix& B, BFMatrix& AB) const
+{
+ if (B.Ncols() && Ncols() != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
+
+ SparseBFMatrix<T> *pAB = dynamic_cast<SparseBFMatrix<T> *>(&AB);
+ if (pAB) { // Means that output is sparse
+ *pAB = *this;
+ pAB->VertConcatBelowMe(B);
+ }
+ else {
+ FullBFMatrix *fpAB = dynamic_cast<FullBFMatrix *>(&AB);
+ if (fpAB) { // Means that output is full
+ *fpAB = FullBFMatrix(this->AsMatrix());
+ fpAB->VertConcatBelowMe(B);
+ }
+ else throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error");
+ }
+}
+
+template<class T>
+void SparseBFMatrix<T>::VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
+{
+ if (B.Ncols() && int(Ncols()) != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
+
+ SparseBFMatrix<T> *pAB = dynamic_cast<SparseBFMatrix<T> *>(&AB);
+ if (pAB) { // Means that output is sparse
+ *pAB = *this;
+ pAB->VertConcatBelowMe(B);
+ }
+ else {
+ FullBFMatrix *fpAB = dynamic_cast<FullBFMatrix *>(&AB);
+ if (fpAB) { // Means that output is full
+ *fpAB = FullBFMatrix(this->AsMatrix());
+ fpAB->VertConcatBelowMe(B);
+ }
+ else throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error");
+ }
+}
+
+//
+// Concatenate another matrix to *this
+//
+template<class T>
+void SparseBFMatrix<T>::HorConcat2MyRight(const BFMatrix& B)
+{
+ if (!B.Nrows()) return;
+
+ if (Nrows() != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
+
+ const SparseBFMatrix<T> *pB = dynamic_cast<const SparseBFMatrix<T> *>(&B);
+ if (pB) { // Means that we want to concatenate a sparse matrix
+ *mp |= *(pB->mp);
+ }
+ else {
+ const FullBFMatrix *fpB = dynamic_cast<const FullBFMatrix *>(&B);
+ if (fpB) { // Means that we want to concatenate a full
+ this->HorConcat2MyRight(fpB->AsMatrix());
+ }
+ else throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: dynamic cast error");
+ }
+}
+
+template<class T>
+void SparseBFMatrix<T>::HorConcat2MyRight(const NEWMAT::Matrix& B)
+{
+ if (!B.Nrows()) return;
+
+ if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
+ *mp |= B;
+}
+
+template<class T>
+void SparseBFMatrix<T>::VertConcatBelowMe(const BFMatrix& B)
+{
+ if (!B.Ncols()) return;
+
+ if (Ncols() != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
+
+ const SparseBFMatrix<T> *pB = dynamic_cast<const SparseBFMatrix<T> *>(&B);
+ if (pB) { // Means that we want to concatenate a sparse matrix
+ *mp &= *(pB->mp);
+ }
+ else {
+ const FullBFMatrix *fpB = dynamic_cast<const FullBFMatrix *>(&B);
+ if (fpB) { // Means that we want to concatenate a full
+ this->VertConcatBelowMe(fpB->AsMatrix());
+ }
+ else throw BFMatrixException("SparseBFMatrix::VertConcatBelowMe: dynamic cast error");
+ }
+}
+
+template<class T>
+void SparseBFMatrix<T>::VertConcatBelowMe(const NEWMAT::Matrix& B)
+{
+ if (!B.Ncols()) return;
+
+ if (int(Ncols()) != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
+ *mp &= B;
+}
+
+// Multiply by vector
+template<class T>
+NEWMAT::ReturnMatrix SparseBFMatrix<T>::MulByVec(const NEWMAT::ColumnVector& invec) const
+{
+ if (invec.Nrows() != int(Ncols())) {throw BFMatrixException("Matrix-vector size mismatch");}
+ NEWMAT::ColumnVector outvec = *mp * invec;
+ outvec.Release();
+ return(outvec);
+}
+
+// Add another matrix to this one
+template<class T>
+void SparseBFMatrix<T>::AddToMe(const BFMatrix& M, double s)
+{
+ if (Ncols() != M.Ncols() || Nrows() != M.Nrows()) {
+ throw BFMatrixException("SparseBFMatrix::AddToMe: Matrix size mismatch");
+ }
+
+ const SparseBFMatrix<T> *pM = dynamic_cast<const SparseBFMatrix<T> *>(&M);
+ if (pM) { // Add sparse matrix to this sparse matrix
+ if (s == 1.0) *mp += *(pM->mp);
+ else *mp += s * *(pM->mp);
+ }
+ else {
+ const FullBFMatrix *fpM = dynamic_cast<const FullBFMatrix *>(&M);
+ if (fpM) { // Add full matrix to this sparse matrix
+ if (s == 1.0) *mp += SpMat<T>(fpM->ReadAsMatrix());
+ else *mp += s * SpMat<T>(fpM->ReadAsMatrix());
+ }
+ else throw BFMatrixException("SparseBFMatrix::AddToMe: dynamic cast error");
+ }
+}
+
+// Given A*x=b, solve for x
+template<class T>
+NEWMAT::ReturnMatrix SparseBFMatrix<T>::SolveForx(const NEWMAT::ColumnVector& b,
+ MISCMATHS::MatrixType type,
+ double tol,
+ int miter) const
+{
+ if (b.Nrows() != int(Nrows())) {
+ throw BFMatrixException("SparseBFMatrix::SolveForx: Matrix-vector size mismatch");
+ }
+ NEWMAT::ColumnVector x = mp->SolveForx(b,type,tol,miter);
+ x.Release();
+ return(x);
+}
+
} // End namespace MISCMATHS
#endif // End #ifndef BFMatrix_h
--
GitLab