Luigi Tibaldo / gammalib

Commit 06c2bd6316a2697ae2135e5cce0e2f5cf652e5f2

Authored by Jürgen Knödlseder
1 parent 44eb3c38

Removed

Inline Side-by-side
Showing 1 changed file with 0 additions and 693 deletions
doc/source/user_manual/using.rst deleted
View file @44eb3c3
1   -Users manual
2   -============
3   -
4   -This section provides a detailed list ...
5   -
6   -Linear algebra
7   ---------------
8   -
9   -Vectors
10   -~~~~~~~
11   -
12   -General
13   -^^^^^^^
14   -
15   -A vector is a one-dimensional array of successive double type values.
16   -Vectors are handled in  by GVector objects. On construction, the
17   -dimension of the vector has to be specified. In other words
18   -
19   -.. code-block:: cpp
20   -
21   - GVector vector; // WRONG: constructor needs dimension
22   -
23   -is not allowed. The minimum dimension of a vector is 1, i.e. there is no
24   -such thing like an empty vector:
25   -
26   -.. code-block:: cpp
27   -
28   - GVector vector(0); // WRONG: empty vector not allowed
29   -
30   -The correct allocation of a vector is done using
31   -
32   -.. code-block:: cpp
33   -
34   - GVector vector(10); // Allocates a vector with 10 elements
35   -
36   -On allocation, all elements of a vector are set to 0. Vectors may also
37   -be allocated by copying from another vector
38   -
39   -.. code-block:: cpp
40   -
41   - GVector vector(10); // Allocates a vector with 10 elements
42   - GVector another = vector; // Allocates another vector with 10 elements
43   -
44   -or by using
45   -
46   -.. code-block:: cpp
47   -
48   - GVector vector = GVector(10); // Allocates a vector with 10 elements
49   -
50   -Vector elements are accessed using the ``( )`` operator:
51   -
52   -.. code-block:: cpp
53   -
54   - GVector vector(10); // Allocates a vector with 10 elements
55   - for (int i = 0; i < 10; ++i)
56   - vector(i) = (i+1)*10.0; // Set elements 10, 20, ..., 100
57   - for (int i = 0; i < 10; ++i)
58   - cout << vector(i) << endl; // Dump all elements, one by row
59   -
60   -The content of a vector may also be dumped using
61   -
62   -.. code-block:: cpp
63   -
64   - cout << vector << endl; // Dump entire vector
65   -
66   -which in the above example will put the sequence ::
67   -
68   - 10 20 30 40 50 60 70 80 90 100
69   -
70   -on the screen.
71   -
72   -Vector arithmetics
73   -^^^^^^^^^^^^^^^^^^
74   -
75   -Vectors can be very much handled like double type variables with the
76   -difference that operations are performed on each element of the vector.
77   -The complete list of fundamental vector operators is:
78   -
79   -.. code-block:: cpp
80   -
81   - c = a + b; // Vector + Vector addition
82   - c = a + s; // Vector + Scalar addition
83   - c = s + b; // Scalar + Vector addition
84   - c = a - b; // Vector - Vector subtraction
85   - c = a - s; // Vector - Scalar subtraction
86   - c = s - b; // Scalar - Vector subtraction
87   - s = a * b; // Vector * Vector multiplication (dot product)
88   - c = a * s; // Vector * Scalar multiplication
89   - c = s * b; // Scalar * Vector multiplication
90   - c = a / s; // Vector * Scalar division
91   -
92   -where ``a``, ``b`` and ``c`` are of type ``GVector`` and ``s`` is of type ``double``. Note in
93   -particular the combination of ``GVector`` and double type objects in
94   -addition, subtraction, multiplication and division. In these cases the
95   -specified operation is applied to each of the vector elements. It is
96   -also obvious that only vector of identicial dimension can occur in
97   -vector operations. Dimension errors can be catched by the try - catch
98   -functionality
99   -
100   -.. code-block:: cpp
101   -
102   - try {
103   - GVector a(10);
104   - GVector b(11);
105   - GVector c = a + b; // WRONG: Vectors have incompatible dimensions
106   - }
107   - catch (GVector::vector_mismatch &e) {
108   - cout << e.what() << endl; // Dimension exception is catched here
109   - throw;
110   - }
111   -
112   -Further vector operations are
113   -
114   -.. code-block:: cpp
115   -
116   - c = a; // Vector assignment
117   - c = s; // Scalar assignment
118   - s = c(index); // Vector element access
119   - c += a; // c = c + a;
120   - c -= a; // c = c - a;
121   - c += s; // c = c + s;
122   - c -= s; // c = c - s;
123   - c *= s; // c = c * s;
124   - c /= s; // c = c / s;
125   - c = -a; // Vector negation
126   -
127   -Finally, the comparison operators
128   -
129   -.. code-block:: cpp
130   -
131   - int equal = (a == b); // True if all elements equal
132   - int unequal = (a != b); // True if at least one elements unequal
133   -
134   -allow to compare all elements of a vector. If all elements are
135   -identical, the ``==`` operator returns ``true``, otherwise ``false``. If at least
136   -one element differs, the ``!=`` operator returns ``true``, is all elements are
137   -identical it returns ``false``.
138   -
139   -In addition to the operators, the following mathematical functions can
140   -be applied to vectors::
141   -
142   - acos atan exp sin tanh
143   - acosh atanh fabs sinh
144   - asin cos log sqrt
145   - asinh cosh log10 tan
146   -
147   -Again, these functions should be understood to be applied element wise.
148   -They all take a vector as argument and produce a vector as result. For
149   -example
150   -
151   -.. code-block:: cpp
152   -
153   - c = sin(a);
154   -
155   -attributes the sine of each element of vector a to vector ``c``. Additional
156   -implemented functions are
157   -
158   -.. code-block:: cpp
159   -
160   - c = cross(a, b); // Vector cross product (for 3d only)
161   - s = norm(a); // Vector norm |a|
162   - s = min(a); // Minimum element of vector
163   - s = max(a); // Maximum element of vector
164   - s = sum(a); // Sum of vector elements
165   -
166   -Finally, a small number of vector methods have been implemented:
167   -
168   -.. code-block:: cpp
169   -
170   - int n = a.size(); // Returns dimension of vector
171   - int n = a.non_zeros(); // Returns number of non-zero elements in vector
172   -
173   -Matrixes
174   -~~~~~~~~
175   -
176   -General
177   -^^^^^^^
178   -
179   -A matrix is a two-dimensional array of double type values, arranged in
180   -rows and columns. Matrixes are handled in  by ``GMatrix`` objects and the
181   -derived classes ``GSymMatrix`` and ``GSparseMatrix`` (see section
182   -[sec:matrix:storage]). On construction, the dimension of the matrix has
183   -to be specified
184   -
185   -.. code-block:: cpp
186   -
187   - GMatrix matrix(10,20); // Allocates 10 rows and 20 columns
188   -
189   -Similar to vectors, there is no such thing as a matrix without
190   -dimensions in GammaLib.
191   -
192   -Matrix storage classes
193   -^^^^^^^^^^^^^^^^^^^^^^
194   -
195   -[sec:matrix:storage]
196   -
197   -In the most general case, the rows and columns of a matrix are stored in
198   -a continuous array of :math:`{\tt rows} \times {\tt columns}` memory
199   -locations. This storage type is referred to as a *full matrix*, and is
200   -implemented by the class ``GMatrix``. Operations on full matrixes are in
201   -general relatively fast, but memory requirements may be important to
202   -hold all the elements. In general matrixes are stored by  column-wise
203   -(or in column-major format). For example, the matrix ::
204   -
205   - 1 2 3 4 5
206   - 6 7 8 9 10
207   - 11 12 13 14 15
208   -
209   -is stored in memory as ::
210   -
211   - | 1 6 11 | 2 7 12 | 3 8 13 | 4 9 14 | 5 10 15 |
212   -
213   -Many physical or mathematical problems treat with a subclass of matrixes
214   -that is symmetric, i.e. for which the element (row,col) is identical to
215   -the element (col,row). In this case, the duplicated elements need not to
216   -be stored. The derived class ``GSymMatrix`` implements such a storage type.
217   -``GSymMatrix`` stores the lower-left triangle of the matrix in column-major
218   -format. For illustration, the matrix ::
219   -
220   - 1 2 3 4
221   - 2 5 6 7
222   - 3 6 8 9
223   - 4 7 9 10
224   -
225   -is stored in memory as ::
226   -
227   - | 1 2 3 4 | 5 6 7 | 8 9 | 10 |
228   -
229   -This divides the storage requirements to hold the matrix elements by
230   -almost a factor of two.
231   -
232   -Finally, quite often one has to deal with matrixes that contain a large
233   -number of zeros. Such matrixes are called *sparse matrixes*. If only the
234   -non-zero elements of a sparse matrix are stored the memory requirements
235   -are considerably reduced. This goes however at the expense of matrix
236   -element access, which has become now more complex. In particular,
237   -filling efficiently a sparse matrix is a non-trivial problem (see
238   -section [sec:matrix:filling]). Sparse matrix storage is implemented in
239   - by the derived class ``GSparseMatrix``. A ``GSparseMatrix`` object contains
240   -three one-dimensional arrays to store the matrix elements: a double type
241   -array that contains in continuous column-major order all non-zero
242   -elements, an int type array that contains for each non-zero element the
243   -row number of its location, and an int type array that contains the
244   -storage location of the first non-zero element for each matrix column.
245   -To illustrate this storage format, the matrix ::
246   -
247   - 1 0 0 7
248   - 2 5 0 0
249   - 3 0 6 0
250   - 4 0 0 8
251   -
252   -is stored in memory as ::
253   -
254   - | 1 2 3 4 | 5 | 6 | 7 8 | Matrix elements
255   - | 0 1 2 3 | 1 | 2 | 0 3 | Row indices for all elements
256   - | 0 | 4 | 5 | 6 | Storage location of first element of each column
257   -
258   -This example is of course not very economic, since the total number of
259   -Bytes used to store the matrix is
260   -:math:`8 \times 8 + (8 + 4) \times 4 = 112` Bytes, while a full
261   -:math:`4 \times 4` matrix is stored in
262   -:math:`(4 \times 4) \times 8 = 128` Bytes (recall: a double type values
263   -takes 8 Bytes, an int type value takes 4 Bytes). For realistic large
264   -systems, however, the gain in memory space can be dramatical.
265   -
266   -The usage of the ``GMatrix``, ``GSymMatrix`` and ``GSparseMatrix`` classes is
267   -analoguous in that they implement basically all functions and methods in
268   -an identical way. So from the semantics the user has not to worry about
269   -the storage class. However, matrix element access speeds are not
270   -identical for all storage types, and if performance is an issue (as it
271   -certainly always will be), the user has to consider matrix access more
272   -carefully (see section [sec:matrix:filling]).
273   -
274   -Matrix allocation is performed using the constructors:
275   -
276   -.. code-block:: cpp
277   -
278   - GMatrix A(10,20); // Full 10 x 20 matrix
279   - GSymMatrix B(10,10); // Symmetric 10 x 10 matrix
280   - GSparseMatrix C(1000,10000); // Sparse 1000 x 10000 matrix
281   -
282   - GMatrix A(0,0); // WRONG: empty matrix not allowed
283   - GSymMatrix B(20,22); // WRONG: symmetric matrix requested
284   -
285   -In the constructor, the first argument specifies the number of rows, the
286   -second the number of columns: ``A(row,column)``. A symmetric matrix needs of
287   -course an equal number of rows and columns. And an empty matrix is not
288   -allowed. All matrix elements are initialised to 0 by the matrix
289   -allocation.
290   -
291   -Matrix elements are accessed by the ``A(row,col)`` function, where row and
292   -col start from 0 for the first row or column and run up to the number of
293   -rows or columns minus 1:
294   -
295   -.. code-block:: cpp
296   -
297   - for (int row = 0; row < n_rows; ++row) {
298   - for (int col = 0; col < n_cols; ++col)
299   - A(row,col) = (row+col)/2.0; // Set value of matrix element
300   - }
301   - ...
302   - double sum2 = 0.0;
303   - for (int row = 0; row < n_rows; ++row) {
304   - for (int col = 0; col < n_cols; ++col)
305   - sum2 *= A(row,col) * A(row,col); // Get value of matrix element
306   - }
307   -
308   -The content of a matrix can be visualised using
309   -
310   -.. code-block:: cpp
311   -
312   - cout << A << endl; // Dump matrix
313   -
314   -Matrix arithmetics
315   -^^^^^^^^^^^^^^^^^^
316   -
317   -The following description of matrix arithmetics applies to all storage
318   -classes (see section [sec:matrix:storage]). The following matrix
319   -operators have been implemented in
320   -
321   -.. code-block:: cpp
322   -
323   - C = A + B; // Matrix Matrix addition
324   - C = A - B; // Matrix Matrix subtraction
325   - C = A * B; // Matrix Matrix multiplication
326   - C = A * v; // Matrix Vector multiplication
327   - C = A * s; // Matrix Scalar multiplication
328   - C = s * A; // Scalar Matrix multiplication
329   - C = A / s; // Matrix Scalar division
330   - C = -A; // Negation
331   - A += B; // Matrix inplace addition
332   - A -= B; // Matrix inplace subtraction
333   - A *= B; // Matrix inplace multiplications
334   - A *= s; // Matrix inplace scalar multiplication
335   - A /= s; // Matrix inplace scalar division
336   -
337   -The comparison operators
338   -
339   -.. code-block:: cpp
340   -
341   - int equal = (A == B); // True if all elements equal
342   - int unequal = (A != B); // True if at least one elements unequal
343   -
344   -allow to compare all elements of a matrix. If all elements are
345   -identical, the ``==`` operator returns ``true``, otherwise ``false``. If at least
346   -one element differs, the ``!=`` operator returns true, is all elements are
347   -identical it returns false.
348   -
349   -Matrix methods and functions
350   -^^^^^^^^^^^^^^^^^^^^^^^^^^^^
351   -
352   -A number of methods has been implemented to manipulate matrixes. The
353   -method
354   -
355   -.. code-block:: cpp
356   -
357   - A.clear(); // Set all elements to 0
358   -
359   -sets all elements to 0. The methods
360   -
361   -.. code-block:: cpp
362   -
363   - int rows = A.rows(); // Returns number of rows in matrix
364   - int cols = A.cols(); // Returns number of columns in matrix
365   -
366   -provide access to the matrix dimensions, the methods
367   -
368   -.. code-block:: cpp
369   -
370   - double sum = A.sum(); // Sum of all elements in matrix
371   - double min = A.min(); // Returns minimum element of matrix
372   - double max = A.max(); // Returns maximum element of matrix
373   -
374   -inform about some matrix properties. The methods
375   -
376   -.. code-block:: cpp
377   -
378   - GVector v_row = A.extract_row(row); // Puts row in vector
379   - GVector v_column = A.extract_col(col); // Puts column in vector
380   -
381   -extract entire rows and columns from a matrix. Extraction of lower or
382   -upper triangle parts of a matrix into another is performed using
383   -
384   -.. code-block:: cpp
385   -
386   - B = A.extract_lower_triangle(); // B holds lower triangle
387   - B = A.extract_upper_triangle(); // B holds upper triangle
388   -
389   -``B`` is of the same storage class as ``A``, except for the case that ``A`` is a
390   -``GSymMatrix`` object. In this case, ``B`` will be a full matrix of type
391   -``GMatrix``.
392   -
393   -The methods
394   -
395   -.. code-block:: cpp
396   -
397   - A.insert_col(v_col,col); // Puts vector in column
398   - A.add_col(v_col,col); // Add vector to column
399   -
400   -inserts or adds the elements of a vector into a matrix column. Note that
401   -no row insertion routines have been implemented (so far) since they
402   -would be less efficient (recall that all matrix types are stored in
403   -column-major format).
404   -
405   -Conversion from one storage type to another is performed using
406   -
407   -.. code-block:: cpp
408   -
409   - B = A.convert_to_full(); // Converts A -> GMatrix
410   - B = A.convert_to_sym(); // Converts A -> GSymMatrix
411   - B = A.convert_to_sparse(); // Converts A -> GSparseMatrix
412   -
413   -Note that ``convert_to_sym()`` can only be applied to a matrix that is
414   -indeed symmetric.
415   -
416   -The transpose of a matrix can be obtained by using one of
417   -
418   -.. code-block:: cpp
419   -
420   - A.transpose(); // Transpose method
421   - B = transpose(A); // Transpose function
422   -
423   -The absolute value of a matrix is provided by
424   -
425   -.. code-block:: cpp
426   -
427   - B = fabs(A); // B = |A|
428   -
429   -Matrix factorisations
430   -^^^^^^^^^^^^^^^^^^^^^
431   -
432   -A general tool of numeric matrix calculs is factorisation.
433   -
434   -Solve linear equation ``Ax = b``. Inverse a matrix (by solving successively
435   -``Ax = e``, where ``e`` are the unit vectors for all dimensions).
436   -
437   -For symmetric and positive definite matrices the most efficient
438   -factorisation is the Cholesky decomposition. The following code fragment
439   -illustrates the usage:
440   -
441   -.. code-block:: cpp
442   -
443   - GMatrix A(n_rows, n_cols};
444   - GVector x(n_rows};
445   - GVector b(n_rows};
446   - ...
447   - A.cholesky_decompose(); // Perform Cholesky factorisation
448   - x = A.cholesky_solver(b); // Solve Ax=b for x
449   -
450   -Note that once the function ``A.cholesky_decompose()`` has been applied,
451   -the original matrix content has been replaced by its Cholesky
452   -decomposition. Since the Cholesky decomposition can be performed inplace
453   -(i.e. without the allocation of additional memory to hold the result),
454   -the matrix replacement is most memory economic. In case that the
455   -original matrix should be kept, one may either copy it before into
456   -another ``GMatrix`` object or use the function
457   -
458   -.. code-block:: cpp
459   -
460   - GMatrix L = cholesky_decompose(A);
461   - x = L.cholesky_solver(b);
462   -
463   -A symmetric and positive definite matrix can be inverted using the
464   -Cholesky decomposition using
465   -
466   -.. code-block:: cpp
467   -
468   - A.cholesky_invert(); // Inverse matrix using Cholesky fact.
469   -
470   -Alternatively, the function
471   -
472   -.. code-block:: cpp
473   -
474   - GMatrix A_inv = cholesky_invert(A);
475   -
476   -may be used.
477   -
478   -The Cholesky decomposition, solver and inversion routines may also be
479   -applied to matrices that contain rows or columns that are filled by
480   -zeros. In this case the functions provide the option to (logically)
481   -compress the matrices by skipping the zero rows and columns during the
482   -calculation.
483   -
484   -For compressed matrix Cholesky factorisation, only the non-zero rows and
485   -columns have to be symmetric and positive definite. In particular, the
486   -full matrix may even be non-symmetric.
487   -
488   -Sparse matrixes
489   -^^^^^^^^^^^^^^^
490   -
491   -The only exception that does not work is
492   -
493   -.. code-block:: cpp
494   -
495   - GSparseMatrix A(10,10);
496   - A(0,0) = A(1,1) = A(2,2) = 1.0; // WRONG: Cannot assign multiple at once
497   -
498   -In this case the value ``1.0`` is only assigned to the last element, i.e.
499   -``A(2,2)``, the other elements will remain ``0``. This feature has to do with
500   -the way how the compiler translates the code and how  implements sparse
501   -matrix filling. ``GSparseMatrix`` provides a pointer for a new element to be
502   -filled. Since there is only one such *fill pointer*, only one element
503   -can be filled at once in a statement. **So it is strongly advised to
504   -avoid multiple matrix element assignment in a single row.** Better write
505   -the above code like
506   -
507   -.. code-block:: cpp
508   -
509   - GSparseMatrix A;
510   - A(0,0) = 1.0;
511   - A(1,1) = 1.0;
512   - A(2,2) = 1.0;
513   -
514   -This way, element assignment works fine.
515   -
516   -Inverting a sparse matrix produces in general a full matrix, so the
517   -inversion function should be used with caution. Note that a full matrix
518   -that is stored in sparse format takes roughly twice the memory than a
519   -normal ``GMatrix`` object. If nevertheless the inverse of a sparse matrix
520   -should be examined, it is recommended to perform the analysis
521   -column-wise
522   -
523   -.. code-block:: cpp
524   -
525   - GSparseMatrix A(rows,cols); // Allocate sparse matrix
526   - GVector unit(rows); // Allocate vector
527   - ...
528   - A.cholesky_decompose(); // Factorise matrix
529   -
530   - // Column-wise solving the matrix equation
531   - for (int col = 0; col < cols; ++col) {
532   - unit(col) = 1.0; // Set unit vector
533   - GVector x = cholesky_solver(unit); // Get column x of inverse
534   - ...
535   - unit(col) = 0.0; // Clear unit vector for next round
536   - }
537   -
538   -Filling sparse matrixes
539   -^^^^^^^^^^^^^^^^^^^^^^^
540   -
541   -[sec:matrix:filling]
542   -
543   -The filling of a sparse matrix is a tricky issue since the storage of
544   -the elements depends on their distribution in the matrix. If one would
545   -know beforehand this distribution, sparse matrix filling would be easy
546   -and fast. In general, however, the distribution is not known a priori,
547   -and matrix filling may become a quite time consuming task.
548   -
549   -If a matrix has to be filled element by element, the access through the
550   -operator
551   -
552   -.. code-block:: cpp
553   -
554   - m(row,col) = value;
555   -
556   -may be mandatory. In principle, if a new element is inserted into a
557   -matrix a new memory cell has to be allocated for this element, and other
558   -elements may be moved. Memory allocation is quite time consuming, and to
559   -reduce the overhead, ``GSparseMatrix`` can be configured to allocate memory
560   -in bunches. By default, each time more matrix memory is needed,
561   -``GSparseMatrix`` allocates 512 cells at once (or 6144 Bytes since each
562   -element requires a double and a int storage location). If this amount of
563   -memory is not adequat one may change this value by using
564   -
565   -.. code-block:: cpp
566   -
567   - m.set_mem_block(size);
568   -
569   -where size is the number of matrix elements that should be allocated at
570   -once (corresponding to a total memory of :math:`12 \times {\tt size}`
571   -Bytes).
572   -
573   -Alternatively, a matrix may be filled column-wise using the functions
574   -
575   -.. code-block:: cpp
576   -
577   - m.insert_col(vector,col); // Insert a vector in column
578   - m.add_col(vector,col); // Add content of a vector to column
579   -
580   -While ``insert_col`` sets the values of column ``col`` (deleting thus any
581   -previously existing entries), ``add_col`` adds the content of vector to all
582   -elements of column col. Using these functions is considerably more rapid
583   -than filling individual values.
584   -
585   -Still, if the matrix is big (i.e. several thousands of rows and
586   -columns), filling individual columns may still be slow. To speed-up
587   -dynamical matrix filling, an internal fill-stack has been implemented in
588   -``GSparseMatrix``. Instead of inserting values column-by-column, the columns
589   -are stored in a stack and filled into the matrix once the stack is full.
590   -This reduces the number of dynamic memory allocations to let the matrix
591   -grow as it is built. By default, the internal stack is disabled. The
592   -stack can be enabled and used as follows:
593   -
594   -.. code-block:: cpp
595   -
596   - m.stack_init(size, entries); // Initialise stack
597   - ...
598   - m.add_col(vector,col); // Add columns
599   - ...
600   - m.stack_destroy(); // Flush and destory stack
601   -
602   -The method ``stack_init`` initialises a stack with a number of size
603   -elements and a maximum of entries columns. The larger the values size
604   -and entries are chosen, the more efficient the stack works. The total
605   -amount of memory of the stack can be estimated as
606   -:math:`12 \times {\tt size} + 8 \times {\tt entries}` Bytes. If a rough
607   -estimate of the total number of non-zero elements is available it is
608   -recommended to set size to this value. As a rule of thumb, size should
609   -be at least of the dimension of either the number of rows or the number
610   -of columns of the matrix (take the maximum of both). entries is best set
611   -to the number of columns of the matrix. If memory limits are an issue
612   -smaller values may be set, but if the values are too small, the speed
613   -increase may become negligible (or stack-filling may even become slower
614   -than normal filling).
615   -
616   -Stack-filling only works with the method ``add_col``. Note also that
617   -filling sub-sequently the same column leads to stack flushing. In the
618   -code
619   -
620   -.. code-block:: cpp
621   -
622   - for (int col = 0; col < 100; ++col) {
623   - column = 0.0; // Reset column
624   - column(col) = col; // Set column
625   - m.add_col(column,col); // Add column
626   - }
627   -
628   -stack flushing occurs in each loop, and consequently, the stack-filling
629   -approach will be not very efficient (it would probably be even slover
630   -than normal filling). If successive operations are to be performed on
631   -columns, it is better to perform them before adding. The code
632   -
633   -.. code-block:: cpp
634   -
635   - column = 0.0; // Reset column
636   - for (int col = 0; col < 100; ++col)
637   - column(col) = col; // Set column
638   - m.add_col(column,col); // Add column
639   -
640   -would be far more efficient.
641   -
642   -A avoidable overhead occurs for the case that the column to be added is
643   -sparse. The vector may contain many zeros, and ``GSparseMatrix`` has to
644   -filter them out. If the sparsity of the column is known, this overhead
645   -can be avoided by directly passing a compressed array to ``add_col``:
646   -
647   -.. code-block:: cpp
648   -
649   - int number = 5; // 5 elements in array
650   - double* values = new double[number]; // Allocate values
651   - int* rows = new int[number]; // Allocate row index
652   - ...
653   - m.stack_init(size, entries); // Initialise stack
654   - ...
655   - for (int i = 0; i < number; ++i) { // Initialise array
656   - values[i] = ... // ... set values
657   - rows[i] = ... // ... set row indices
658   - }
659   - ...
660   - m.add_col(values,rows,number,col); // Add array
661   - ...
662   - m.stack_destroy(); // Flush and destory stack
663   - ...
664   - delete [] values; // Free array
665   - delete [] rows;
666   -
667   -The method ``add_col`` calls the method ``stack_push_column`` for stack
668   -filling. ``add_col`` is more general than ``stack_push_column`` in that it
669   -decides which of stack- or direct filling is more adequate. In
670   -particular, ``stack_push_column`` may refuse pushing a column onto the
671   -stack if there is not enough space. In that case, ``stack_push_column``
672   -returns a non-zero value that corresponds to the number of non-zero
673   -elements in the vector that should be added. However, it is recommended
674   -to not use ``stack_push_column`` and call instead ``add_col``.
675   -
676   -The method ``stack_destroy`` is used to flush and destroy the stack. After
677   -this call the stack memory is liberated. If the stack should be flushed
678   -without destroying it, the method ``stack_flush`` may be used:
679   -
680   -.. code-block:: cpp
681   -
682   - m.stack_init(size, entries); // Initialise stack
683   - ...
684   - m.add_col(vector,col); // Add columns
685   - ...
686   - m.stack_flush(); // Simply flush stack
687   -
688   -Once flushed, the stack can be filled anew.
689   -
690   -Note that stack flushing is not automatic! This means, if one trys to
691   -use a matrix for calculs without flushing, the calculs may be wrong.
692   -**If a stack is used for filling, always flush the stack before using
693   -the matrix.**