GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GSparseNumeric.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GSparseNumeric.cpp - Sparse matrix numeric analysis class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2006-2013 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GSparseNumeric.cpp
23  * @brief Sparse matrix numeric analysis class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GSparseNumeric.hpp"
32 
33 /* __ Macros _____________________________________________________________ */
34 #define CS_FLIP(i) (-(i)-2)
35 #define CS_MARK(w,j) { w [j] = CS_FLIP (w [j]) ; }
36 #define CS_MARKED(w,j) (w [j] < 0)
37 
38 
39 /*==========================================================================
40  = =
41  = Constructors/destructors =
42  = =
43  ==========================================================================*/
44 
45 /***************************************************************************
46  * GSparseNumeric constructor *
47  ***************************************************************************/
49 {
50  // Initialise private members for clean destruction
51  m_L = NULL;
52  m_U = NULL;
53  m_pinv = NULL;
54  m_B = NULL;
55  m_n_pinv = 0;
56  m_n_B = 0;
57 
58  // Return
59  return;
60 }
61 
62 
63 /***************************************************************************
64  * GSparseNumeric destructor *
65  ***************************************************************************/
67 {
68  // De-allocate only if memory has indeed been allocated
69  if (m_L != NULL) delete m_L;
70  if (m_U != NULL) delete m_U;
71  if (m_pinv != NULL) delete [] m_pinv;
72  if (m_B != NULL) delete [] m_B;
73 
74  // Return
75  return;
76 }
77 
78 
79 /*==========================================================================
80  = =
81  = GSparseNumeric operators =
82  = =
83  ==========================================================================*/
84 
85 /***************************************************************************
86  * GSparseNumeric assignment operator *
87  ***************************************************************************/
89 {
90  // Execute only if object is not identical
91  if (this != &n) {
92 
93  // De-allocate only if memory has indeed been allocated
94  if (m_L != NULL) delete m_L;
95  if (m_U != NULL) delete m_U;
96  if (m_pinv != NULL) delete [] m_pinv;
97  if (m_B != NULL) delete [] m_B;
98 
99  // Initialise private members for clean destruction
100  m_L = NULL;
101  m_U = NULL;
102  m_pinv = NULL;
103  m_B = NULL;
104  m_n_pinv = 0;
105  m_n_B = 0;
106 
107  // Copy m_L if it exists
108  if (n.m_L != NULL)
109  m_L = new GMatrixSparse(*n.m_L);
110 
111  // Copy m_U if it exists
112  if (n.m_U != NULL)
113  m_U = new GMatrixSparse(*n.m_U);
114 
115  // Copy m_pinv array if it exists
116  if (n.m_pinv != NULL && n.m_n_pinv > 0) {
117  m_pinv = new int[n.m_n_pinv];
118  for (int i = 0; i < n.m_n_pinv; ++i)
119  m_pinv[i] = n.m_pinv[i];
120  m_n_pinv = n.m_n_pinv;
121  }
122 
123  // Copy m_B array if it exists
124  if (n.m_B != NULL && n.m_n_B > 0) {
125  m_B = new double[n.m_n_B];
126  for (int i = 0; i < n.m_n_B; ++i)
127  m_B[i] = n.m_B[i];
128  m_n_B = n.m_n_B;
129  }
130 
131  } // endif: object was not identical
132 
133  // Return
134  return *this;
135 }
136 
137 
138 /*==========================================================================
139  = =
140  = GSparseNumeric member functions =
141  = =
142  ==========================================================================*/
143 
144 /***************************************************************************
145  * cholesky_numeric_analysis *
146  * ----------------------------------------------------------------------- *
147  * L = chol (A, [pinv parent cp]), pinv is optional. This function *
148  * allocates the memory to hold the information. *
149  * ----------------------------------------------------------------------- *
150  * Input: A Sparse matrix *
151  * S Symbolic analysis of sparse matrix *
152  ***************************************************************************/
154  const GSparseSymbolic& S)
155 {
156  // De-allocate memory that has indeed been previously allocated
157  if (m_L != NULL) delete m_L;
158  if (m_U != NULL) delete m_U;
159  if (m_pinv != NULL) delete [] m_pinv;
160  if (m_B != NULL) delete [] m_B;
161 
162  // Initialise members
163  m_L = NULL;
164  m_U = NULL;
165  m_pinv = NULL;
166  m_B = NULL;
167  m_n_pinv = 0;
168  m_n_B = 0;
169 
170  // Return if arrays in the symbolic analysis have not been allocated
171  if (!S.m_cp || !S.m_parent) return;
172 
173  // Declare
174  double lki;
175  int top, i, p, k;
176 
177  // Assign input matrix attributes
178  int n = A.m_cols;
179 
180  // Allocate int workspace
181  int wrk_size = 2*n;
182  int* wrk_int = new int[wrk_size];
183 
184  // Allocate double workspace
185  wrk_size = n;
186  double* wrk_double = new double[wrk_size];
187 
188  // Assign pointers
189  int* cp = S.m_cp;
190  int* pinv = S.m_pinv;
191  int* parent = S.m_parent;
192 
193  // Assign C = A(p,p) where A and C are symmetric and the upper part stored
194  GMatrixSparse C = (pinv) ? cs_symperm(A, pinv) : (A);
195 
196  // Assign workspace pointer
197  int* c = wrk_int;
198  int* s = wrk_int + n;
199  double* x = wrk_double;
200 
201  // Assign C matrix pointers
202  int* Cp = C.m_colstart;
203  int* Ci = C.m_rowinx;
204  double* Cx = C.m_data;
205 
206  // Allocate L matrix
207  m_L = new GMatrixSparse(n, n, cp[n]);
208 
209  // Assign L matrix pointers
210  int* Lp = m_L->m_colstart;
211  int* Li = m_L->m_rowinx;
212  double* Lx = m_L->m_data;
213 
214  // Initialise column pointers of L and c
215  for (k = 0; k < n; ++k)
216  Lp[k] = c[k] = cp[k];
217 
218  // Compute L(:,k) for L*L' = C
219  for (k = 0; k < n; k++) {
220 
221  // Nonzero pattern of L(k,:).
222  // Returns -1 if parent = NULL, s = NULL or c = NULL
223  top = cs_ereach(&C, k, parent, s, c); // find pattern of L(k,:)
224  x[k] = 0; // x (0:k) is now zero
225 
226  // x = full(triu(C(:,k)))
227  for (p = Cp[k]; p < Cp[k+1]; p++) {
228  if (Ci[p] <= k)
229  x[Ci[p]] = Cx[p];
230  }
231 
232  // d = C(k,k)
233  double d = x[k];
234 
235  // Clear x for k+1st iteration
236  x[k] = 0;
237 
238  // Triangular solve: Solve L(0:k-1,0:k-1) * x = C(:,k)
239  for ( ; top < n; top++) {
240  i = s[top]; // s [top..n-1] is pattern of L(k,:)
241  lki = x[i]/Lx[Lp[i]]; // L(k,i) = x (i) / L(i,i)
242  x[i] = 0; // clear x for k+1st iteration
243  for (p = Lp[i]+1; p < c[i]; p++)
244  x[Li[p]] -= Lx[p] * lki;
245  d -= lki * lki; // d = d - L(k,i)*L(k,i)
246  p = c[i]++;
247  Li[p] = k; // store L(k,i) in column i
248  Lx[p] = lki;
249  }
250 
251  // Compute L(k,k)
252  if (d <= 0)
254  "GSparseNumeric::cholesky_numeric_analysis(GMatrixSparse&, const GSparseSymbolic&)",
255  k, d);
256 
257  // Store L(k,k) = sqrt (d) in column k
258  p = c[k]++;
259  Li[p] = k;
260  Lx[p] = sqrt(d);
261  }
262 
263  // Finalize L
264  Lp[n] = cp[n];
265 
266  // Free workspace
267  delete [] wrk_int;
268  delete [] wrk_double;
269 
270  // Return void
271  return;
272 }
273 
274 
275 /*==========================================================================
276  = =
277  = GSparseNumeric private functions =
278  = =
279  ==========================================================================*/
280 
281 /***************************************************************************
282  * cs_ereach *
283  * ----------------------------------------------------------------------- *
284  * Find nonzero pattern of Cholesky L(k,1:k-1) using etree and *
285  * triu(A(:,k)) *
286  * ----------------------------------------------------------------------- *
287  * Input: A Sparse matrix *
288  * k Node *
289  * parent ... *
290  * s ... *
291  * w ... *
292  * Input: top Index, where s[top..n-1] contains pattern *
293  * of L(k,:); -1 if arrays have not been *
294  * allocated *
295  ***************************************************************************/
297  const int* parent, int* s, int* w)
298 {
299  // Return -1 if arrays have not been allocated
300  if (!parent || !s || !w) return (-1);
301 
302  // Assign A matrix attributes and pointers
303  int n = A->m_cols;
304  int top = n;
305  int* Ap = A->m_colstart;
306  int* Ai = A->m_rowinx;
307 
308  // Mark node k as visited
309  CS_MARK(w, k);
310 
311  // Loop over elements of node
312  int i;
313  for (int p = Ap[k]; p < Ap[k+1]; ++p) {
314  i = Ai[p]; // A(i,k) is nonzero
315  if (i > k) continue; // only use upper triangular part of A
316 
317  // Traverse up etree
318  int len;
319  for (len = 0; !CS_MARKED(w,i); i = parent[i]) {
320  s[len++] = i; // L(k,i) is nonzero
321  CS_MARK(w, i); // mark i as visited
322  }
323 
324  // Push path onto stack
325  while (len > 0)
326  s[--top] = s[--len];
327  }
328 
329  // Unmark all nodes
330  for (int p = top; p < n; ++p)
331  CS_MARK(w, s[p]);
332 
333  // Unmark node k
334  CS_MARK(w, k);
335 
336  // Return index top, where s[top..n-1] contains pattern of L(k,:)
337  return top;
338 }
339 
340 
341 /*==========================================================================
342  = =
343  = GSparseNumeric static member functions =
344  = =
345  ==========================================================================*/
346 
347 /*==========================================================================
348  = =
349  = GSparseNumeric friends =
350  = =
351  ==========================================================================*/
352 
353 /***************************************************************************
354  * Output operator *
355  ***************************************************************************/
356 std::ostream& operator<< (std::ostream& os, const GSparseNumeric& n)
357 {
358  // Put header is stream
359  os << "=== GSparseNumeric ===";
360 
361  // Show L or V matrix
362  if (n.m_L != NULL) {
363  os << std::endl <<
364  " === L or V ";
365  os << *n.m_L << std::endl;
366  }
367 
368  // Show U or R matrix
369  if (n.m_U != NULL) {
370  os << std::endl <<
371  " === U or R matrix ";
372  os << *n.m_U << std::endl;
373  }
374 
375  // Show inverse permutation if it exists
376  if (n.m_pinv != NULL) {
377  os << std::endl << " Inverse permutation (of " << n.m_n_pinv << " elements)" << std::endl;
378  for (int i = 0; i < n.m_n_pinv; ++i)
379  os << " " << n.m_pinv[i];
380  }
381 
382  // Show beta array if it exists
383  if (n.m_B != NULL) {
384  os << std::endl << " Beta[0.." << n.m_n_B << "]" << std::endl;
385  for (int i = 0; i < n.m_n_B; ++i)
386  os << " " << n.m_B[i];
387  }
388 
389  // Return output stream
390  return os;
391 }
392 
393 /*==========================================================================
394  = =
395  = GSparseNumeric exception classes =
396  = =
397  ==========================================================================*/
GMatrixSparse cs_symperm(const GMatrixSparse &matrix, const int *pinv)
cs_symperm
GMatrixSparse * m_L
Sparse matrix numeric analysis class definition.
int * m_parent
elimination tree for Cholesky and QR
int * m_cp
column pointers for Cholesky, row counts for QR
Sparse matrix class interface definition.
#define CS_MARKED(w, j)
GMatrixSparse * m_U
#define CS_MARK(w, j)
Sparse matrix symbolic analysis class.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
double * m_data
Matrix data.
int * m_pinv
Inverse row permutation for QR, fill reduce permutation for Cholesky.
int m_cols
Number of columns.
std::ostream & operator<<(std::ostream &os, const GBase &base)
Output operator.
Definition: GBase.cpp:57
int * m_rowinx
Row-indices of all elements.
Sparse matrix numeric analysis class.
GSparseNumeric & operator=(const GSparseNumeric &n)
int cs_ereach(const GMatrixSparse *A, int k, const int *parent, int *s, int *w)
virtual ~GSparseNumeric(void)
friend class GMatrixSparse
void cholesky_numeric_analysis(const GMatrixSparse &m, const GSparseSymbolic &s)
int * m_colstart
Column start indices (m_cols+1)