Francois Brun / gammalib

Commit 150a17e86e4433e8cdbc4614c51a2e754b450673

Authored by knodlseder
1 parent 14f661df

start conform to coding conventions

Inline Side-by-side
Showing 1 changed file with 670 additions and 688 deletions
src/linalg/GSparseMatrix.cpp
  View file @150a17e
1 /*************************************************************************** 1 /***************************************************************************
2 * GSparseMatrix.cpp - sparse matrix class * 2 * GSparseMatrix.cpp - sparse matrix class *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 - * copyright (C) 2006-2009 by Jurgen Knodlseder * 4 + * copyright (C) 2006-2010 by Jurgen Knodlseder *
5 * ----------------------------------------------------------------------- * 5 * ----------------------------------------------------------------------- *
6 * * 6 * *
7 * This program is free software; you can redistribute it and/or modify * 7 * This program is free software; you can redistribute it and/or modify *
@@ -18,42 +18,41 @@ @@ -18,42 +18,41 @@
18 * Column 'col' covers therefore [m_colstart[j], ..., m_colstart[j+1]-1] * 18 * Column 'col' covers therefore [m_colstart[j], ..., m_colstart[j+1]-1] *
19 * Note that 'm_colstart' has m_cols+1 elements. * 19 * Note that 'm_colstart' has m_cols+1 elements. *
20 ***************************************************************************/ 20 ***************************************************************************/
  21 +/**
  22 + * @file GSparseMatrix.cpp
  23 + * @brief GSparseMatrix class implementation
  24 + * @author J. Knodlseder
  25 + */
21 26
22 /* __ Includes ___________________________________________________________ */ 27 /* __ Includes ___________________________________________________________ */
23 -#include "GException.hpp"  
24 #include "GVector.hpp" 28 #include "GVector.hpp"
25 -#include "GMatrix.hpp"  
26 -#include "GSymMatrix.hpp"  
27 #include "GSparseMatrix.hpp" 29 #include "GSparseMatrix.hpp"
28 #include "GSparseSymbolic.hpp" 30 #include "GSparseSymbolic.hpp"
29 #include "GSparseNumeric.hpp" 31 #include "GSparseNumeric.hpp"
30 32
31 -/* __ Namespaces _________________________________________________________ */  
32 -using namespace std;  
33 -  
34 /* __ Method name definitions ____________________________________________ */ 33 /* __ Method name definitions ____________________________________________ */
35 -#define G_OP_MUL_VEC "GSparseMatrix::operator* (const GVector&) const"  
36 -#define G_OP_ADD "GSparseMatrix::operator+= (const GSparseMatrix&)"  
37 -#define G_OP_SUB "GSparseMatrix::operator-= (const GSparseMatrix&)"  
38 -#define G_OP_MAT_MUL "GSparseMatrix::operator*= (const GSparseMatrix&)"  
39 -#define G_ADD_COL "GSparseMatrix::add_col(const GVector&, int)"  
40 -#define G_ADD_COL2 "GSparseMatrix::add_col(const double*, const int*, int, int)"  
41 -#define G_CHOL_DECOMP "GSparseMatrix::cholesky_decompose(int)"  
42 -#define G_CHOL_SOLVE "GSparseMatrix::cholesky_solver(const GVector&, int)"  
43 -#define G_EXTRACT_ROW "GSparseMatrix::extract_row(int) const"  
44 -#define G_EXTRACT_COL "GSparseMatrix::extract_col(int) const"  
45 -#define G_INSERT_COL "GSparseMatrix::insert_col(const GVector&, int)"  
46 -#define G_INSERT_COL2 "GSparseMatrix::insert_col(const double*, const int*, int, int)"  
47 -#define G_STACK_INIT "GSparseMatrix::stack_init(int, int)"  
48 -#define G_STACK_PUSH "GSparseMatrix::stack_push_column(const double*, const int*, int, int)"  
49 -#define G_STACK_FLUSH "GSparseMatrix::stack_flush()"  
50 -#define G_CONSTRUCTOR "GSparseMatrix::constructor(int, int, int)"  
51 -#define G_COPY_MEMBERS "GSparseMatrix::copy_members(const GSparseMatrix&)"  
52 -#define G_ALLOC "GSparseMatrix::alloc_elements(int, int)"  
53 -#define G_FREE "GSparseMatrix::free_elements(int, int)"  
54 -#define G_REMOVE_ZERO "GSparseMatrix::remove_zero_row_col()"  
55 -#define G_SYMPERM "cs_symperm(GSparseMatrix*, const int*, int)"  
56 -#define G_TRANSPOSE "cs_transpose(GSparseMatrix*, int)" 34 +#define G_OP_MUL_VEC "GSparseMatrix::operator* (GVector)"
  35 +#define G_OP_ADD "GSparseMatrix::operator+= (GSparseMatrix)"
  36 +#define G_OP_SUB "GSparseMatrix::operator-= (GSparseMatrix)"
  37 +#define G_OP_MAT_MUL "GSparseMatrix::operator*= (GSparseMatrix)"
  38 +#define G_ADD_COL "GSparseMatrix::add_col(GVector,int)"
  39 +#define G_ADD_COL2 "GSparseMatrix::add_col(double*,int*int,int)"
  40 +#define G_CHOL_DECOMP "GSparseMatrix::cholesky_decompose(int)"
  41 +#define G_CHOL_SOLVE "GSparseMatrix::cholesky_solver(GVector,int)"
  42 +#define G_EXTRACT_ROW "GSparseMatrix::extract_row(int)"
  43 +#define G_EXTRACT_COL "GSparseMatrix::extract_col(int)"
  44 +#define G_INSERT_COL "GSparseMatrix::insert_col(GVector,int)"
  45 +#define G_INSERT_COL2 "GSparseMatrix::insert_col(double*,int*,int,int)"
  46 +#define G_STACK_INIT "GSparseMatrix::stack_init(int,int)"
  47 +#define G_STACK_PUSH "GSparseMatrix::stack_push_column(double*,int*,int,int)"
  48 +#define G_STACK_FLUSH "GSparseMatrix::stack_flush(void)"
  49 +#define G_CONSTRUCTOR "GSparseMatrix::constructor(int,int,int)"
  50 +#define G_COPY_MEMBERS "GSparseMatrix::copy_members(GSparseMatrix)"
  51 +#define G_ALLOC "GSparseMatrix::alloc_elements(int,int)"
  52 +#define G_FREE "GSparseMatrix::free_elements(int,int)"
  53 +#define G_REMOVE_ZERO "GSparseMatrix::remove_zero_row_col(void)"
  54 +#define G_SYMPERM "cs_symperm(GSparseMatrix*,int*,int)"
  55 +#define G_TRANSPOSE "cs_transpose(GSparseMatrix*,int)"
57 56
58 /* __ Macros _____________________________________________________________ */ 57 /* __ Macros _____________________________________________________________ */
59 #define G_MIN(a,b) (((a) < (b)) ? (a) : (b)) 58 #define G_MIN(a,b) (((a) < (b)) ? (a) : (b))
@@ -80,126 +79,62 @@ using namespace std; @@ -80,126 +79,62 @@ using namespace std;
80 = = 79 = =
81 ==========================================================================*/ 80 ==========================================================================*/
82 81
83 -/***************************************************************************  
84 - * Constructor *  
85 - * ----------------------------------------------------------------------- *  
86 - * First calls the void base class constructor that initialises all base *  
87 - * class members, then call the GSparseMatrix initialisation method and *  
88 - * finally construct the object. * 82 +/***********************************************************************//**
  83 + * @brief Matrix constructor
  84 + *
  85 + * @param[in] rows Number of rows in matrix.
  86 + * @param[in] cols Number of columns in matrix.
  87 + * @param[in] elements Number of elements for which memory should be allocated.
  88 + *
  89 + * First calls the void base class constructor that initialises all base
  90 + * class members, then call the GSparseMatrix initialisation method and
  91 + * finally construct the object.
89 ***************************************************************************/ 92 ***************************************************************************/
90 GSparseMatrix::GSparseMatrix(int rows, int cols, int elements) : GMatrixBase() 93 GSparseMatrix::GSparseMatrix(int rows, int cols, int elements) : GMatrixBase()
91 { 94 {
92 - // Initialise private members for clean destruction  
93 - init_members(); 95 + // Initialise private members for clean destruction
  96 + init_members();
94 97
95 - // Construct sparse matrix  
96 - constructor(rows, cols, elements); 98 + // Construct sparse matrix
  99 + constructor(rows, cols, elements);
97 100
98 - // Return  
99 - return; 101 + // Return
  102 + return;
100 } 103 }
101 104
102 105
103 -/***************************************************************************  
104 - * Copy constructor *  
105 - * ----------------------------------------------------------------------- *  
106 - * First calls the void base class copy constructor that copys all base *  
107 - * class members, then call the GSparseMatrix member copy method. *  
108 - * NOTE: The copy constructor does NOT fill a pending element. * 106 +/***********************************************************************//**
  107 + * @brief Copy constructor
  108 + *
  109 + * @param[in] m Matrix from which class should be instantiated.
  110 + *
  111 + * First calls the void base class copy constructor that copys all base
  112 + * class members, then call the GSparseMatrix member copy method.
  113 + * NOTE: The copy constructor does NOT fill a pending element.
109 ***************************************************************************/ 114 ***************************************************************************/
110 GSparseMatrix::GSparseMatrix(const GSparseMatrix& m) : GMatrixBase(m) 115 GSparseMatrix::GSparseMatrix(const GSparseMatrix& m) : GMatrixBase(m)
111 { 116 {
112 - // Initialise private members for clean destruction  
113 - init_members(); 117 + // Initialise private members for clean destruction
  118 + init_members();
114 119
115 - // Copy members  
116 - copy_members(m); 120 + // Copy members
  121 + copy_members(m);
117 122
118 - // Return  
119 - return; 123 + // Return
  124 + return;
120 } 125 }
121 126
122 127
123 -/***************************************************************************  
124 - * GMatrix -> GSparseMatrix storage class conversion *  
125 - * ----------------------------------------------------------------------- *  
126 - * First calls the void base class constructor that initialises all base *  
127 - * class members, then call the GSparseMatrix initialisation method and *  
128 - * finally construct the object and fill it with the GMatrix data. * 128 +/***********************************************************************//**
  129 + * @brief Destructor
129 ***************************************************************************/ 130 ***************************************************************************/
130 -/*  
131 -GSparseMatrix::GSparseMatrix(const GMatrix& m) : GMatrixBase()  
132 -{  
133 - // Initialise private members for clean destruction  
134 - init_members();  
135 -  
136 - // Construct sparse matrix  
137 - constructor(m.rows(), m.cols());  
138 -  
139 - // Use a fill-stack for better performance  
140 - stack_init(m_cols, G_SPARSE_MATRIX_DEFAULT_STACK_ENTRIES);  
141 -  
142 - // Loop over all columns  
143 - GVector column(m_cols);  
144 - for (int col = 0; col < m_cols; ++col) {  
145 - column = m.extract_col(col);  
146 - insert_col(column, col);  
147 - }  
148 -  
149 - // Destroy fill stack  
150 - stack_destroy();  
151 -  
152 - // Return result  
153 - return;  
154 -}  
155 -*/  
156 -  
157 -/***************************************************************************  
158 - * GSymMatrix -> GSparseMatrix storage class conversion *  
159 - * ----------------------------------------------------------------------- *  
160 - * First calls the void base class constructor that initialises all base *  
161 - * class members, then call the GSparseMatrix initialisation method and *  
162 - * finally construct the object and fill it with the GSymMatrix data. *  
163 - ***************************************************************************/  
164 -/*  
165 -GSparseMatrix::GSparseMatrix(const GSymMatrix& m) : GMatrixBase()  
166 -{  
167 - // Initialise private members for clean destruction  
168 - init_members();  
169 -  
170 - // Construct sparse matrix  
171 - constructor(m.rows(), m.cols());  
172 -  
173 - // Use a fill-stack for better performance  
174 - stack_init(m_cols, G_SPARSE_MATRIX_DEFAULT_STACK_ENTRIES);  
175 -  
176 - // Loop over all columns  
177 - GVector column(m_cols);  
178 - for (int col = 0; col < m_cols; ++col) {  
179 - column = m.extract_col(col);  
180 - insert_col(column, col);  
181 - }  
182 -  
183 - // Destroy fill stack  
184 - stack_destroy();  
185 -  
186 - // Return result  
187 - return;  
188 -}  
189 -*/  
190 -  
191 -/***************************************************************************  
192 - * Destructor *  
193 - * ----------------------------------------------------------------------- *  
194 - * First destroys class members, then destroy base class members. *  
195 - ***************************************************************************/  
196 -GSparseMatrix::~GSparseMatrix() 131 +GSparseMatrix::~GSparseMatrix(void)
197 { 132 {
198 - // Free members  
199 - free_members(); 133 + // Free members
  134 + free_members();
200 135
201 - // Return  
202 - return; 136 + // Return
  137 + return;
203 } 138 }
204 139
205 140
@@ -209,237 +144,255 @@ GSparseMatrix::~GSparseMatrix() @@ -209,237 +144,255 @@ GSparseMatrix::~GSparseMatrix()
209 = = 144 = =
210 ==========================================================================*/ 145 ==========================================================================*/
211 146
212 -/***************************************************************************  
213 - * Assignment operator *  
214 - * ----------------------------------------------------------------------- * 147 +/***********************************************************************//**
  148 + * @brief Assignment operator
  149 + *
  150 + * @param[in] m GSparseMatrix instance to be assigned
215 ***************************************************************************/ 151 ***************************************************************************/
216 GSparseMatrix& GSparseMatrix::operator= (const GSparseMatrix& m) 152 GSparseMatrix& GSparseMatrix::operator= (const GSparseMatrix& m)
217 { 153 {
218 - // Execute only if object is not identical  
219 - if (this != &m) { 154 + // Execute only if object is not identical
  155 + if (this != &m) {
220 156
221 - // Copy base class members  
222 - this->GMatrixBase::operator=(m); 157 + // Copy base class members
  158 + this->GMatrixBase::operator=(m);
223 159
224 - // Free members  
225 - free_members(); 160 + // Free members
  161 + free_members();
226 162
227 - // Initialise private members for clean destruction  
228 - init_members(); 163 + // Initialise private members for clean destruction
  164 + init_members();
229 165
230 - // Copy members  
231 - copy_members(m); 166 + // Copy members
  167 + copy_members(m);
232 168
233 - } // endif: object was not identical 169 + } // endif: object was not identical
234 170
235 - // Return  
236 - return *this; 171 + // Return
  172 + return *this;
237 } 173 }
238 174
239 175
240 -/***************************************************************************  
241 - * Vector multiplication *  
242 - * ----------------------------------------------------------------------- *  
243 - * Performs sparse matrix*vector multiplication including an eventually *  
244 - * pending fill. * 176 +/***********************************************************************//**
  177 + * @brief Vector multiplication
  178 + *
  179 + * @param[in] v GVector with which matrix is to be multiplied
  180 + *
  181 + * Performs sparse matrix*vector multiplication including any pending fill.
245 ***************************************************************************/ 182 ***************************************************************************/
246 GVector GSparseMatrix::operator* (const GVector& v) const 183 GVector GSparseMatrix::operator* (const GVector& v) const
247 { 184 {
248 - // Raise an exception if the matrix and vector dimensions are not compatible  
249 - if (m_cols != v.m_num)  
250 - throw GException::matrix_vector_mismatch(G_OP_MUL_VEC, v.m_num, m_rows, m_cols);  
251 -  
252 - // Initialise result vector  
253 - GVector result(m_rows);  
254 -  
255 - // Multiply only if there are elements in matrix  
256 - if (m_elements > 0) {  
257 -  
258 - // Perform vector multiplication  
259 - for (int col = 0; col < m_cols; ++col) {  
260 - int i_start = m_colstart[col];  
261 - int i_stop = m_colstart[col+1];  
262 - double* ptr_data = m_data + i_start;  
263 - int* ptr_rowinx = m_rowinx + i_start;  
264 - for (int i = i_start; i < i_stop; ++i)  
265 - result(*ptr_rowinx++) += *ptr_data++ * v(col);  
266 - } 185 + // Raise an exception if the matrix and vector dimensions are not compatible
  186 + if (m_cols != v.m_num)
  187 + throw GException::matrix_vector_mismatch(G_OP_MUL_VEC, v.m_num,
  188 + m_rows, m_cols);
  189 +
  190 + // Initialise result vector
  191 + GVector result(m_rows);
  192 +
  193 + // Multiply only if there are elements in matrix
  194 + if (m_elements > 0) {
  195 +
  196 + // Perform vector multiplication
  197 + for (int col = 0; col < m_cols; ++col) {
  198 + int i_start = m_colstart[col];
  199 + int i_stop = m_colstart[col+1];
  200 + double* ptr_data = m_data + i_start;
  201 + int* ptr_rowinx = m_rowinx + i_start;
  202 + for (int i = i_start; i < i_stop; ++i)
  203 + result(*ptr_rowinx++) += *ptr_data++ * v(col);
  204 + }
267 205
268 - // If fill is pending then add-in also this element into the product  
269 - if (m_fill_val != 0.0) {  
270 - result(m_fill_row) += m_fill_val * v(m_fill_col);  
271 - #if defined(G_DEBUG_SPARSE_PENDING)  
272 - cout << G_OP_MUL_VEC << ": pending value " << m_fill_val <<  
273 - " for location (" << m_fill_row << "," << m_fill_col <<  
274 - ") has been used." << endl;  
275 - #endif  
276 - } 206 + // If fill is pending then add-in also this element into the product
  207 + if (m_fill_val != 0.0) {
  208 + result(m_fill_row) += m_fill_val * v(m_fill_col);
  209 + #if defined(G_DEBUG_SPARSE_PENDING)
  210 + std::cout << G_OP_MUL_VEC << ": pending value " << m_fill_val
  211 + << " for location (" << m_fill_row << "," << m_fill_col
  212 + << ") has been used." << std::endl;
  213 + #endif
  214 + }
277 215
278 - } // endif: there were elements in matrix 216 + } // endif: there were elements in matrix
279 217
280 - // Return result  
281 - return result; 218 + // Return result
  219 + return result;
282 } 220 }
283 221
284 222
285 -/***************************************************************************  
286 - * Equalty operator *  
287 - * ----------------------------------------------------------------------- *  
288 - * Two matrixes are considered equal if they have the same dimensions and *  
289 - * identicial elements. *  
290 - * !!! NOT YET IN NATIVE CODING !!! * 223 +/***********************************************************************//**
  224 + * @brief Equalty operator
  225 + *
  226 + * @param[in] m Matrix to compare with.
  227 + *
  228 + * Two matrixes are considered equal if they have the same dimensions and
  229 + * identicial elements. The return value is 1 if the matrices are identical,
  230 + * 0 otherwise.
  231 + * TODO: Implement native sparse code
291 ***************************************************************************/ 232 ***************************************************************************/
292 int GSparseMatrix::operator== (const GSparseMatrix &m) const 233 int GSparseMatrix::operator== (const GSparseMatrix &m) const
293 { 234 {
294 - // Initalise the result to 'equal matrices'  
295 - int result = 1;  
296 -  
297 - // Perform comparison (only if matrix dimensions are identical)  
298 - if (m_rows == m.m_rows && m_cols == m.m_cols) {  
299 -  
300 - // Loop over all matrix elements  
301 - for (int row = 0; row < m_rows; ++row) {  
302 - for (int col = 0; col < m_cols; ++col) {  
303 - if ((*this)(row,col) != m(row,col)) {  
304 - result = 0;  
305 - break; 235 + // Initalise the result to 'equal matrices'
  236 + int result = 1;
  237 +
  238 + // Perform comparison (only if matrix dimensions are identical)
  239 + if (m_rows == m.m_rows && m_cols == m.m_cols) {
  240 +
  241 + // Loop over all matrix elements
  242 + for (int row = 0; row < m_rows; ++row) {
  243 + for (int col = 0; col < m_cols; ++col) {
  244 + if ((*this)(row,col) != m(row,col)) {
  245 + result = 0;
  246 + break;
  247 + }
  248 + }
  249 + if (!result)
  250 + break;
306 } 251 }
307 - }  
308 - if (!result)  
309 - break;  
310 } 252 }
311 - }  
312 - else  
313 - result = 0; 253 + else
  254 + result = 0;
314 255
315 - // Return result  
316 - return result; 256 + // Return result
  257 + return result;
317 } 258 }
318 259
319 260
320 -/***************************************************************************  
321 - * Equalty operator *  
322 - * ----------------------------------------------------------------------- *  
323 - * Two matrixes are considered equal if they have the same dimensions and *  
324 - * identicial elements. *  
325 - * !!! NOT YET IN NATIVE CODING !!! * 261 +/***********************************************************************//**
  262 + * @brief Unequalty operator
  263 + *
  264 + * @param[in] m Matrix to compare with.
  265 + *
  266 + * Two matrixes are considered unequal if they have different dimensions or
  267 + * differing elements. The return value is 1 if the matrices are not
  268 + * identical and 0 if they are identical.
  269 + * TODO: Implement native sparse code
326 ***************************************************************************/ 270 ***************************************************************************/
327 int GSparseMatrix::operator!= (const GSparseMatrix &m) const 271 int GSparseMatrix::operator!= (const GSparseMatrix &m) const
328 { 272 {
329 - // Initalise the result to 'equal matrices'  
330 - int result = 0;  
331 -  
332 - // Perform comparison (only if matrix dimensions are identical)  
333 - if (m_rows == m.m_rows && m_cols == m.m_cols) {  
334 -  
335 - // Loop over all matrix elements  
336 - for (int row = 0; row < m_rows; ++row) {  
337 - for (int col = 0; col < m_cols; ++col) {  
338 - if ((*this)(row,col) != m(row,col)) {  
339 - result = 1;  
340 - break; 273 + // Initalise the result to 'equal matrices'
  274 + int result = 0;
  275 +
  276 + // Perform comparison (only if matrix dimensions are identical)
  277 + if (m_rows == m.m_rows && m_cols == m.m_cols) {
  278 +
  279 + // Loop over all matrix elements
  280 + for (int row = 0; row < m_rows; ++row) {
  281 + for (int col = 0; col < m_cols; ++col) {
  282 + if ((*this)(row,col) != m(row,col)) {
  283 + result = 1;
  284 + break;
  285 + }
  286 + }
  287 + if (!result)
  288 + break;
341 } 289 }
342 - }  
343 - if (!result)  
344 - break;  
345 } 290 }
346 - }  
347 - else  
348 - result = 1; 291 + else
  292 + result = 1;
349 293
350 - // Return result  
351 - return result; 294 + // Return result
  295 + return result;
352 } 296 }
353 297
354 298
355 -/***************************************************************************  
356 - * Unary matrix addition operator *  
357 - * ----------------------------------------------------------------------- *  
358 - * !!! NOT YET IN NATIVE CODING !!! * 299 +/***********************************************************************//**
  300 + * @brief Unary matrix addition operator
  301 + *
  302 + * @param[in] m Matrix to be added.
  303 + *
  304 + * TODO: Implement native sparse code
359 ***************************************************************************/ 305 ***************************************************************************/
360 GSparseMatrix& GSparseMatrix::operator+= (const GSparseMatrix& m) 306 GSparseMatrix& GSparseMatrix::operator+= (const GSparseMatrix& m)
361 { 307 {
362 - // Raise an exception if the matrix dimensions are not compatible  
363 - if (m_rows != m.m_rows || m_cols != m.m_cols)  
364 - throw GException::matrix_mismatch(G_OP_ADD, m_rows, m_cols, m.m_rows, m.m_cols); 308 + // Raise an exception if the matrix dimensions are not compatible
  309 + if (m_rows != m.m_rows || m_cols != m.m_cols)
  310 + throw GException::matrix_mismatch(G_OP_ADD, m_rows, m_cols,
  311 + m.m_rows, m.m_cols);
365 312
366 - // Perform inplace matrix addition  
367 - for (int col = 0; col < m_cols; ++col) {  
368 - GVector v_result = extract_col(col);  
369 - GVector v_operand = m.extract_col(col);  
370 - v_result += v_operand;  
371 - insert_col(v_result, col);  
372 - } 313 + // Perform inplace matrix addition using vectors
  314 + for (int col = 0; col < m_cols; ++col) {
  315 + GVector v_result = extract_col(col);
  316 + GVector v_operand = m.extract_col(col);
  317 + v_result += v_operand;
  318 + insert_col(v_result, col);
  319 + }
373 320
374 - // Return result  
375 - return *this; 321 + // Return result
  322 + return *this;
376 } 323 }
377 324
378 325
379 -/***************************************************************************  
380 - * Unary matrix subtraction operator *  
381 - * ----------------------------------------------------------------------- *  
382 - * !!! NOT YET IN NATIVE CODING !!! * 326 +/***********************************************************************//**
  327 + * @brief Unary matrix subtraction operator
  328 + *
  329 + * @param[in] m Matrix to be subtracted.
  330 + *
  331 + * TODO: Implement native sparse code
383 ***************************************************************************/ 332 ***************************************************************************/
384 GSparseMatrix& GSparseMatrix::operator-= (const GSparseMatrix& m) 333 GSparseMatrix& GSparseMatrix::operator-= (const GSparseMatrix& m)
385 { 334 {
386 - // Raise an exception if the matrix dimensions are not compatible  
387 - if (m_rows != m.m_rows || m_cols != m.m_cols)  
388 - throw GException::matrix_mismatch(G_OP_SUB, m_rows, m_cols, m.m_rows, m.m_cols); 335 + // Raise an exception if the matrix dimensions are not compatible
  336 + if (m_rows != m.m_rows || m_cols != m.m_cols)
  337 + throw GException::matrix_mismatch(G_OP_SUB, m_rows, m_cols,
  338 + m.m_rows, m.m_cols);
389 339
390 - // Perform inplace matrix subtraction  
391 - for (int col = 0; col < m_cols; ++col) {  
392 - GVector v_result = extract_col(col);  
393 - GVector v_operand = m.extract_col(col);  
394 - v_result -= v_operand;  
395 - insert_col(v_result, col);  
396 - } 340 + // Perform inplace matrix subtraction
  341 + for (int col = 0; col < m_cols; ++col) {
  342 + GVector v_result = extract_col(col);
  343 + GVector v_operand = m.extract_col(col);
  344 + v_result -= v_operand;
  345 + insert_col(v_result, col);
  346 + }
397 347
398 - // Return result  
399 - return *this; 348 + // Return result
  349 + return *this;
400 } 350 }
401 351
402 352
403 -/***************************************************************************  
404 - * Matrix multiplication (unary multiplication operator) *  
405 - * ----------------------------------------------------------------------- *  
406 - * !!! NOT YET IN NATIVE CODING !!! * 353 +/***********************************************************************//**
  354 + * @brief Unary matrix multiplication operator
  355 + *
  356 + * @param[in] m Matrix to be multiplied.
  357 + *
  358 + * TODO: Implement native sparse code
407 ***************************************************************************/ 359 ***************************************************************************/
408 GSparseMatrix& GSparseMatrix::operator*= (const GSparseMatrix& m) 360 GSparseMatrix& GSparseMatrix::operator*= (const GSparseMatrix& m)
409 { 361 {
410 - // Raise an exception if the matrix dimensions are not compatible  
411 - if (m_cols != m.m_rows || m_rows != m.m_cols)  
412 - throw GException::matrix_mismatch(G_OP_MAT_MUL, m_rows, m_cols, m.m_rows, m.m_cols);  
413 -  
414 - // Allocate result matrix  
415 - GSparseMatrix result(m_rows, m.m_cols);  
416 -  
417 - // Multiply only if there are elements in both matrices  
418 - if (m_elements > 0 && m.m_elements > 0) { 362 + // Raise an exception if the matrix dimensions are not compatible
  363 + if (m_cols != m.m_rows || m_rows != m.m_cols)
  364 + throw GException::matrix_mismatch(G_OP_MAT_MUL, m_rows, m_cols,
  365 + m.m_rows, m.m_cols);
  366 +
  367 + // Allocate result matrix
  368 + GSparseMatrix result(m_rows, m.m_cols);
  369 +
  370 + // Multiply only if there are elements in both matrices
  371 + if (m_elements > 0 && m.m_elements > 0) {
  372 +
  373 + // Loop over all elements of result matrix
  374 + for (int row = 0; row < m_rows; ++row) {
  375 + for (int col = 0; col < m.m_cols; ++col) {
  376 + double sum = 0.0;
  377 + for (int i = 0; i < m_cols; ++i)
  378 + sum += (*this)(row,i) * m(i,col);
  379 + result(row,col) = sum;
  380 + }
  381 + }
419 382
420 - // Loop over all elements of result matrix  
421 - for (int row = 0; row < m_rows; ++row) {  
422 - for (int col = 0; col < m.m_cols; ++col) {  
423 - double sum = 0.0;  
424 - for (int i = 0; i < m_cols; ++i)  
425 - sum += (*this)(row,i) * m(i,col);  
426 - result(row,col) = sum;  
427 - }  
428 } 383 }
429 384
430 - }  
431 -  
432 - // Assign result  
433 - *this = result; 385 + // Assign result
  386 + *this = result;
434 387
435 - // Return result  
436 - return *this; 388 + // Return result
  389 + return *this;
437 } 390 }
438 391
439 392
440 /*========================================================================== 393 /*==========================================================================
441 = = 394 = =
442 - = GSparseMatrix methods = 395 + = GSparseMatrix public methods =
443 = = 396 = =
444 ==========================================================================*/ 397 ==========================================================================*/
445 398
@@ -714,214 +667,215 @@ void GSparseMatrix::cholesky_decompose(int compress) @@ -714,214 +667,215 @@ void GSparseMatrix::cholesky_decompose(int compress)
714 } 667 }
715 668
716 669
717 -/***************************************************************************  
718 - * GSparseMatrix Cholesky solver *  
719 - * ----------------------------------------------------------------------- *  
720 - * Solves the linear equation A*x=b using a Cholesky decomposition of A. *  
721 - * This function is to be applied on a GSparseMatrix matrix for which a *  
722 - * Choleksy factorization has been produced using 'cholesky_decompose'. * 670 +/***********************************************************************//**
  671 + * @brief Cholesky solver
  672 + *
  673 + * @param[in] v Solution vector.
  674 + * @param[in] compress Request matrix compression.
  675 + *
  676 + * Solves the linear equation A*x=b using a Cholesky decomposition of A.
  677 + * This function is to be applied on a GSparseMatrix matrix for which a
  678 + * Choleksy factorization has been produced using 'cholesky_decompose'.
723 ***************************************************************************/ 679 ***************************************************************************/
724 GVector GSparseMatrix::cholesky_solver(const GVector& v, int compress) 680 GVector GSparseMatrix::cholesky_solver(const GVector& v, int compress)
725 { 681 {
726 - // Declare loop variables  
727 - int row;  
728 - int col;  
729 - int c_row;  
730 - int c_col;  
731 -  
732 - // Dump header  
733 - #if defined(G_DEBUG_SPARSE_COMPRESSION)  
734 - cout << "GSparseMatrix::cholesky_solver" << endl;  
735 - cout << " Input vector .....: " << v << endl;  
736 - #endif  
737 -  
738 - // Raise an exception if the matrix and vector dimensions are incompatible  
739 - if (m_rows != v.m_num)  
740 - throw GException::matrix_vector_mismatch(G_CHOL_SOLVE, v.m_num, m_rows, m_cols);  
741 -  
742 - // Raise an exception if there is no symbolic pointer  
743 - if (!m_symbolic)  
744 - throw GException::matrix_not_factorised(G_CHOL_SOLVE, "Cholesky decomposition");  
745 -  
746 - // Get symbolic pointer  
747 - GSparseSymbolic* symbolic = (GSparseSymbolic*)m_symbolic;  
748 -  
749 - // Raise an exception if there is no permutation  
750 - if (!symbolic->m_pinv)  
751 - throw GException::matrix_not_factorised(G_CHOL_SOLVE, "Cholesky decomposition");  
752 -  
753 - // Flag row and column compression  
754 - int row_compressed = (m_rowsel != NULL && m_num_rowsel < m_rows);  
755 - int col_compressed = (m_colsel != NULL && m_num_colsel < m_cols);  
756 -  
757 - // Decide if we need a compression algorithm or not  
758 - int no_zero = !(compress && (row_compressed || col_compressed));  
759 -  
760 - // Allocate vector for permutation and result vector  
761 - GVector result(m_cols);  
762 -  
763 - // Setup pointers to L matrix and x vector  
764 - int* Lp = m_colstart;  
765 - int* Li = m_rowinx;  
766 - double* Lx = m_data;  
767 -  
768 - // Case A: no zero-row/col compression needed  
769 - if (no_zero) {  
770 -  
771 - // Declare working vector  
772 - GVector x(m_rows);  
773 -  
774 - // Perform inverse vector permutation  
775 - for (int i = 0; i < v.m_num; ++i)  
776 - x.m_data[symbolic->m_pinv[i]] = v.m_data[i]; 682 + // Declare loop variables
  683 + int row;
  684 + int col;
  685 + int c_row;
  686 + int c_col;
777 687
778 - // Inplace solve L\x=x  
779 - for (col = 0; col < m_cols; ++col) { // loop over columns  
780 - x.m_data[col] /= Lx[Lp[col]]; // divide by diag.  
781 - for (int p = Lp[col]+1; p < Lp[col+1]; p++) // loop over elements  
782 - x.m_data[Li[p]] -= Lx[p] * x.m_data[col];  
783 - } 688 + // Dump header
  689 + #if defined(G_DEBUG_SPARSE_COMPRESSION)
  690 + std::cout << "GSparseMatrix::cholesky_solver" << std::endl;
  691 + std::cout << " Input vector .....: " << v << std::endl;
  692 + #endif
784 693
785 - // Inplace solve L'\x=x  
786 - for (col = m_cols-1; col >= 0; --col) { // loop over columns  
787 - for (int p = Lp[col]+1; p < Lp[col+1]; p++) // loop over elements  
788 - x.m_data[col] -= Lx[p] * x.m_data[Li[p]];  
789 - x.m_data[col] /= Lx[Lp[col]];  
790 - } 694 + // Raise an exception if the matrix and vector dimensions are incompatible
  695 + if (m_rows != v.m_num)
  696 + throw GException::matrix_vector_mismatch(G_CHOL_SOLVE, v.m_num,
  697 + m_rows, m_cols);
  698 + // Raise an exception if there is no symbolic pointer
  699 + if (!m_symbolic)
  700 + throw GException::matrix_not_factorised(G_CHOL_SOLVE,
  701 + "Cholesky decomposition");
  702 +
  703 + // Get symbolic pointer
  704 + GSparseSymbolic* symbolic = (GSparseSymbolic*)m_symbolic;
  705 +
  706 + // Raise an exception if there is no permutation
  707 + if (!symbolic->m_pinv)
  708 + throw GException::matrix_not_factorised(G_CHOL_SOLVE,
  709 + "Cholesky decomposition");
  710 +
  711 + // Flag row and column compression
  712 + int row_compressed = (m_rowsel != NULL && m_num_rowsel < m_rows);
  713 + int col_compressed = (m_colsel != NULL && m_num_colsel < m_cols);
  714 +
  715 + // Decide if we need a compression algorithm or not
  716 + int no_zero = !(compress && (row_compressed || col_compressed));
  717 +
  718 + // Allocate vector for permutation and result vector
  719 + GVector result(m_cols);
  720 +
  721 + // Setup pointers to L matrix and x vector
  722 + int* Lp = m_colstart;
  723 + int* Li = m_rowinx;
  724 + double* Lx = m_data;
  725 +
  726 + // Case A: no zero-row/col compression needed
  727 + if (no_zero) {
  728 +
  729 + // Declare working vector
  730 + GVector x(m_rows);
  731 +
  732 + // Perform inverse vector permutation
  733 + for (int i = 0; i < v.m_num; ++i)
  734 + x.m_data[symbolic->m_pinv[i]] = v.m_data[i];
  735 +
  736 + // Inplace solve L\x=x
  737 + for (col = 0; col < m_cols; ++col) { // loop over columns
  738 + x.m_data[col] /= Lx[Lp[col]]; // divide by diag.
  739 + for (int p = Lp[col]+1; p < Lp[col+1]; p++) // loop over elements
  740 + x.m_data[Li[p]] -= Lx[p] * x.m_data[col];
  741 + }
791 742
792 - // Perform vector permutation  
793 - for (int i = 0; i < m_cols; ++i)  
794 - result.m_data[i] = x.m_data[symbolic->m_pinv[i]]; 743 + // Inplace solve L'\x=x
  744 + for (col = m_cols-1; col >= 0; --col) { // loop over columns
  745 + for (int p = Lp[col]+1; p < Lp[col+1]; p++) // loop over elements
  746 + x.m_data[col] -= Lx[p] * x.m_data[Li[p]];
  747 + x.m_data[col] /= Lx[Lp[col]];
  748 + }
795 749
796 - } // endif: Case A 750 + // Perform vector permutation
  751 + for (int i = 0; i < m_cols; ++i)
  752 + result.m_data[i] = x.m_data[symbolic->m_pinv[i]];
797 753
798 - // Case B: zero-row/column compression requested  
799 - else { 754 + } // endif: Case A
800 755
801 - // Allocate row and column mapping arrays  
802 - int* row_map = new int[m_rows];  
803 - if (row_map == NULL)  
804 - throw GException::mem_alloc(G_CHOL_SOLVE, m_rows);  
805 - int* col_map = new int[m_cols];  
806 - if (col_map == NULL)  
807 - throw GException::mem_alloc(G_CHOL_SOLVE, m_cols);  
808 -  
809 - // Setup row mapping array that maps original matrix rows into compressed  
810 - // matrix rows. An entry of -1 indicates that the row should be dropped.  
811 - // If no selection exists then setup an identity map.  
812 - if (row_compressed) {  
813 - for (row = 0; row < m_rows; ++row)  
814 - row_map[row] = -1;  
815 - for (c_row = 0; c_row < m_num_rowsel; ++c_row)  
816 - row_map[m_rowsel[c_row]] = c_row;  
817 - } 756 + // Case B: zero-row/column compression requested
818 else { 757 else {
819 - for (row = 0; row < m_rows; ++row)  
820 - row_map[row] = row;  
821 - }  
822 - #if defined(G_DEBUG_SPARSE_COMPRESSION)  
823 - cout << " Row mapping ......:";  
824 - for (row = 0; row < m_rows; ++row)  
825 - cout << " " << row_map[row];  
826 - cout << endl;  
827 - #endif  
828 758
829 - // Setup column mapping array that maps original matrix column into compressed  
830 - // matrix columns. An entry of -1 indicates that the column should be dropped  
831 - // If no selection exists then setup an identity map.  
832 - if (col_compressed) {  
833 - for (col = 0; col < m_cols; ++col)  
834 - col_map[col] = -1;  
835 - for (c_col = 0; c_col < m_num_colsel; ++c_col)  
836 - col_map[m_colsel[c_col]] = c_col;  
837 - }  
838 - else {  
839 - for (col = 0; col < m_cols; ++col)  
840 - col_map[col] = col;  
841 - }  
842 - #if defined(G_DEBUG_SPARSE_COMPRESSION)  
843 - cout << " Column mapping ...:";  
844 - for (col = 0; col < m_cols; ++col)  
845 - cout << " " << col_map[col];  
846 - cout << endl;  
847 - #endif 759 + // Allocate row and column mapping arrays
  760 + int* row_map = new int[m_rows];
  761 + int* col_map = new int[m_cols];
  762 +
  763 + // Setup row mapping array that maps original matrix rows into compressed
  764 + // matrix rows. An entry of -1 indicates that the row should be dropped.
  765 + // If no selection exists then setup an identity map.
  766 + if (row_compressed) {
  767 + for (row = 0; row < m_rows; ++row)
  768 + row_map[row] = -1;
  769 + for (c_row = 0; c_row < m_num_rowsel; ++c_row)
  770 + row_map[m_rowsel[c_row]] = c_row;
  771 + }
  772 + else {
  773 + for (row = 0; row < m_rows; ++row)
  774 + row_map[row] = row;
  775 + }
  776 + #if defined(G_DEBUG_SPARSE_COMPRESSION)
  777 + std::cout << " Row mapping ......:";
  778 + for (row = 0; row < m_rows; ++row)
  779 + std::cout << " " << row_map[row];
  780 + std::cout << std::endl;
  781 + #endif
848 782
849 - // Declare working vector  
850 - GVector x(row_compressed ? m_num_rowsel : m_rows); 783 + // Setup column mapping array that maps original matrix column into compressed
  784 + // matrix columns. An entry of -1 indicates that the column should be dropped
  785 + // If no selection exists then setup an identity map.
  786 + if (col_compressed) {
  787 + for (col = 0; col < m_cols; ++col)
  788 + col_map[col] = -1;
  789 + for (c_col = 0; c_col < m_num_colsel; ++c_col)
  790 + col_map[m_colsel[c_col]] = c_col;
  791 + }
  792 + else {
  793 + for (col = 0; col < m_cols; ++col)
  794 + col_map[col] = col;
  795 + }
  796 + #if defined(G_DEBUG_SPARSE_COMPRESSION)
  797 + std::cout << " Column mapping ...:";
  798 + for (col = 0; col < m_cols; ++col)
  799 + std::cout << " " << col_map[col];
  800 + std::cout << std::endl;
  801 + #endif
851 802
852 - // Compress input vector v -> c_v if required  
853 - if (m_rowsel != NULL && m_num_rowsel < m_rows) {  
854 - for (c_row = 0; c_row < m_num_rowsel; ++c_row)  
855 - x.m_data[c_row] = v.m_data[m_rowsel[c_row]];  
856 - }  
857 - else  
858 - x = v;  
859 - #if defined(G_DEBUG_SPARSE_COMPRESSION)  
860 - cout << " Compressed vector : " << x << endl;  
861 - #endif 803 + // Declare working vector
  804 + GVector x(row_compressed ? m_num_rowsel : m_rows);
862 805
863 - // Perform inverse permutation  
864 - x = iperm(x, symbolic->m_pinv);  
865 - #if defined(G_DEBUG_SPARSE_COMPRESSION)  
866 - cout << " Permutated vector : " << x << endl;  
867 - #endif 806 + // Compress input vector v -> c_v if required
  807 + if (m_rowsel != NULL && m_num_rowsel < m_rows) {
  808 + for (c_row = 0; c_row < m_num_rowsel; ++c_row)
  809 + x.m_data[c_row] = v.m_data[m_rowsel[c_row]];
  810 + }
  811 + else
  812 + x = v;
  813 + #if defined(G_DEBUG_SPARSE_COMPRESSION)
  814 + std::cout << " Compressed vector : " << x << std::endl;
  815 + #endif
868 816
869 - // Inplace solve L\x=x. The column and row maps are just use to see which  
870 - // columns or rows should be skipped in the calculations.  
871 - for (col = 0; col < m_cols; ++col) { // loop over columns  
872 - if ((c_col = col_map[col]) >= 0) { // use only non-zero cols  
873 - x.m_data[c_col] /= Lx[Lp[col]]; // divide by diag.  
874 - for (int p = Lp[col]+1; p < Lp[col+1]; p++) { // loop over elements  
875 - if ((c_row = row_map[Li[p]]) >= 0) // use only non-zero rows  
876 - x.m_data[c_row] -= Lx[p] * x.m_data[c_col]; 817 + // Perform inverse permutation
  818 + x = iperm(x, symbolic->m_pinv);
  819 + #if defined(G_DEBUG_SPARSE_COMPRESSION)
  820 + std::cout << " Permutated vector : " << x << std::endl;
  821 + #endif
  822 +
  823 + // Inplace solve L\x=x. The column and row maps are just use to see which
  824 + // columns or rows should be skipped in the calculations.
  825 + for (col = 0; col < m_cols; ++col) { // loop over columns
  826 + if ((c_col = col_map[col]) >= 0) { // use only non-zero cols
  827 + x.m_data[c_col] /= Lx[Lp[col]]; // divide by diag.
  828 + for (int p = Lp[col]+1; p < Lp[col+1]; p++) { // loop over elements
  829 + if ((c_row = row_map[Li[p]]) >= 0) // use only non-zero rows
  830 + x.m_data[c_row] -= Lx[p] * x.m_data[c_col];
  831 + }
  832 + }
877 } 833 }
878 - }  
879 - }  
880 - #if defined(G_DEBUG_SPARSE_COMPRESSION)  
881 - cout << " Solve Lx=x .......: " << x << endl;  
882 - #endif 834 + #if defined(G_DEBUG_SPARSE_COMPRESSION)
  835 + std::cout << " Solve Lx=x .......: " << x << std::endl;
  836 + #endif
883 837
884 - // Inplace solve L'\x=x. The column and row maps are just use to see which  
885 - // columns or rows should be skipped in the calculations.  
886 - for (col = m_cols-1; col >= 0; --col) { // loop over columns  
887 - if ((c_col = col_map[col]) >= 0) { // use only non-zero cols  
888 - for (int p = Lp[col]+1; p < Lp[col+1]; p++) { // loop over elements  
889 - if ((c_row = row_map[Li[p]]) >= 0) // use only non-zero rows  
890 - x.m_data[c_col] -= Lx[p] * x.m_data[c_row]; 838 + // Inplace solve L'\x=x. The column and row maps are just use to see which
  839 + // columns or rows should be skipped in the calculations.
  840 + for (col = m_cols-1; col >= 0; --col) { // loop over columns
  841 + if ((c_col = col_map[col]) >= 0) { // use only non-zero cols
  842 + for (int p = Lp[col]+1; p < Lp[col+1]; p++) { // loop over elements
  843 + if ((c_row = row_map[Li[p]]) >= 0) // use only non-zero rows
  844 + x.m_data[c_col] -= Lx[p] * x.m_data[c_row];
  845 + }
  846 + x.m_data[c_col] /= Lx[Lp[col]];
  847 + }
891 } 848 }
892 - x.m_data[c_col] /= Lx[Lp[col]];  
893 - }  
894 - }  
895 - #if defined(G_DEBUG_SPARSE_COMPRESSION)  
896 - cout << " Solve L'x=x ......: " << x << endl;  
897 - #endif 849 + #if defined(G_DEBUG_SPARSE_COMPRESSION)
  850 + std::cout << " Solve L'x=x ......: " << x << std::endl;
  851 + #endif
898 852
899 - // Perform vector permutation  
900 - x = perm(x, symbolic->m_pinv);  
901 - #if defined(G_DEBUG_SPARSE_COMPRESSION)  
902 - cout << " Permutated vector : " << x << endl;  
903 - #endif 853 + // Perform vector permutation
  854 + x = perm(x, symbolic->m_pinv);
  855 + #if defined(G_DEBUG_SPARSE_COMPRESSION)
  856 + std::cout << " Permutated vector : " << x << std::endl;
  857 + #endif
904 858
905 - // If column compression has been performed the expand the result vector  
906 - // accordingly  
907 - if (m_colsel != NULL && m_num_colsel < m_cols) {  
908 - for (c_col = 0; c_col < m_num_colsel; ++c_col)  
909 - result.m_data[m_colsel[c_col]] = x.m_data[c_col];  
910 - }  
911 - else  
912 - result = x;  
913 - #if defined(G_DEBUG_SPARSE_COMPRESSION)  
914 - cout << " Restored vector ..: " << result << endl;  
915 - #endif 859 + // If column compression has been performed the expand the result vector
  860 + // accordingly
  861 + if (m_colsel != NULL && m_num_colsel < m_cols) {
  862 + for (c_col = 0; c_col < m_num_colsel; ++c_col)
  863 + result.m_data[m_colsel[c_col]] = x.m_data[c_col];
  864 + }
  865 + else
  866 + result = x;
  867 + #if defined(G_DEBUG_SPARSE_COMPRESSION)
  868 + std::cout << " Restored vector ..: " << result << std::endl;
  869 + #endif
916 870
917 - // Free mapping arrays  
918 - delete [] row_map;  
919 - delete [] col_map; 871 + // Free mapping arrays
  872 + delete [] row_map;
  873 + delete [] col_map;
920 874
921 - } // endelse: Case B 875 + } // endelse: Case B
922 876
923 - // Return result vector  
924 - return result; 877 + // Return result vector
  878 + return result;
925 } 879 }
926 880
927 881
@@ -963,21 +917,20 @@ void GSparseMatrix::cholesky_invert(int compress) @@ -963,21 +917,20 @@ void GSparseMatrix::cholesky_invert(int compress)
963 } 917 }
964 918
965 919
966 -/***************************************************************************  
967 - * Set all matrix elements to 0 *  
968 - * ----------------------------------------------------------------------- * 920 +/***********************************************************************//**
  921 + * @brief Set all matrix elements to 0
969 ***************************************************************************/ 922 ***************************************************************************/
970 -void GSparseMatrix::clear() 923 +void GSparseMatrix::clear(void)
971 { 924 {
972 - // Free memory  
973 - free_elements(0, m_elements); 925 + // Free memory
  926 + free_elements(0, m_elements);
974 927
975 - // Initialise column start indices to 0  
976 - for (int col = 0; col <= m_cols; ++col)  
977 - m_colstart[col] = 0; 928 + // Initialise column start indices to 0
  929 + for (int col = 0; col <= m_cols; ++col)
  930 + m_colstart[col] = 0;
978 931
979 - // Return  
980 - return; 932 + // Return
  933 + return;
981 } 934 }
982 935
983 936
@@ -1057,17 +1010,35 @@ GVector GSparseMatrix::extract_col(int col) const @@ -1057,17 +1010,35 @@ GVector GSparseMatrix::extract_col(int col) const
1057 } 1010 }
1058 1011
1059 1012
1060 -/***************************************************************************  
1061 - * Get matrix fill *  
1062 - * ----------------------------------------------------------------------- * 1013 +/***********************************************************************//**
  1014 + * @brief Returns fill of matrix
  1015 + *
  1016 + * The fill of a matrix is defined as the number non-zero elements devided
  1017 + * by the number of total elements. By definiton, the fill is comprised
  1018 + * in the interval [0,..,1]. The fill of an undefined matrix is defined to
  1019 + * be 0.
1063 ***************************************************************************/ 1020 ***************************************************************************/
1064 double GSparseMatrix::fill() const 1021 double GSparseMatrix::fill() const
1065 { 1022 {
1066 - // Determine number of elements in matrix  
1067 - int num = (m_fill_val == 0.0) ? m_elements : m_elements + 1; 1023 + // Initialise result
  1024 + double result = 0.0;
  1025 +
  1026 + // Compute matrix size
  1027 + int size = m_rows*m_cols;
1068 1028
1069 - // Return fill  
1070 - return (double(num) / double(m_rows*m_cols)); 1029 + // Continue only if matrix has elements
  1030 + if (size > 0) {
  1031 +
  1032 + // Determine number of elements in matrix
  1033 + int num = (m_fill_val == 0.0) ? m_elements : m_elements + 1;
  1034 +
  1035 + // Compute fill
  1036 + result = double(num) / double(size);
  1037 +
  1038 + } // endif: there were elements in matrix
  1039 +
  1040 + // Return fill
  1041 + return result;
1071 } 1042 }
1072 1043
1073 1044
@@ -1842,182 +1813,187 @@ double GSparseMatrix::sum() const @@ -1842,182 +1813,187 @@ double GSparseMatrix::sum() const
1842 1813
1843 /*========================================================================== 1814 /*==========================================================================
1844 = = 1815 = =
1845 - = GSparseMatrix private functions = 1816 + = GSparseMatrix private methods =
1846 = = 1817 = =
1847 ==========================================================================*/ 1818 ==========================================================================*/
1848 1819
1849 -/***************************************************************************  
1850 - * Constructor method *  
1851 - * ----------------------------------------------------------------------- *  
1852 - * This is the main constructor code, without any initialisation in it. *  
1853 - * It is used as service function to a number of different constructors of *  
1854 - * the GSparseMatrix class. * 1820 +/***********************************************************************//**
  1821 + * @brief Allocate matrix
  1822 + *
  1823 + * @param[in] rows Number of rows to be allocated.
  1824 + * @param[in] cols Number of columns to be allocated.
  1825 + * @param[in] elements Number of matrix elements to be physically allocated.
  1826 + *
  1827 + * This is the main constructor code that allocates and initialises memory
  1828 + * for matrix elements.
1855 ***************************************************************************/ 1829 ***************************************************************************/
1856 void GSparseMatrix::constructor(int rows, int cols, int elements) 1830 void GSparseMatrix::constructor(int rows, int cols, int elements)
1857 { 1831 {
1858 - // Allocate column start array. This is the only array that we can  
1859 - // allocate at this time. The other arrays can only be allocated during  
1860 - // filling of the matrix  
1861 - m_colstart = new int[cols+1];  
1862 - if (m_colstart == NULL)  
1863 - throw GException::mem_alloc(G_CONSTRUCTOR, cols+1);  
1864 -  
1865 - // Store (logical) matrix size  
1866 - m_rows = rows;  
1867 - m_cols = cols; 1832 + // Continue only if matrix is not empty
  1833 + if (rows > 0 && cols > 0) {
1868 1834
1869 - // Initialise column start indices to 0  
1870 - for (int col = 0; col <= m_cols; ++col)  
1871 - m_colstart[col] = 0; 1835 + // Allocate column start array. This is the only array that we can
  1836 + // allocate at this time. The other arrays can only be allocated during
  1837 + // filling of the matrix
  1838 + m_colstart = new int[cols+1];
1872 1839
1873 - // Optionally allocate elements  
1874 - if (elements > 0)  
1875 - alloc_elements(0, elements); 1840 + // Store (logical) matrix size
  1841 + m_rows = rows;
  1842 + m_cols = cols;
1876 1843
1877 - // Return  
1878 - return;  
1879 -} 1844 + // Initialise column start indices to 0
  1845 + for (int col = 0; col <= m_cols; ++col)
  1846 + m_colstart[col] = 0;
1880 1847
  1848 + // Optionally allocate memory for matrix elements
  1849 + if (elements > 0)
  1850 + alloc_elements(0, elements);
1881 1851
1882 -/***************************************************************************  
1883 - * Initialise class members *  
1884 - * ----------------------------------------------------------------------- *  
1885 - ***************************************************************************/  
1886 -void GSparseMatrix::init_members(void)  
1887 -{  
1888 - // Initialise sparse matrix members  
1889 - m_rowinx = NULL;  
1890 - m_mem_block = G_SPARSE_MATRIX_DEFAULT_MEM_BLOCK;  
1891 - m_zero = 0.0;  
1892 - m_fill_val = 0.0;  
1893 - m_fill_row = 0;  
1894 - m_fill_col = 0;  
1895 - m_symbolic = NULL;  
1896 - m_numeric = NULL;  
1897 -  
1898 - // Initialise stack members  
1899 - init_stack_members(); 1852 + } // endif: matrix was not empty
1900 1853
1901 - // Return  
1902 - return; 1854 + // Return
  1855 + return;
1903 } 1856 }
1904 1857
1905 1858
1906 -/***************************************************************************  
1907 - * Initialise fill-stack *  
1908 - * ----------------------------------------------------------------------- * 1859 +/***********************************************************************//**
  1860 + * @brief Initialise class mambers
1909 ***************************************************************************/ 1861 ***************************************************************************/
1910 -void GSparseMatrix::init_stack_members(void) 1862 +void GSparseMatrix::init_members(void)
1911 { 1863 {
1912 - // Initialise stack members  
1913 - m_stack_max_entries = 0;  
1914 - m_stack_size = 0;  
1915 - m_stack_entries = 0;  
1916 - m_stack_colinx = NULL;  
1917 - m_stack_start = NULL;  
1918 - m_stack_data = NULL;  
1919 - m_stack_rowinx = NULL;  
1920 - m_stack_work = NULL;  
1921 - m_stack_buffer = NULL;  
1922 -  
1923 - // Return  
1924 - return; 1864 + // Initialise sparse matrix members
  1865 + m_rowinx = NULL;
  1866 + m_mem_block = G_SPARSE_MATRIX_DEFAULT_MEM_BLOCK;
  1867 + m_zero = 0.0;
  1868 + m_fill_val = 0.0;
  1869 + m_fill_row = 0;
  1870 + m_fill_col = 0;
  1871 + m_symbolic = NULL;
  1872 + m_numeric = NULL;
  1873 +
  1874 + // Initialise stack members
  1875 + init_stack_members();
  1876 +
  1877 + // Return
  1878 + return;
1925 } 1879 }
1926 1880
1927 1881
1928 -/***************************************************************************  
1929 - * Copy class members *  
1930 - * ----------------------------------------------------------------------- * 1882 +/***********************************************************************//**
  1883 + * @brief Allocate matrix
  1884 + *
  1885 + * @param[in] m Matrix to be copied.
1931 ***************************************************************************/ 1886 ***************************************************************************/
1932 void GSparseMatrix::copy_members(const GSparseMatrix& m) 1887 void GSparseMatrix::copy_members(const GSparseMatrix& m)
1933 { 1888 {
1934 - // Copy GSparseMatrix members  
1935 - m_mem_block = m.m_mem_block;  
1936 - m_zero = m.m_zero;  
1937 - m_fill_val = m.m_fill_val;  
1938 - m_fill_row = m.m_fill_row;  
1939 - m_fill_col = m.m_fill_col;  
1940 -  
1941 - // Copy row indices  
1942 - if (m.m_rowinx != NULL && m.m_alloc > 0) {  
1943 - m_rowinx = new int[m.m_alloc];  
1944 - if (m_rowinx == NULL)  
1945 - throw GException::mem_alloc(G_COPY_MEMBERS, m.m_alloc);  
1946 - for (int i = 0; i < m.m_elements; ++i)  
1947 - m_rowinx[i] = m.m_rowinx[i];  
1948 - } 1889 + // Copy GSparseMatrix members
  1890 + m_mem_block = m.m_mem_block;
  1891 + m_zero = m.m_zero;
  1892 + m_fill_val = m.m_fill_val;
  1893 + m_fill_row = m.m_fill_row;
  1894 + m_fill_col = m.m_fill_col;
  1895 +
  1896 + // Copy row indices if they exist
  1897 + if (m.m_rowinx != NULL && m.m_alloc > 0) {
  1898 + m_rowinx = new int[m.m_alloc];
  1899 + for (int i = 0; i < m.m_elements; ++i)
  1900 + m_rowinx[i] = m.m_rowinx[i];
  1901 + }
1949 1902
1950 - // Copy symbolic decomposition  
1951 - if (m.m_symbolic != NULL) {  
1952 - GSparseSymbolic* symbolic = new GSparseSymbolic();  
1953 - if (symbolic == NULL)  
1954 - throw GException::mem_alloc(G_COPY_MEMBERS, 1);  
1955 - *symbolic = *((GSparseSymbolic*)m.m_symbolic);  
1956 - m_symbolic = (void*)symbolic;  
1957 - } 1903 + // Copy symbolic decomposition if it exists
  1904 + if (m.m_symbolic != NULL) {
  1905 + GSparseSymbolic* symbolic = new GSparseSymbolic();
  1906 + *symbolic = *((GSparseSymbolic*)m.m_symbolic);
  1907 + m_symbolic = (void*)symbolic;
  1908 + }
1958 1909
1959 - // Copy numeric decomposition  
1960 - if (m.m_numeric != NULL) {  
1961 - GSparseNumeric* numeric = new GSparseNumeric();  
1962 - if (numeric == NULL)  
1963 - throw GException::mem_alloc(G_COPY_MEMBERS, 1);  
1964 - *numeric = *((GSparseNumeric*)m.m_numeric);  
1965 - m_numeric = (void*)numeric;  
1966 - } 1910 + // Copy numeric decomposition if it exists
  1911 + if (m.m_numeric != NULL) {
  1912 + GSparseNumeric* numeric = new GSparseNumeric();
  1913 + *numeric = *((GSparseNumeric*)m.m_numeric);
  1914 + m_numeric = (void*)numeric;
  1915 + }
1967 1916
1968 - // Return  
1969 - return; 1917 + // Return
  1918 + return;
1970 } 1919 }
1971 1920
1972 1921
1973 -/***************************************************************************  
1974 - * Delete class members *  
1975 - * ----------------------------------------------------------------------- * 1922 +/***********************************************************************//**
  1923 + * @brief Delete class members
1976 ***************************************************************************/ 1924 ***************************************************************************/
1977 void GSparseMatrix::free_members(void) 1925 void GSparseMatrix::free_members(void)
1978 { 1926 {
1979 - // De-allocate only if memory has indeed been allocated  
1980 - if (m_numeric != NULL) delete (GSparseNumeric*)m_numeric;  
1981 - if (m_symbolic != NULL) delete (GSparseSymbolic*)m_symbolic;  
1982 - if (m_rowinx != NULL) delete [] m_rowinx; 1927 + // De-allocate only if memory has indeed been allocated
  1928 + if (m_numeric != NULL) delete (GSparseNumeric*)m_numeric;
  1929 + if (m_symbolic != NULL) delete (GSparseSymbolic*)m_symbolic;
  1930 + if (m_rowinx != NULL) delete [] m_rowinx;
1983 1931
1984 - // Properly mark members as free  
1985 - m_rowinx = NULL;  
1986 - m_symbolic = NULL;  
1987 - m_numeric = NULL; 1932 + // Properly mark members as free
  1933 + m_rowinx = NULL;
  1934 + m_symbolic = NULL;
  1935 + m_numeric = NULL;
1988 1936
1989 - // Free stack members  
1990 - free_stack_members(); 1937 + // Free stack members
  1938 + free_stack_members();
1991 1939
1992 - // Return  
1993 - return; 1940 + // Return
  1941 + return;
1994 } 1942 }
1995 1943
1996 1944
1997 -/***************************************************************************  
1998 - * Delete fill-stack class members *  
1999 - * ----------------------------------------------------------------------- * 1945 +/***********************************************************************//**
  1946 + * @brief Initialise fill stack
  1947 + *
  1948 + * The fill stack is used to fill the sparse matrix without any prior know-
  1949 + * ledge about the number of location of the non-zero matrix elements.
  1950 + * Values are first filled into the stack and flushed into the sparse
  1951 + * matrix once the stack is full.
2000 ***************************************************************************/ 1952 ***************************************************************************/
2001 -void GSparseMatrix::free_stack_members(void) 1953 +void GSparseMatrix::init_stack_members(void)
2002 { 1954 {
2003 - // Free stack members  
2004 - if (m_stack_colinx != NULL) delete [] m_stack_colinx;  
2005 - if (m_stack_start != NULL) delete [] m_stack_start;  
2006 - if (m_stack_data != NULL) delete [] m_stack_data;  
2007 - if (m_stack_rowinx != NULL) delete [] m_stack_rowinx;  
2008 - if (m_stack_work != NULL) delete [] m_stack_work;  
2009 - if (m_stack_buffer != NULL) delete [] m_stack_buffer;  
2010 -  
2011 - // Properly mark members as free  
2012 - m_stack_colinx = NULL;  
2013 - m_stack_start = NULL;  
2014 - m_stack_data = NULL;  
2015 - m_stack_rowinx = NULL;  
2016 - m_stack_work = NULL;  
2017 - m_stack_buffer = NULL; 1955 + // Initialise stack members
  1956 + m_stack_max_entries = 0;
  1957 + m_stack_size = 0;
  1958 + m_stack_entries = 0;
  1959 + m_stack_colinx = NULL;
  1960 + m_stack_start = NULL;
  1961 + m_stack_data = NULL;
  1962 + m_stack_rowinx = NULL;
  1963 + m_stack_work = NULL;
  1964 + m_stack_buffer = NULL;
  1965 +
  1966 + // Return
  1967 + return;
  1968 +}
2018 1969
2019 - // Return  
2020 - return; 1970 +
  1971 +/***********************************************************************//**
  1972 + * @brief Delete fill-stack class members
  1973 + *
  1974 + * Deletes all memory that has been allocated using the init_stack_members()
  1975 + * method.
  1976 + ***************************************************************************/
  1977 +void GSparseMatrix::free_stack_members(void)
  1978 +{
  1979 + // Free stack members
  1980 + if (m_stack_colinx != NULL) delete [] m_stack_colinx;
  1981 + if (m_stack_start != NULL) delete [] m_stack_start;
  1982 + if (m_stack_data != NULL) delete [] m_stack_data;
  1983 + if (m_stack_rowinx != NULL) delete [] m_stack_rowinx;
  1984 + if (m_stack_work != NULL) delete [] m_stack_work;
  1985 + if (m_stack_buffer != NULL) delete [] m_stack_buffer;
  1986 +
  1987 + // Properly mark members as free
  1988 + m_stack_colinx = NULL;
  1989 + m_stack_start = NULL;
  1990 + m_stack_data = NULL;
  1991 + m_stack_rowinx = NULL;
  1992 + m_stack_work = NULL;
  1993 + m_stack_buffer = NULL;
  1994 +
  1995 + // Return
  1996 + return;
2021 } 1997 }
2022 1998
2023 1999
@@ -2704,75 +2680,81 @@ void GSparseMatrix::mix_column(const double* src1_data, const int* src1_row, @@ -2704,75 +2680,81 @@ void GSparseMatrix::mix_column(const double* src1_data, const int* src1_row,
2704 = = 2680 = =
2705 ==========================================================================*/ 2681 ==========================================================================*/
2706 2682
2707 -/***************************************************************************  
2708 - * Output operator *  
2709 - * ----------------------------------------------------------------------- * 2683 +/***********************************************************************//**
  2684 + * @brief Output operator
  2685 + *
  2686 + * @param[in] os Output stream.
  2687 + * @param[in] m Matrix to put in output stream.
2710 ***************************************************************************/ 2688 ***************************************************************************/
2711 -ostream& operator<< (ostream& os, const GSparseMatrix& m) 2689 +std::ostream& operator<< (std::ostream& os, const GSparseMatrix& m)
2712 { 2690 {
2713 - // Put header in stream  
2714 - os << "=== GSparseMatrix ===" << endl;  
2715 - if (m.m_rowsel != NULL)  
2716 - os << " Number of rows ............: " << m.m_rows << " (compressed " <<  
2717 - m.m_num_rowsel << ")" << endl;  
2718 - else  
2719 - os << " Number of rows ............: " << m.m_rows << endl;  
2720 - if (m.m_colsel != NULL)  
2721 - os << " Number of columns .........: " << m.m_cols << " (compressed " <<  
2722 - m.m_num_colsel << ")" << endl;  
2723 - else  
2724 - os << " Number of columns .........: " << m.m_cols << endl;  
2725 - if (m.m_fill_val == 0.0)  
2726 - os << " Number of non-zero elements: " << m.m_colstart[m.m_cols] << " (" <<  
2727 - m.m_elements << ")" << endl;  
2728 - else {  
2729 - os << " Number of non-zero elements: " << m.m_colstart[m.m_cols]+1 << " (" <<  
2730 - m.m_elements+1 << ")" << endl;  
2731 - os << " Pending element ...........: (" << m.m_fill_row << "," << m.m_fill_col <<  
2732 - ")=" << m.m_fill_val << endl;  
2733 - }  
2734 - os << " Number of allocated cells .: " << m.m_alloc << endl;  
2735 - os << " Memory block size .........: " << m.m_mem_block << endl;  
2736 - os << " Sparse matrix fill ........: " << m.fill() << endl;  
2737 -  
2738 - // Dump elements and compression schemes  
2739 - m.dump_elements(os);  
2740 - m.dump_row_comp(os);  
2741 - m.dump_col_comp(os);  
2742 -  
2743 - // If there is a symbolic decomposition then put it also in the stream  
2744 - if (m.m_symbolic != NULL) {  
2745 - os << endl << *((GSparseSymbolic*)m.m_symbolic);  
2746 - } 2691 + // Determine number of elements
  2692 + int nonzero;
  2693 + int elements = m.m_elements;
  2694 + if (m.m_fill_val == 0.0) {
  2695 + nonzero = (m.m_colstart != NULL) ? m.m_colstart[m.m_cols] : 0;
  2696 + }
  2697 + else {
  2698 + nonzero = (m.m_colstart != NULL) ? m.m_colstart[m.m_cols]+1 : 0;
  2699 + elements++;
  2700 + }
2747 2701
2748 - // If there is a numeric decomposition then put it also in the stream  
2749 - if (m.m_numeric != NULL) {  
2750 - os << endl << *((GSparseNumeric*)m.m_numeric);  
2751 - } 2702 + // Put header in stream
  2703 + os << "=== GSparseMatrix ===" << std::endl;
  2704 + if (m.m_rowsel != NULL) {
  2705 + os << " Number of rows ............: " << m.m_rows
  2706 + << " (compressed " << m.m_num_rowsel << ")" << std::endl;
  2707 + }
  2708 + else
  2709 + os << " Number of rows ............: " << m.m_rows << std::endl;
  2710 + if (m.m_colsel != NULL) {
  2711 + os << " Number of columns .........: " << m.m_cols
  2712 + << " (compressed " << m.m_num_colsel << ")" << std::endl;
  2713 + }
  2714 + else
  2715 + os << " Number of columns .........: " << m.m_cols << std::endl;
  2716 + os << " Number of non-zero elements: " << nonzero
  2717 + << " (" << elements << ")" << std::endl;
  2718 + if (m.m_fill_val != 0.0) {
  2719 + os << " Pending element ...........: (" << m.m_fill_row << ","
  2720 + << m.m_fill_col << ")=" << m.m_fill_val << std::endl;
  2721 + }
  2722 + os << " Number of allocated cells .: " << m.m_alloc << std::endl;
  2723 + os << " Memory block size .........: " << m.m_mem_block << std::endl;
  2724 + os << " Sparse matrix fill ........: " << m.fill() << std::endl;
  2725 +
  2726 + // If there is a symbolic decomposition then put it also in the stream
  2727 + if (m.m_symbolic != NULL)
  2728 + os << std::endl << *((GSparseSymbolic*)m.m_symbolic);
  2729 +
  2730 + // If there is a numeric decomposition then put it also in the stream
  2731 + if (m.m_numeric != NULL)
  2732 + os << std::endl << *((GSparseNumeric*)m.m_numeric);
2752 2733
2753 - // Return output stream  
2754 - return os; 2734 + // Return output stream
  2735 + return os;
2755 } 2736 }
2756 2737
2757 2738
2758 -/***************************************************************************  
2759 - * Matrix absolute values *  
2760 - * ----------------------------------------------------------------------- * 2739 +/***********************************************************************//**
  2740 + * @brief Return matrix with absolute values of all elements
  2741 + *
  2742 + * @param[in] m Matrix for which absolute values are to be returned.
2761 ***************************************************************************/ 2743 ***************************************************************************/
2762 GSparseMatrix fabs(const GSparseMatrix& m) 2744 GSparseMatrix fabs(const GSparseMatrix& m)
2763 { 2745 {
2764 - // Define result matrix  
2765 - GSparseMatrix result = m; 2746 + // Define result matrix
  2747 + GSparseMatrix result = m;
2766 2748
2767 - // Fill pending element  
2768 - result.fill_pending(); 2749 + // Fill pending element
  2750 + result.fill_pending();
2769 2751
2770 - // Convert all elements to absolute values  
2771 - for (int i = 0; i < result.m_elements; ++i)  
2772 - result.m_data[i] = fabs(result.m_data[i]); 2752 + // Convert all elements to absolute values
  2753 + for (int i = 0; i < result.m_elements; ++i)
  2754 + result.m_data[i] = fabs(result.m_data[i]);
2773 2755
2774 - // Return result  
2775 - return result; 2756 + // Return result
  2757 + return result;
2776 } 2758 }
2777 2759
2778 2760