From 305067486e2cf831fc2edfa05960a4ed28936f07 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sun, 14 Apr 2013 12:14:49 +0200 Subject: [PATCH] Improve Doxygen documentation, add at() methods and value assignment operator. --- include/GMatrixBase.hpp | 143 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------- include/GMatrixSparse.hpp | 101 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------ pyext/GMatrixBase.i | 40 +++++++++++++++++++++------------------- pyext/GMatrixSparse.i | 14 ++++++++++---- src/linalg/GMatrixSparse.cpp | 212 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------ 5 files changed, 389 insertions(+), 121 deletions(-) diff --git a/include/GMatrixBase.hpp b/include/GMatrixBase.hpp index 26ceff5..58094f7 100644 --- a/include/GMatrixBase.hpp +++ b/include/GMatrixBase.hpp @@ -38,9 +38,110 @@ * * @brief Abstract matrix base class interface defintion * - * This class is an abstract base class from which all matrix classes will - * be derived. It defines the common interface for all matrix classes and - * implements common methods and data members. + * This class provides the abstract interface for all matrix classes. A + * GammaLib matrix is a two dimensional array of double precision values + * that is commonly used for linear algebra computations. + * + * The matrix classes were designed in a manner so that matrices can be used + * in a natural syntax that is comparable to the syntax used for scalar + * values. In particular, the following matrix operations are available: + * + * assigned_matrix = matrix; + * assigned_matrix = value; + * added_matrix = matrix1 + matrix2; + * added_matrix += matrix; + * subtracted_matrix = matrix1 - matrix2; + * subtracted_matrix -= matrix; + * multiplied_matrix = matrix1 * matrix2; + * multiplied_matrix *= matrix; + * scaled_matrix = matrix * scale; + * scaled_matrix = scale * matrix; + * scaled_matrix *= scale; + * divided_matrix = matrix / factor; + * divided_matrix /= factor; + * negated_matrix = -matrix; + * + * where @p matrix, @p matrix1 and @p matrix2 are matrices, and + * @p value, @p scale and @p factor are double precision values. For the + * value assignment operator, each matrix element is assigned the specified + * value. For the scaling and division operators, every matrix element is + * multiplied or divided by the specified value. + * + * Matrices can also be compared using + * + * if (matrix == another_matrix) ... + * if (matrix != another_matrix) ... + * + * where identify is defined as having the same number of rows and columns + * and every element having the same value. + * + * Matrix elements are accessed using the access operators or the at() methods + * + * double element = matrix(row, column); + * double element = matrix.at(row, column); + * + * the difference being that the access operators do NOT but the at() + * methods do perform validity checking of the @p row and @p column + * arguments. Access operators and at() methods exist for non-const and + * const references. + * + * For more efficient matrix access, methods operating on matrix rows and + * columns are available. The + * + * GVector row = matrix.row(index); + * GVector column = matrix.column(index); + * + * methods return a row or a column in a GVector object, while the + * + * GVector row; + * GVector column; + * ... + * matrix.row(index, row); + * matrix.column(index, column); + * + * store a row or a column provided by a GVector object into the matrix. + * Obviously, the vector needs to be of the same size as the number of rows + * or columns available in the matrix. Adding of matrix elements row-by-row + * or column-by-column is done using the + * + * matrix.add_to_row(index, row); + * matrix.add_to_column(index, column); + * + * methods. + * + * Matrix multiplication with a vector is done using + * + * GVector vector; + * ... + * GVector product = matrix * vector; + * + * where the size of @p vector needs to be identical to the number of + * columns in the matrix. + * + * Furthermore, the following methods exist for matrix operations: + * + * matrix.transpose(); // Transposes matrix + * matrix.invert(); // Inverts matrix + * matrix.negate(); // Negates all matrix elements + * matrix.abs(); // Sets all matrix elements to their absolute values + * + * Friend functions are implemented to perform the same operations also in a + * more natural syntax that will generate a new matrix + * + * transposed_matrix = transpose(matrix); + * inverted_matrix = invert(matrix); + * negated_matrix = -matrix; + * absolute_matrix = abs(matrix); + * + * The following methods allow to access matrix attributes: + * + * matrix.fill(); // Percentage of non-zero elements in matrix [0,1] + * matrix.min(); // Minimum matrix element + * matrix.max(); // Maximum matrix element + * matrix.sum(); // Sum of all matrix elements + * matrix.rows(); // Number of rows in matrix + * matrix.columns(); // Number of columns in matrix + * matrix.size(); // Number of elements in matrix ***************************************************************************/ class GMatrixBase : public GBase { @@ -61,23 +162,25 @@ public: virtual bool operator!=(const GMatrixBase& matrix) const; // Pure virtual methods - virtual void clear(void) = 0; - virtual GMatrixBase* clone(void) const = 0; - virtual GVector row(const int& row) const = 0; - virtual void row(const int& row, const GVector& vector) = 0; - virtual GVector column(const int& column) const = 0; - virtual void column(const int& column, const GVector& vector) = 0; - virtual void add_to_row(const int& row, const GVector& vector) = 0; - virtual void add_to_column(const int& column, const GVector& vector) = 0; - virtual void transpose(void) = 0; - virtual void invert(void) = 0; - virtual void negate(void) = 0; - virtual void abs(void) = 0; - virtual double fill(void) const = 0; - virtual double min(void) const = 0; - virtual double max(void) const = 0; - virtual double sum(void) const = 0; - virtual std::string print(const GChatter& chatter = NORMAL) const = 0; + virtual void clear(void) = 0; + virtual GMatrixBase* clone(void) const = 0; + virtual double& at(const int& row, const int& column) = 0; + virtual const double& at(const int& row, const int& column) const = 0; + virtual GVector row(const int& row) const = 0; + virtual void row(const int& row, const GVector& vector) = 0; + virtual GVector column(const int& column) const = 0; + virtual void column(const int& column, const GVector& vector) = 0; + virtual void add_to_row(const int& row, const GVector& vector) = 0; + virtual void add_to_column(const int& column, const GVector& vector) = 0; + virtual void transpose(void) = 0; + virtual void invert(void) = 0; + virtual void negate(void) = 0; + virtual void abs(void) = 0; + virtual double fill(void) const = 0; + virtual double min(void) const = 0; + virtual double max(void) const = 0; + virtual double sum(void) const = 0; + virtual std::string print(const GChatter& chatter = NORMAL) const = 0; // Base class methods const int& size(void) const; diff --git a/include/GMatrixSparse.hpp b/include/GMatrixSparse.hpp index 218cdb0..34926f0 100644 --- a/include/GMatrixSparse.hpp +++ b/include/GMatrixSparse.hpp @@ -49,16 +49,20 @@ class GSparseNumeric; * * @brief Sparse matrix class interface defintion * - * This class implements a sparse matrix class. The class only stores - * non-zero elements, which can considerably reduce the memory requirements - * for large systems that are sparsly filled. + * This class implements a sparse matrix. The class only stores non-zero + * elements, which can considerably reduce the memory requirements for large + * systems that are sparsely filled. Also the computations will be faster + * because only non-zero elements will be used in the operations. + * + * For a description of the common matrix methods, please refer to the + * GMatrixBase class. * * The class has been inspired from the code CSparse that can be downloaded - * from http://www.cise.ufl.edu/research/sparse/CSparse/ - * CSparse has been written by Timothy A. Davis to accompany his book - * "Direct Methods for Sparse Linear Systems". In the CSparse code, a sparse - * matrix is implemented by the structure cs_sparse. The members of this - * structure map as follows to the member of the GMatrixSparse class: + * from http://www.cise.ufl.edu/research/sparse/CSparse/ which has been + * written by Timothy A. Davis to accompany his book "Direct Methods for + * Sparse Linear Systems". In the CSparse code, a sparse matrix is + * implemented by the structure cs_sparse. The members of this + * structure map as follows to the members of the GMatrixSparse class: * * * @@ -106,6 +110,76 @@ class GSparseNumeric; * *
* + * Matrix allocation is done using the constructors + * + * GMatrixSparse sparsematrix(rows, columns); + * GMatrixSparse sparsematrix(rows, columns, elements); + * GMatrixSparse sparsematrix(matrix); + * GMatrixSparse sparsematrix(sparsematrix); + * GMatrixSparse sparsematrix(symmetricmatrix); + * + * where @p rows and @p columns specify the number of rows and columns of + * the matrix and @p elements is the size of the memory that should be + * pre-allocated for matrix storage (optional). Storage conversion + * constructors exist that allow allocating a sparse matrix by copying from + * a general matrix of type GMatrix and a symmetrix matrix of type + * GMatrixSymmetric. + * + * Because the location of a specific matrix element is not necessarily known + * in advance, the class implements a so called pending element for matrix + * element setting. A pending element is an element that is scheduled for + * insertion in the matrix, but which has not yet been inserted. The matrix + * element access operators () and at() return the reference to this pending + * element if a memory location for a given row and column does not yet + * exist. The protected method fill_pending() is then called before all + * operations that are done on the matrix to insert the element before the + * operation is executed. + * + * Another mechanism that has been implemented to speed up the fill of the + * sparse matrix is a "fill stack" that is used for insertion or addition + * of matrix columns by the + * + * matrix.column(index, vector); + * matrix.column(index, values, rows, number); + * matrix.add_to_column(index, vector); + * matrix.add_to_column(index, values, rows, number); + * + * methods. The fill stack is a buffer that implements a queue for columns + * that are to be inserted into the matrix. This is more efficient than + * insertion of elements one by one, as every insertion of an element may + * require to shift all elements back by one memory location. + * + * The fill stack is used as follows: + * + * matrix.stack_init(size, entries); + * ... + * matrix.column(index, vector); // (or one of the other column methods) + * ... + * matrix.stack_flush(); + * ... + * matrix.stack_destroy(); + * + * The stack_init(size, entries) method initialises the fill stack, where + * @p size is the size of the allocated memory buffer and @p entries is the + * maximum number of entries (i.e. columns) that will be held by the buffer. + * If @p size is set to 0 (the default value), a default @p size value of + * 512 is used. If @p entries is set to 0 (the default value), the number of + * matrix columns is taken as default @p entries value. Note that a too large + * number of elements will produce some overhead due to fill stack + * management, hence @p entries should not exceeed a value of the order of + * 10-100. + * + * The stack_flush() method flushes the stack, which is mandatory before any + * usage of the matrix. Note that the fill stack IS NOT INSERTED AUTOMATICALLY + * before any matrix operation, hence manual stack flushing is needed to make + * all filled matrix elements available for usage. The stack_destroy() method + * will flush the stack and free all stack elements. This method + * should be called once no filling is required anymore. If stack_destroy() + * is called immediately after filling, no call to stack_flush() is needed as + * the stack_destroy() method flushes the stack before destroying it. The + * matrix stack is also destroyed by the sparse matrix destructor, hence + * manual stack destruction is not mandatory. + * * Except for *m_rowinx which is implemented on the level of GMatrixSparse, * all other members are implemented by the base class GMatrixBase. ***************************************************************************/ @@ -144,6 +218,7 @@ public: // Other operators virtual GMatrixSparse& operator=(const GMatrixSparse& matrix); + virtual GMatrixSparse& operator=(const double& value); virtual GMatrixSparse operator+(const GMatrixSparse& matrix) const; virtual GMatrixSparse operator-(const GMatrixSparse& matrix) const; virtual GMatrixSparse operator*(const GMatrixSparse& matrix) const; @@ -157,6 +232,8 @@ public: // Implemented pure virtual base class methods virtual void clear(void); virtual GMatrixSparse* clone(void) const; + virtual double& at(const int& row, const int& column); + virtual const double& at(const int& row, const int& column) const; virtual GVector row(const int& row) const; virtual void row(const int& row, const GVector& vector); virtual GVector column(const int& column) const; @@ -174,10 +251,10 @@ public: virtual std::string print(const GChatter& chatter = NORMAL) const; // Other methods - void column(const double* values, const int* rows, - int number, const int& col); - void add_to_column(const double* values, const int* rows, - int number, const int& col); + void column(const int& column, const double* values, + const int* rows, int number); + void add_to_column(const int& column, const double* values, + const int* rows, int number); void cholesky_decompose(bool compress = true); GVector cholesky_solver(const GVector& vector, bool compress = true); void cholesky_invert(bool compress = true); diff --git a/pyext/GMatrixBase.i b/pyext/GMatrixBase.i index 7bea5f3..689de5e 100755 --- a/pyext/GMatrixBase.i +++ b/pyext/GMatrixBase.i @@ -42,29 +42,31 @@ public: virtual ~GMatrixBase(void); // Pure virtual operators - virtual GVector operator*(const GVector& vector) const = 0; + virtual GVector operator*(const GVector& vector) const = 0; // Implemented base class operators - virtual bool operator==(const GMatrixBase& matrix) const; - virtual bool operator!=(const GMatrixBase& matrix) const; + virtual bool operator==(const GMatrixBase& matrix) const; + virtual bool operator!=(const GMatrixBase& matrix) const; // Pure virtual methods - virtual void clear(void) = 0; - virtual GMatrixBase* clone(void) const = 0; - virtual GVector row(const int& row) const = 0; - virtual void row(const int& row, const GVector& vector) = 0; - virtual GVector column(const int& column) const = 0; - virtual void column(const int& column, const GVector& vector) = 0; - virtual void add_to_row(const int& row, const GVector& vector) = 0; - virtual void add_to_column(const int& column, const GVector& vector) = 0; - virtual void transpose(void) = 0; - virtual void invert(void) = 0; - virtual void negate(void) = 0; - virtual void abs(void) = 0; - virtual double fill(void) const = 0; - virtual double min(void) const = 0; - virtual double max(void) const = 0; - virtual double sum(void) const = 0; + virtual void clear(void) = 0; + virtual GMatrixBase* clone(void) const = 0; + virtual double& at(const int& row, const int& column) = 0; + virtual const double& at(const int& row, const int& column) const = 0; + virtual GVector row(const int& row) const = 0; + virtual void row(const int& row, const GVector& vector) = 0; + virtual GVector column(const int& column) const = 0; + virtual void column(const int& column, const GVector& vector) = 0; + virtual void add_to_row(const int& row, const GVector& vector) = 0; + virtual void add_to_column(const int& column, const GVector& vector) = 0; + virtual void transpose(void) = 0; + virtual void invert(void) = 0; + virtual void negate(void) = 0; + virtual void abs(void) = 0; + virtual double fill(void) const = 0; + virtual double min(void) const = 0; + virtual double max(void) const = 0; + virtual double sum(void) const = 0; // Base class methods const int& size(void) const; diff --git a/pyext/GMatrixSparse.i b/pyext/GMatrixSparse.i index 706827e..b16374e 100644 --- a/pyext/GMatrixSparse.i +++ b/pyext/GMatrixSparse.i @@ -73,6 +73,8 @@ public: // Implemented pure virtual base class methods virtual void clear(void); virtual GMatrixSparse* clone(void) const; + virtual double& at(const int& row, const int& column); + virtual const double& at(const int& row, const int& column) const; virtual GVector row(const int& row) const; virtual void row(const int& row, const GVector& vector); virtual GVector column(const int& column) const; @@ -89,10 +91,10 @@ public: virtual double sum(void) const; // Other methods - void add_to_column(const double* values, const int* rows, - int number, const int& column); - void column(const double* values, const int* rows, - int number, const int& column); + void column(const int& column, const double* values, + const int* rows, int number); + void add_to_column(const int& column, const double* values, + const int* rows, int number); void cholesky_decompose(bool compress = true); GVector cholesky_solver(const GVector& vector, bool compress = true); void cholesky_invert(bool compress = true); @@ -125,4 +127,8 @@ public: GMatrixSparse copy() { return (*self); } + GMatrixSparse set(const double& value) { + (*self) = value; + return (*self); + } }; diff --git a/src/linalg/GMatrixSparse.cpp b/src/linalg/GMatrixSparse.cpp index c850cc7..438a9dc 100644 --- a/src/linalg/GMatrixSparse.cpp +++ b/src/linalg/GMatrixSparse.cpp @@ -17,19 +17,6 @@ * You should have received a copy of the GNU General Public License * * along with this program. If not, see . * * * - * ----------------------------------------------------------------------- * - * This class implements the compressed sparse column format. The * - * following arrays are allocated: * - * * - * m_rows (n) Number of rows * - * m_cols (m) Number of columns * - * m_data (Ax) Holds all 'elements' non-zero values, in column order * - * m_colstart (Ap) Holds the index of the first element of each column * - * m_rowinx (Ai) Holds the row indices for all elements * - * * - * Column 'col' covers therefore [m_colstart[j], ..., m_colstart[j+1]-1] * - * Note that 'm_colstart' has m_cols+1 elements. * - * * ***************************************************************************/ /** * @file GMatrixSparse.cpp @@ -53,20 +40,20 @@ /* __ Method name definitions ____________________________________________ */ #define G_CONSTRUCTOR "GMatrixSparse::GMatrixSparse(int&, int&, int&)" -#define G_ACCESS "GMatrixSparse::operator(int&, int&)" #define G_OP_MUL_VEC "GMatrixSparse::operator*(GVector&)" #define G_OP_ADD "GMatrixSparse::operator+=(GMatrixSparse&)" #define G_OP_SUB "GMatrixSparse::operator-=(GMatrixSparse&)" #define G_OP_MAT_MUL "GMatrixSparse::operator*=(GMatrixSparse&)" +#define G_AT "GMatrixSparse::at(int&, int&)" #define G_EXTRACT_ROW "GMatrixSymmetric::row(int&)" #define G_SET_ROW "GMatrixSymmetric::row(int&, GVector&)" #define G_EXTRACT_COLUMN "GMatrixSymmetric::column(int&)" #define G_SET_COLUMN "GMatrixSymmetric::column(int&, GVector&)" -#define G_SET_COLUMN2 "GMatrixSymmetric::column(double*, int*, int, int&)" +#define G_SET_COLUMN2 "GMatrixSymmetric::column(int&, double*, int*, int)" #define G_ADD_TO_ROW "GMatrixSymmetric::add_to_row(int&, GVector&)" #define G_ADD_TO_COLUMN "GMatrixSymmetric::add_to_column(int&, GVector&)" -#define G_ADD_TO_COLUMN2 "GMatrixSymmetric::add_to_column(double*, int*,"\ - " int, int&)" +#define G_ADD_TO_COLUMN2 "GMatrixSymmetric::add_to_column(int&, double*," + " int*, int)" #define G_INVERT "GMatrixSparse::invert(void)" #define G_CHOL_DECOMP "GMatrixSparse::cholesky_decompose(bool)" #define G_CHOL_SOLVE "GMatrixSparse::cholesky_solver(GVector&, bool)" @@ -137,6 +124,9 @@ GMatrixSparse::GMatrixSparse(void) : GMatrixBase() * of rows and columns. The elements parameter allows to specify how much * physical memory should be allocated initially. By default, no memory * will be allocated. + * + * @todo Is there a real need to throw an empty exception? The class should + * be able to operate on empty matrices. ***************************************************************************/ GMatrixSparse::GMatrixSparse(const int& rows, const int& columns, @@ -257,9 +247,12 @@ GMatrixSparse::~GMatrixSparse(void) ==========================================================================*/ /***********************************************************************//** - * @brief Assignment operator + * @brief Matrix assignment operator * * @param[in] matrix Matrix. + * @return Matrix. + * + * Assigns the content of another matrix to the actual matrix instance. ***************************************************************************/ GMatrixSparse& GMatrixSparse::operator=(const GMatrixSparse& matrix) { @@ -288,25 +281,65 @@ GMatrixSparse& GMatrixSparse::operator=(const GMatrixSparse& matrix) /***********************************************************************//** - * @brief Access operator + * @brief Value assignment operator + * + * @param[in] value Value. + * @return Matrix. + * + * Assigns the specified @p value to all elements of the matrix. + ***************************************************************************/ +GMatrixSparse& GMatrixSparse::operator=(const double& value) +{ + // Fill any pending element to have a non-pending state + fill_pending(); + + // If value is 0 then simply reinitialize column start indices + if (value == 0) { + + // Initialise column start indices to 0 + for (int col = 0; col <= m_cols; ++col) { + m_colstart[col] = 0; + } + + } + + // ... otherwise fill column-wise + else { + + // Set column vector + GVector column(m_rows); + column = value; + + // Column-wise setting + for (int col = 0; col < m_cols; ++col) { + this->column(col, column); + } + + } + + // Return this object + return *this; +} + + +/***********************************************************************//** + * @brief Return reference to matrix element * * @param[in] row Matrix row [0,...,rows()-1]. * @param[in] column Matrix column [0,...,columns()-1]. + * @return Reference to matrix element. * - * @exception GException::out_of_range - * Row or column index out of range. + * Returns the reference to the matrix element at @p row and @p column. If + * the matrix element does not yet exist, a reference to the pending element + * with a value of 0.0 is returned. ***************************************************************************/ double& GMatrixSparse::operator()(const int& row, const int& column) { - // Compile option: perform range check - #if defined(G_RANGE_CHECK) - if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) { - throw GException::out_of_range(G_ACCESS, row, column, m_rows, m_cols); - } - #endif + // Fill pending element. This will set the value of the pending element + // to 0.0. + fill_pending(); // Get element - fill_pending(); int inx = get_index(row, column); double* value; if (inx < 0) { @@ -324,29 +357,24 @@ double& GMatrixSparse::operator()(const int& row, const int& column) /***********************************************************************//** - * @brief Access operator (const version) + * @brief Return reference to matrix element (const version) * * @param[in] row Matrix row [0,...,rows()-1]. * @param[in] column Matrix column [0,...,columns()-1]. + * @return Const reference to matrix element. * - * @exception GException::out_of_range - * Row or column index out of range. - * - * We need here the zero element to return also a pointer for 0.0 entry that - * is not stored. Since we have the const version we don't have to care about - * modification of this zero value. + * Returns a const reference to the matrix element at @p row and @p column. + * If the matrix element does not yet exist, a reference to the zero element + * is returned. If the matrix element corresponds to the pending element, + * a reference to the pending element is returned. Otherwise, a reference + * to the matrix elements is returned. ***************************************************************************/ const double& GMatrixSparse::operator()(const int& row, const int& column) const { - // Compile option: perform range check - #if defined(G_RANGE_CHECK) - if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) { - throw GException::out_of_range(G_ACCESS, row, column, m_rows, m_cols); - } - #endif - - // Get element + // Get element. We need here the zero element to return also a pointer + // for 0.0 entry that is not stored. Since we have the const version we + // don't have to care about modification of this zero value. int inx = get_index(row, column); double* value; if (inx < 0) { @@ -635,6 +663,50 @@ GMatrixSparse* GMatrixSparse::clone(void) const /***********************************************************************//** + * @brief Return reference to matrix element + * + * @param[in] row Matrix row [0,...,rows()-1]. + * @param[in] column Matrix column [0,...,columns()-1]. + * @return Reference to matrix element. + * + * @exception GException::out_of_range + * Row or column index out of range. + ***************************************************************************/ +double& GMatrixSparse::at(const int& row, const int& column) +{ + // Raise exception if row or column index is out of range + if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) { + throw GException::out_of_range(G_AT, row, column, m_rows, m_cols); + } + + // Return element + return ((*this)(row, column)); +} + + +/***********************************************************************//** + * @brief Return reference to matrix element (const version) + * + * @param[in] row Matrix row [0,...,rows()-1]. + * @param[in] column Matrix column [0,...,columns()-1]. + * @return Reference to matrix element. + * + * @exception GException::out_of_range + * Row or column index out of range. + ***************************************************************************/ +const double& GMatrixSparse::at(const int& row, const int& column) const +{ + // Raise exception if row or column index is out of range + if (row < 0 || row >= m_rows || column < 0 || column >= m_cols) { + throw GException::out_of_range(G_AT, row, column, m_rows, m_cols); + } + + // Return element + return ((*this)(row, column)); +} + + +/***********************************************************************//** * @brief Extract row as vector from matrix * * @param[in] row Row to be extracted (starting from 0). @@ -888,10 +960,10 @@ void GMatrixSparse::column(const int& column, const GVector& vector) /***********************************************************************//** * @brief Insert compressed array into matrix column * + * @param[in] column Column index [0,...,columns()-1]. * @param[in] values Compressed array. * @param[in] rows Row indices of array. * @param[in] number Number of elements in array. - * @param[in] column Column index [0,...,columns()-1]. * * @exception GException::out_of_range * Invalid column index specified. @@ -904,12 +976,13 @@ void GMatrixSparse::column(const int& column, const GVector& vector) * This is the main driver routine to insert data into a matrix. Note that * there is another instance of this function that takes a vector. ***************************************************************************/ -void GMatrixSparse::column(const double* values, const int* rows, - int number, const int& column) +void GMatrixSparse::column(const int& column, const double* values, + const int* rows, int number) { // Debug header #if defined(G_DEBUG_SPARSE_INSERTION) - std::cout << "GMatrixSparse::insert_col(v, i, n, " << column << "):" << std::endl; + std::cout << "GMatrixSparse::column("; + std::cout << column << ", values, rows, " << number << "):" << std::endl; std::cout << " Matrix Data : "; for (int i = 0; i < m_elements; ++i) { std::cout << m_data[i] << " "; @@ -1166,10 +1239,10 @@ void GMatrixSparse::add_to_column(const int& column, const GVector& vector) /***********************************************************************//** * @brief Add compressed array into matrix column * + * @param[in] column Column index [0,...,columns()-1]. * @param[in] values Compressed array. * @param[in] rows Row indices of array. * @param[in] number Number of elements in array. - * @param[in] col Column index (starting from 0). * * @exception GException::out_of_range * Invalid column index specified. @@ -1182,12 +1255,13 @@ void GMatrixSparse::add_to_column(const int& column, const GVector& vector) * normal and stack-based filled. Note that there is another instance of this * method that takes a vector. ***************************************************************************/ -void GMatrixSparse::add_to_column(const double* values, const int* rows, - int number, const int& col) +void GMatrixSparse::add_to_column(const int& column, const double* values, + const int* rows, int number) { // Debug header #if defined(G_DEBUG_SPARSE_ADDITION) - std::cout << "GMatrixSparse::add_col(v, i, n, " << col << "):" << std::endl; + std::cout << "GMatrixSparse::add_col("; + std::cout << column << ", values, rows, " << number << "):" << std::endl; std::cout << " Matrix Data : "; for (int i = 0; i < m_elements; ++i) { std::cout << m_data[i] << " "; @@ -1216,7 +1290,7 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows, // stack_push_column does its own argument verifications, so to avoid // double checking we don't do anything before this call ... if (m_stack_data != NULL) { - number = stack_push_column(values, rows, number, col); + number = stack_push_column(values, rows, number, column); if (number == 0) { return; } @@ -1231,8 +1305,8 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows, // Raise an exception if the column index is invalid #if defined(G_RANGE_CHECK) - if (col < 0 || col >= m_cols) { - throw GException::out_of_range(G_ADD_TO_COLUMN2, col, 0, m_cols-1); + if (column < 0 || column >= m_cols) { + throw GException::out_of_range(G_ADD_TO_COLUMN2, column, 0, m_cols-1); } #endif @@ -1245,8 +1319,8 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows, } // endelse: there was no stack // Get indices of column in matrix - int i_start = m_colstart[col]; - int i_stop = m_colstart[col+1]; + int i_start = m_colstart[column]; + int i_stop = m_colstart[column+1]; // Case A: the column exists in the matrix, so mix new elements with existing // data @@ -1267,7 +1341,7 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows, wrk_double, wrk_int, &num_mix); // Insert mixed column - this->column(wrk_double, wrk_int, num_mix, col); + this->column(column, wrk_double, wrk_int, num_mix); // Free workspace delete [] wrk_int; @@ -1277,20 +1351,23 @@ void GMatrixSparse::add_to_column(const double* values, const int* rows, // Case B: the column does not yet exist in the matrix, so just insert it else { - this->column(values, rows, number, col); + this->column(column, values, rows, number); } // Debugging: show sparse matrix after insertion #if defined(G_DEBUG_SPARSE_ADDITION) std::cout << " Out Data: "; - for (int i = 0; i < m_elements; ++i) + for (int i = 0; i < m_elements; ++i) { std::cout << m_data[i] << " "; + } std::cout << std::endl << " Out Row : "; - for (int i = 0; i < m_elements; ++i) + for (int i = 0; i < m_elements; ++i) { std::cout << m_rowinx[i] << " "; + } std::cout << std::endl << " Out Col : "; - for (int i = 0; i < m_cols+1; ++i) + for (int i = 0; i < m_cols+1; ++i) { std::cout << m_colstart[i] << " "; + } std::cout << std::endl; #endif @@ -1321,15 +1398,18 @@ void GMatrixSparse::transpose(void) /***********************************************************************//** * @brief Invert matrix * - * @exception GException::feature_not_implemented - * Feature not yet implemented. + * Inverts the matrix using the Cholesky decomposition. This does not work + * on any kind of matrix. * - * @todo Needs to be implemented. + * @todo Specify in documentation for which kind of matrix the method works. ***************************************************************************/ void GMatrixSparse::invert(void) { - // Throw exception - throw GException::feature_not_implemented(G_INVERT); + // Perform Chloesky decomposition + cholesky_decompose(true); + + // Invert matrix + cholesky_invert(true); // Return return; -- libgit2 0.21.4