From 5b02f4fbe8586c82b5f12b5a236a5c19e3b9c50a Mon Sep 17 00:00:00 2001
From: Jesper Andersson <jesper@fmrib.ox.ac.uk>
Date: Mon, 4 Jun 2007 14:03:00 +0000
Subject: [PATCH] Implements generic full/sparse matrix for use in nonlin
---
bfmatrix.cpp | 321 +++++++++++++++++++++++++++++++++++++++++++++++++++
bfmatrix.h | 270 +++++++++++++++++++++++++++++++++++++++++++
2 files changed, 591 insertions(+)
create mode 100644 bfmatrix.cpp
create mode 100644 bfmatrix.h
diff --git a/bfmatrix.cpp b/bfmatrix.cpp
new file mode 100644
index 0000000..480d351
--- /dev/null
+++ b/bfmatrix.cpp
@@ -0,0 +1,321 @@
+//
+// Definitions for class BFMatrix
+//
+
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <boost/shared_ptr.hpp>
+#include "newmat.h"
+#include "newmatio.h"
+#include "bfmatrix.h"
+
+namespace MISCMATHS {
+
+//
+// Member functions for BFMatrix
+//
+
+void BFMatrix::print(const NEWMAT::Matrix& m,
+ const std::string& fname) const
+{
+ if (!fname.length()) {
+ cout << endl << m << endl;
+ }
+ else {
+ try {
+ ofstream fout(fname.c_str());
+ fout << setprecision(10) << m;
+ }
+ catch(...) {
+ std::string errmsg("BFMatrix::print: Failed to write to file " + fname);
+ throw BFMatrixException(errmsg);
+ }
+ }
+}
+
+//
+// Member functions for SparseBFMatrix
+//
+
+//
+// Concatenation of two matrices returning a third
+//
+void SparseBFMatrix::HorConcat(const BFMatrix& B, BFMatrix& AB) const
+{
+ try {
+ const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
+ SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
+ if (Nrows() != lB.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
+ *(lAB.mp) = *mp | *(lB.mp);
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error");
+ }
+}
+
+void SparseBFMatrix::HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
+{
+ try {
+ SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
+ if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat: Matrices must have same # of rows");}
+ *(lAB.mp) = *mp | B;
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("SparseBFMatrix::HorConcat: dynamic cast error");
+ }
+}
+
+void SparseBFMatrix::VertConcat(const BFMatrix& B, BFMatrix& AB) const
+{
+ try {
+ const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
+ SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
+ if (Ncols() != lB.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
+ *(lAB.mp) = *mp & *(lB.mp);
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error");
+ }
+}
+
+void SparseBFMatrix::VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
+{
+ try {
+ SparseBFMatrix& lAB = dynamic_cast<SparseBFMatrix&>(AB);
+ if (int(Ncols()) != B.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcat: Matrices must have same # of columns");}
+ *(lAB.mp) = *mp & B;
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error");
+ }
+}
+
+//
+// Concatenate another matrix to *this
+//
+void SparseBFMatrix::HorConcat2MyRight(const BFMatrix& B)
+{
+ try {
+ const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
+ if (Nrows() != lB.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
+ *mp |= *(lB.mp);
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error");
+ }
+}
+
+void SparseBFMatrix::HorConcat2MyRight(const NEWMAT::Matrix& B)
+{
+ if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("SparseBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
+ *mp |= B;
+}
+
+void SparseBFMatrix::VertConcatBelowMe(const BFMatrix& B)
+{
+ try {
+ const SparseBFMatrix& lB = dynamic_cast<const SparseBFMatrix&>(B);
+ if (Ncols() != lB.Ncols()) {throw BFMatrixException("SparseBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
+ *mp &= *(lB.mp);
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("SparseBFMatrix::VertConcat: dynamic cast error");
+ }
+}
+
+void SparseBFMatrix::VertConcatBelowMe(const NEWMAT::Matrix& B)
+{
+ 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)
+{
+ try {
+ const SparseBFMatrix& lM = dynamic_cast<const SparseBFMatrix&>(M);
+ if (Ncols() != lM.Ncols() || Nrows() != lM.Nrows()) {
+ throw BFMatrixException("SparseBFMatrix::AddToMe: Matrix size mismatch");
+ }
+ if (s == 1.0) *mp += *(lM.mp);
+ else *mp += s * *(lM.mp);
+ }
+ 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
+//
+
+boost::shared_ptr<BFMatrix> FullBFMatrix::Transpose(boost::shared_ptr<BFMatrix>& pA) const
+{
+ boost::shared_ptr<FullBFMatrix> tmp(new FullBFMatrix);
+ *(tmp->mp) = mp->t();
+ return(tmp);
+}
+
+//
+// Concatenate two matrices yielding a third
+//
+
+void FullBFMatrix::HorConcat(const BFMatrix& B, BFMatrix& AB) const
+{
+ try {
+ const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
+ FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);
+ if (Nrows() != lB.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat: Matrices must have same # of rows");}
+ *(lAB.mp) = *mp | *(lB.mp);
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error");
+ }
+}
+
+void FullBFMatrix::HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
+{
+ try {
+ FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);
+ if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat: Matrices must have same # of rows");}
+ *(lAB.mp) = *mp | B;
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("FullBFMatrix::HorConcat: dynamic cast error");
+ }
+}
+
+void FullBFMatrix::VertConcat(const BFMatrix& B, BFMatrix& AB) const
+{
+ try {
+ const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
+ FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);
+ if (Ncols() != lB.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcat: Matrices must have same # of columns");}
+ *(lAB.mp) = *mp & *(lB.mp);
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error");
+ }
+}
+
+void FullBFMatrix::VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const
+{
+ try {
+ FullBFMatrix& lAB = dynamic_cast<FullBFMatrix&>(AB);
+ if (int(Ncols()) != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcat: Matrices must have same # of columns");}
+ *(lAB.mp) = *mp & B;
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("FullBFMatrix::VertConcat: dynamic cast error");
+ }
+}
+
+//
+// Concatenation of another matrix to *this
+//
+void FullBFMatrix::HorConcat2MyRight(const BFMatrix& B)
+{
+ try {
+ const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
+ if (Nrows() != lB.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
+ *mp |= *(lB.mp);
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: dynamic cast error");
+ }
+}
+
+void FullBFMatrix::HorConcat2MyRight(const NEWMAT::Matrix& B)
+{
+ if (int(Nrows()) != B.Nrows()) {throw BFMatrixException("FullBFMatrix::HorConcat2MyRight: Matrices must have same # of rows");}
+ *mp |= B;
+}
+
+void FullBFMatrix::VertConcatBelowMe(const BFMatrix& B)
+{
+ try {
+ const FullBFMatrix& lB = dynamic_cast<const FullBFMatrix&>(B);
+ if (Ncols() != lB.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
+ *mp &= *(lB.mp);
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("FullBFMatrix::HorConcatBelowMe: dynamic cast error");
+ }
+}
+
+void FullBFMatrix::VertConcatBelowMe(const NEWMAT::Matrix& B)
+{
+ if (int(Ncols()) != B.Ncols()) {throw BFMatrixException("FullBFMatrix::VertConcatBelowMe: Matrices must have same # of columns");}
+ *mp &= B;
+}
+
+// Multiply this matrix with scalar
+
+void FullBFMatrix::MulMeByScalar(double s)
+{
+ *mp = s*(*mp);
+}
+
+// Multiply by vector
+NEWMAT::ReturnMatrix FullBFMatrix::MulByVec(const NEWMAT::ColumnVector& invec) const
+{
+ if (invec.Nrows() != int(Ncols())) {throw BFMatrixException("FullBFMatrix::MulByVec: Matrix-vector size mismatch");}
+ NEWMAT::ColumnVector ret;
+ ret = (*mp)*invec;
+ ret.Release();
+ return(ret);
+}
+
+// Add another matrix to this one
+void FullBFMatrix::AddToMe(const BFMatrix& m, double s)
+{
+ try {
+ const FullBFMatrix& lm = dynamic_cast<const FullBFMatrix&>(m);
+ if (Ncols() != lm.Ncols() || Nrows() != lm.Nrows()) {
+ throw BFMatrixException("FullBFMatrix::AddToMe: Matrix size mismatch");
+ }
+ *mp += s*(*lm.mp);
+ }
+ catch (std::bad_cast) {
+ throw BFMatrixException("FullBFMatrix::AddToMe: dynamic cast error");
+ }
+}
+
+// Given A*x=b, solve for x
+NEWMAT::ReturnMatrix FullBFMatrix::SolveForx(const NEWMAT::ColumnVector& b, // Ignoring all parameters except b
+ MISCMATHS::MatrixType type,
+ double tol,
+ int miter) const
+{
+ if (int(Nrows()) != b.Nrows()) {throw BFMatrixException("FullBFMatrix::SolveForx: Matrix-vector size mismatch");}
+ NEWMAT::ColumnVector ret;
+ ret = mp->i()*b;
+ ret.Release();
+ return(ret);
+}
+
+} // End namespace MISCMATHS
diff --git a/bfmatrix.h b/bfmatrix.h
new file mode 100644
index 0000000..87fe792
--- /dev/null
+++ b/bfmatrix.h
@@ -0,0 +1,270 @@
+// Declarations for class BFMatrix.
+//
+// The purpose of class BFmatrix is to have a class from which
+// to derive 2 other classes; FullBFMatrix and SparseBFMatrix.
+// The reason for this is that the two classes SplineField and
+// DCTField will return Hessian matrices that are either Sparse
+// (SplineField) or full (DCTField). By defining a pure virtual
+// class BFMatrix with a minimal (only what is needed for non-
+// linear reg.) functionality I will be able to write code that
+// is independent of type of matrix returned, and hence of type
+// field.
+//
+// The syntax for the (little) functionality is sort of a mixture
+// of Newmat and SparseMatrix. Mostly SparseMatrix actually.
+// I hope this will not complicate the use of the nonlin package
+// for those who are only interested in the full (normal) case.
+//
+// At one point SparseMatrix was replaced by SpMat as the underlying
+// sparse matrix representation in SparseBFMatrix. SpMat was written
+// with an API that largely mimicks that of NEWMAT. This is the
+// "historical" reason why a wrapper class was written, rather than
+// using templatisation which would have been possible given the
+// similarities in API between SpMat and NEWMAT.
+//
+
+#ifndef BFMatrix_h
+#define BFMatrix_h
+
+#include <boost/shared_ptr.hpp>
+#include "newmat.h"
+#include "SpMat.h"
+#include "cg.h"
+#include "bicg.h"
+
+namespace MISCMATHS {
+
+
+class BFMatrixException: public std::exception
+{
+private:
+ std::string m_msg;
+public:
+ BFMatrixException(const std::string& msg) throw(): m_msg(msg) {}
+
+ virtual const char * what() const throw() {
+ return string("BFMatrix::" + m_msg).c_str();
+ }
+
+ ~BFMatrixException() throw() {}
+};
+
+class BFMatrix
+{
+protected:
+virtual void print(const NEWMAT::Matrix& m, const std::string& fname) const;
+public:
+ // Constructors, destructors and stuff
+ BFMatrix() {}
+ BFMatrix(unsigned int m, unsigned int n) {}
+ virtual ~BFMatrix() {}
+
+ // Access as NEWMAT::Matrix
+ virtual NEWMAT::ReturnMatrix AsMatrix() const = 0;
+
+ // Basic properties
+ virtual unsigned int Nrows() const = 0;
+ virtual unsigned int Ncols() const = 0;
+
+ // Print matrix (for debugging)
+ 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 NEWMAT::Matrix& M) = 0;
+ virtual void Clear() = 0;
+ virtual void Resize(unsigned int m, unsigned int n) = 0;
+
+ // Accessing
+ inline double operator()(unsigned int r, unsigned int c) const {return(Peek(r,c));}
+ virtual double Peek(unsigned int r, unsigned int c) const = 0;
+
+ // Assigning
+ virtual void Set(unsigned int x, unsigned int y, double val) = 0;
+ virtual void Insert(unsigned int x, unsigned int y, double val) = 0;
+ virtual void AddTo(unsigned int x, unsigned int y, double val) = 0;
+
+ // Transpose
+ virtual boost::shared_ptr<BFMatrix> Transpose(boost::shared_ptr<BFMatrix>& pA) const = 0;
+
+ // Concatenation. Note that desired polymorphism prevents us from using BFMatrix->NEWMAT::Matrix conversion
+ // Concatenate two matrices yielding a third
+ // AB = [*this B] in Matlab lingo
+ virtual void HorConcat(const BFMatrix& B, BFMatrix& AB) const = 0;
+ virtual void HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const = 0;
+ // AB = [*this; B] in Matlab lingo
+ virtual void VertConcat(const BFMatrix& B, BFMatrix& AB) const = 0;
+ virtual void VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const = 0;
+ // Concatenate another matrix to *this
+ virtual void HorConcat2MyRight(const BFMatrix& B) = 0;
+ virtual void HorConcat2MyRight(const NEWMAT::Matrix& B) = 0;
+ virtual void VertConcatBelowMe(const BFMatrix& B) = 0;
+ virtual void VertConcatBelowMe(const NEWMAT::Matrix& B) = 0;
+
+ // Multiply by scalar
+ virtual void MulMeByScalar(double s) = 0;
+ // Multiply by vector
+ virtual NEWMAT::ReturnMatrix MulByVec(const NEWMAT::ColumnVector& v) const = 0;
+ // Add another matrix to this one
+ virtual void AddToMe(const BFMatrix& m, double s=1.0) = 0;
+ // Given A*x=b, solve for x.
+ virtual NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector& b,
+ MISCMATHS::MatrixType type=SYM_POSDEF,
+ double tol=1e-6,
+ int miter=200) const = 0;
+};
+
+class SparseBFMatrix : public BFMatrix
+{
+private:
+ boost::shared_ptr<MISCMATHS::SpMat<double> > mp;
+
+public:
+ // Constructors, destructor and assignment
+ SparseBFMatrix()
+ : mp(boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>())) {}
+ SparseBFMatrix(unsigned int m, unsigned int n)
+ : mp(boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(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))) {}
+ SparseBFMatrix(const NEWMAT::Matrix& M)
+ : mp(boost::shared_ptr<MISCMATHS::SpMat<double> >(new MISCMATHS::SpMat<double>(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);
+ }
+
+ // Access as NEWMAT::Matrix
+ virtual NEWMAT::ReturnMatrix AsMatrix() const {NEWMAT::Matrix ret; ret = mp->AsNEWMAT(); ret.Release(); return(ret);}
+
+ // Basic properties
+ virtual unsigned int Nrows() const {return(mp->Nrows());}
+ virtual unsigned int Ncols() const {return(mp->Ncols());}
+
+ // Print matrix (for debugging)
+ virtual void Print(const std::string fname=std::string("")) const {print(mp->AsNEWMAT(),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));}
+
+ // Accessing values
+ virtual double Peek(unsigned int r, unsigned int c) const {return(mp->Peek(r,c));}
+
+ // Setting and inserting values
+ virtual void Set(unsigned int x, unsigned int y, double val) {mp->Set(x,y,val);}
+ virtual void Insert(unsigned int x, unsigned int y, double val) {mp->Set(x,y,val);}
+ virtual void AddTo(unsigned int x, unsigned int y, double val) {mp->AddTo(x,y,val);}
+
+ // Transpose.
+ virtual boost::shared_ptr<BFMatrix> Transpose(boost::shared_ptr<BFMatrix>& pA) const {throw BFMatrixException("SparseBFMatrix::Transpose: Not yet implemented");}
+
+ // Concatenation of two matrices returning a third
+ // AB = [*this B] in Matlab lingo
+ virtual void HorConcat(const BFMatrix& B, BFMatrix& AB) const;
+ virtual void HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const;
+ // AB = [*this; B] in Matlab lingo
+ virtual void VertConcat(const BFMatrix& B, BFMatrix& AB) const;
+ virtual void VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const;
+
+ // Concatenation of another matrix to *this
+ virtual void HorConcat2MyRight(const BFMatrix& B);
+ virtual void HorConcat2MyRight(const NEWMAT::Matrix& B);
+ virtual void VertConcatBelowMe(const BFMatrix& B);
+ virtual void VertConcatBelowMe(const NEWMAT::Matrix& B);
+
+ // Multiply by scalar
+ virtual void MulMeByScalar(double s) {(*mp)*=s;}
+
+ // Multiply by vector
+ virtual NEWMAT::ReturnMatrix MulByVec(const NEWMAT::ColumnVector& invec) const;
+
+ // Add another matrix to this one
+ virtual void AddToMe(const BFMatrix& m, double s=1.0);
+
+ // Given A*x=b, solve for x
+ virtual NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector& b,
+ MISCMATHS::MatrixType type,
+ double tol,
+ int miter) const;
+};
+
+class FullBFMatrix : public BFMatrix
+{
+private:
+ boost::shared_ptr<NEWMAT::Matrix> mp;
+public:
+ // Constructors, destructor and assignment
+ FullBFMatrix() {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix());}
+ FullBFMatrix(unsigned int m, unsigned int n) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(m,n));}
+ FullBFMatrix(const MISCMATHS::SpMat<double>& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M.AsNEWMAT()));}
+ FullBFMatrix(const NEWMAT::Matrix& M) {mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(M));}
+ virtual ~FullBFMatrix() {}
+ virtual const FullBFMatrix& operator=(const FullBFMatrix& M) {
+ mp = boost::shared_ptr<NEWMAT::Matrix>(new NEWMAT::Matrix(*(M.mp))); return(*this);
+ }
+
+ virtual NEWMAT::ReturnMatrix AsMatrix() const {NEWMAT::Matrix ret; ret = *mp; ret.Release(); return(ret);}
+
+ // Basic properties
+ virtual unsigned int Nrows() const {return(mp->Nrows());}
+ virtual unsigned int Ncols() const {return(mp->Ncols());}
+
+ // Print matrix (for debugging)
+ virtual void Print(const std::string fname=std::string("")) const {print(*mp,fname);}
+
+ // 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 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);}
+ virtual void Resize(unsigned int m, unsigned int n) {mp->ReSize(m,n);}
+
+ // Accessing values
+ virtual double Peek(unsigned int r, unsigned int c) const {return((*mp)(r,c));}
+
+ // Setting and inserting values.
+ virtual void Set(unsigned int x, unsigned int y, double val) {(*mp)(x,y)=val;}
+ virtual void Insert(unsigned int x, unsigned int y, double val) {(*mp)(x,y)=val;}
+ virtual void AddTo(unsigned int x, unsigned int y, double val) {(*mp)(x,y)+=val;}
+
+ // Transpose.
+ virtual boost::shared_ptr<BFMatrix> Transpose(boost::shared_ptr<BFMatrix>& pA) const;
+
+ // Concatenation of two matrices returning a third
+ virtual void HorConcat(const BFMatrix& B, BFMatrix& AB) const;
+ virtual void HorConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const;
+ virtual void VertConcat(const BFMatrix& B, BFMatrix& AB) const;
+ virtual void VertConcat(const NEWMAT::Matrix& B, BFMatrix& AB) const;
+
+ // Concatenation of another matrix to *this
+ virtual void HorConcat2MyRight(const BFMatrix& B);
+ virtual void HorConcat2MyRight(const NEWMAT::Matrix& B);
+ virtual void VertConcatBelowMe(const BFMatrix& B);
+ virtual void VertConcatBelowMe(const NEWMAT::Matrix& B);
+
+ // Multiply by scalar
+ virtual void MulMeByScalar(double s);
+
+ // Multiply by vector
+ virtual NEWMAT::ReturnMatrix MulByVec(const NEWMAT::ColumnVector& invec) const;
+
+ // Add another matrix to this one
+ virtual void AddToMe(const BFMatrix& m, double s);
+
+ // Given A*x=b, solve for x
+ virtual NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector& b,
+ MISCMATHS::MatrixType type,
+ double tol,
+ int miter) const;
+
+};
+
+} // End namespace MISCMATHS
+
+#endif // End #ifndef BFMatrix_h
--
GitLab