GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GSparseSymbolic.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GSparseSymbolic.cpp - Sparse matrix symbolic analysis class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2006-2021 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 GSparseSymbolic.cpp
23  * @brief Sparse matrix symbolic analysis class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GTools.hpp"
33 #include "GMatrixSparse.hpp"
34 #include "GSparseSymbolic.hpp"
35 
36 /* __ Method name definitions ____________________________________________ */
37 #define G_CHOLESKY "GSparseSymbolic::cholesky_symbolic_analysis(int, "\
38  "GMatrixSparse&)"
39 #define G_CS_AMD "GSparseSymbolic::cs_amd(int, GMatrixSparse*)"
40 
41 /* __ Macros _____________________________________________________________ */
42 #define CS_MIN(a,b) (((a) < (b)) ? (a) : (b))
43 #define CS_MAX(a,b) (((a) > (b)) ? (a) : (b))
44 #define CS_FLIP(i) (-(i)-2)
45 #define HEAD(k,j) (ata ? head [k] : j)
46 #define NEXT(J) (ata ? next [J] : -1)
47 
48 /* __ Definitions ________________________________________________________ */
49 //#define G_DEBUG_SPARSE_SYM_AMD_ORDERING // Debug AMD ordering
50 //#define G_DEBUG_SPARSE_CHOLESKY // Analyse Cholesky decomposition
51 
52 
53 /*==========================================================================
54  = =
55  = Constructors/destructors =
56  = =
57  ==========================================================================*/
58 
59 /***********************************************************************//**
60  * @brief Void constructor
61  ***************************************************************************/
63 {
64  // Initialise private members
65  m_n_pinv = 0;
66  m_n_q = 0;
67  m_n_parent = 0;
68  m_n_cp = 0;
69  m_n_leftmost = 0;
70  m_pinv = NULL;
71  m_q = NULL;
72  m_parent = NULL;
73  m_cp = NULL;
74  m_leftmost = NULL;
75  m_m2 = 0;
76  m_lnz = 0.0;
77  m_unz = 0.0;
78 
79  // Return
80  return;
81 }
82 
83 
84 /***********************************************************************//**
85  * @brief Destructor
86  ***************************************************************************/
88 {
89  // De-allocate only if memory has indeed been allocated
90  if (m_pinv != NULL) delete [] m_pinv;
91  if (m_q != NULL) delete [] m_q;
92  if (m_parent != NULL) delete [] m_parent;
93  if (m_cp != NULL) delete [] m_cp;
94  if (m_leftmost != NULL) delete [] m_leftmost;
95 
96  // Return
97  return;
98 }
99 
100 
101 /*==========================================================================
102  = =
103  = Operators =
104  = =
105  ==========================================================================*/
106 
107 /***************************************************************************
108  * GSparseSymbolic assignment operator *
109  ***************************************************************************/
111 {
112  // Execute only if object is not identical
113  if (this != &s) {
114 
115  // De-allocate only if memory has indeed been allocated
116  if (m_pinv != NULL) delete [] m_pinv;
117  if (m_q != NULL) delete [] m_q;
118  if (m_parent != NULL) delete [] m_parent;
119  if (m_cp != NULL) delete [] m_cp;
120  if (m_leftmost != NULL) delete [] m_leftmost;
121 
122  // Initialise private members for clean destruction
123  m_n_pinv = 0;
124  m_n_q = 0;
125  m_n_parent = 0;
126  m_n_cp = 0;
127  m_n_leftmost = 0;
128  m_pinv = NULL;
129  m_q = NULL;
130  m_parent = NULL;
131  m_cp = NULL;
132  m_leftmost = NULL;
133  m_m2 = 0;
134  m_lnz = 0.0;
135  m_unz = 0.0;
136 
137  // Copy data members
138  m_m2 = s.m_m2;
139  m_lnz = s.m_lnz;
140  m_unz = s.m_unz;
141 
142  // Copy m_pinv array if it exists
143  if (s.m_pinv != NULL && s.m_n_pinv > 0) {
144  m_pinv = new int[s.m_n_pinv];
145  for (int i = 0; i < s.m_n_pinv; ++i) {
146  m_pinv[i] = s.m_pinv[i];
147  }
148  m_n_pinv = s.m_n_pinv;
149  }
150 
151  // Copy m_q array if it exists
152  if (s.m_q != NULL && s.m_n_q > 0) {
153  m_q = new int[s.m_n_q];
154  for (int i = 0; i < s.m_n_q; ++i) {
155  m_q[i] = s.m_q[i];
156  }
157  m_n_q = s.m_n_q;
158  }
159 
160  // Copy m_parent array if it exists
161  if (s.m_parent != NULL && s.m_n_parent > 0) {
162  m_parent = new int[s.m_n_parent];
163  for (int i = 0; i < s.m_n_parent; ++i) {
164  m_parent[i] = s.m_parent[i];
165  }
167  }
168 
169  // Copy m_cp array if it exists
170  if (s.m_cp != NULL && s.m_n_cp > 0) {
171  m_cp = new int[s.m_n_cp];
172  for (int i = 0; i < s.m_n_cp; ++i) {
173  m_cp[i] = s.m_cp[i];
174  }
175  m_n_cp = s.m_n_cp;
176  }
177 
178  // Copy m_leftmost array if it exists
179  if (s.m_leftmost != NULL && s.m_n_leftmost > 0) {
180  m_leftmost = new int[s.m_n_leftmost];
181  for (int i = 0; i < s.m_n_leftmost; ++i) {
182  m_leftmost[i] = s.m_leftmost[i];
183  }
185  }
186 
187  } // endif: object was not identical
188 
189  // Return
190  return *this;
191 }
192 
193 
194 /*==========================================================================
195  = =
196  = GSparseSymbolic member functions =
197  = =
198  ==========================================================================*/
199 
200 /***************************************************************************
201  * @brief Symbolc Cholesky analysis
202  *
203  * @param[in] order Ordering type (0: natural, 1: Cholesky decomposition).
204  * @param[in] m Spare matrix.
205  *
206  * Ordering and symbolic analysis for a Cholesky factorization. This method
207  * allocates the memory to hold the information.
208  ***************************************************************************/
210  const GMatrixSparse& m)
211 {
212  // Debug
213  #if defined(G_DEBUG_SPARSE_CHOLESKY)
214  std::cout << endl << "GSparseSymbolic::cholesky_symbolic_analysis entered" << std::endl;
215  std::cout << " order ......................: " << order << std::endl;
216  std::cout << " number of rows .............: " << m.m_rows << std::endl;
217  std::cout << " number of columns ..........: " << m.m_cols << std::endl;
218  std::cout << " number of non-zero elements : " << m.m_elements << std::endl;
219  #endif
220 
221  // De-allocate memory that has indeed been previously allocated
222  if (m_pinv != NULL) delete [] m_pinv;
223  if (m_q != NULL) delete [] m_q;
224  if (m_parent != NULL) delete [] m_parent;
225  if (m_cp != NULL) delete [] m_cp;
226  if (m_leftmost != NULL) delete [] m_leftmost;
227 
228  // Initialise members
229  m_n_pinv = 0;
230  m_n_q = 0;
231  m_n_parent = 0;
232  m_n_cp = 0;
233  m_n_leftmost = 0;
234  m_pinv = NULL;
235  m_q = NULL;
236  m_parent = NULL;
237  m_cp = NULL;
238  m_leftmost = NULL;
239  m_m2 = 0;
240  m_lnz = 0.0;
241  m_unz = 0.0;
242 
243  // Check if order type is valid
244  if (order < 0 || order > 1) {
245  std::string msg = "Invalid ordering type "+gammalib::str(order)+
246  "specified. The ordering type must be comprised in "
247  "[0,1]. Please specify a valid ordering type.";
249  }
250 
251  // Assign input matrix attributes
252  int n = m.m_cols;
253 
254  // P = amd(A+A'), or natural
255  int* P = cs_amd(order, &m);
256  #if defined(G_DEBUG_SPARSE_CHOLESKY)
257  std::cout << " AMD ordering permutation (" << P << ")" << std::endl;
258  if (P != NULL) {
259  for (int i = 0; i < n; ++i) {
260  std::cout << " " << P[i];
261  }
262  std::cout << std::endl;
263  }
264  #endif
265 
266  // Find inverse permutation and store it in class member 'm_pinv'.
267  // Note that if P = NULL or n < 1 this function returns NULL.
268  m_pinv = cs_pinv(P, n);
269  m_n_pinv = n;
270  #if defined(G_DEBUG_SPARSE_CHOLESKY)
271  std::cout << " Inverse permutation (" << m_pinv << ")" << std::endl;
272  if (m_pinv != NULL) {
273  for (int i = 0; i < n; ++i) {
274  std::cout << " " << m_pinv[i];
275  }
276  std::cout << std::endl;
277  }
278  #endif
279 
280  // Delete workspace
281  if (P != NULL) {
282  delete [] P;
283  }
284 
285  // C = spones(triu(A(P,P)))
287  #if defined(G_DEBUG_SPARSE_CHOLESKY)
288  std::cout << " C = spones(triu(A(P,P))) " << C << std::endl;
289  #endif
290 
291  // Find etree of C and store it in class member 'm_parent'
292  m_parent = cs_etree(&C, 0);
293  m_n_parent = n;
294  #if defined(G_DEBUG_SPARSE_CHOLESKY)
295  std::cout << " Elimination tree (" << m_parent << ")" << std::endl;
296  if (m_parent != NULL) {
297  for (int i = 0; i < n; ++i) {
298  std::cout << " " << m_parent[i];
299  }
300  std::cout << std::endl;
301  }
302  #endif
303 
304  // Post order the etree
305  // Note that if m_parent = NULL or n < 1 this function returns NULL.
306  int* post = cs_post(m_parent, n);
307  #if defined(G_DEBUG_SPARSE_CHOLESKY)
308  std::cout << " Post ordered elimination tree (" << post << ")" << std::endl;
309  if (post != NULL) {
310  for (int i = 0; i < n; ++i) {
311  std::cout << " " << post[i];
312  }
313  std::cout << std::endl;
314  }
315  #endif
316 
317  // Find column counts of chol(C)
318  // Note that if m_parent = NULL or post = NULL this function returns NULL.
319  int* c = cs_counts(&C, m_parent, post, 0);
320  #if defined(G_DEBUG_SPARSE_CHOLESKY)
321  std::cout << " Column counts (" << c << ")" << std::endl;
322  if (c != NULL) {
323  for (int i = 0; i < n; ++i) {
324  std::cout << " " << c[i];
325  }
326  std::cout << std::endl;
327  }
328  #endif
329 
330  // Delete workspace
331  if (post != NULL) {
332  delete [] post;
333  }
334 
335  // Allocate column pointers for Cholesky decomposition
336  m_n_cp = n+1;
337  m_cp = new int[m_n_cp];
338 
339  // Find column pointers for L
340  m_unz = m_lnz = cs_cumsum(m_cp, c, n);
341  #if defined(G_DEBUG_SPARSE_CHOLESKY)
342  std::cout << " Column pointers for L (" << m_cp << ")" << std::endl;
343  if (m_cp != NULL) {
344  for (int i = 0; i < m_n_cp; ++i) {
345  std::cout << " " << m_cp[i];
346  }
347  std::cout << std::endl;
348  }
349  std::cout << " Number of non-zero elements in L: " << m_lnz << std::endl;
350  #endif
351 
352  // Delete workspace
353  if (c != NULL) {
354  delete [] c;
355  }
356 
357  // Delete symbolic analysis if it is not valid
358  if (m_lnz < 0) {
359 
360  // De-allocate memory that has indeed been previously allocated
361  if (m_pinv != NULL) delete [] m_pinv;
362  if (m_q != NULL) delete [] m_q;
363  if (m_parent != NULL) delete [] m_parent;
364  if (m_cp != NULL) delete [] m_cp;
365  if (m_leftmost != NULL) delete [] m_leftmost;
366 
367  // Initialise members
368  m_n_pinv = 0;
369  m_n_q = 0;
370  m_n_parent = 0;
371  m_n_cp = 0;
372  m_n_leftmost = 0;
373  m_pinv = NULL;
374  m_q = NULL;
375  m_parent = NULL;
376  m_cp = NULL;
377  m_leftmost = NULL;
378  m_m2 = 0;
379  m_lnz = 0.0;
380  m_unz = 0.0;
381  }
382 
383  // Debug
384  #if defined(G_DEBUG_SPARSE_CHOLESKY)
385  std::cout << "GSparseSymbolic::cholesky_symbolic_analysis finished" << std::endl;
386  #endif
387 
388  // Return
389  return;
390 }
391 
392 
393 /*==========================================================================
394  = =
395  = GSparseSymbolic private functions =
396  = =
397  ==========================================================================*/
398 
399 /***************************************************************************
400  * @brief cs_amd
401  *
402  * @param[in] order Ordering type.
403  * @param[in] A Spare Matrix.
404  * @return Integer array of n+1 elements
405  *
406  * Applies an approximate minimum degree ordering algorithm to derive the
407  * permutations that minimise the fill-in.
408  * p = amd(A+A') if symmetric is true, or amd(A'A) otherwise
409  *
410  * The following ordering types are supported:
411  * 0: natural
412  * 1: Cholesky decomposition
413  * 2: LU decomposition
414  * 3: QR decomposition
415  ***************************************************************************/
416 int* GSparseSymbolic::cs_amd(int order, const GMatrixSparse* A)
417 {
418  // Throw exception if order is invalid
419  if (order < 0 || order > 3) {
420  std::string msg = "Invalid ordering type "+gammalib::str(order)+
421  "specified. The ordering type must be comprised in "
422  "[0,3]. Please specify a valid ordering type.";
424  }
425 
426  // Debug
427  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
428  std::cout << std::endl << "GSparseSymbolic::cs_amd entered" << std::endl;
429  std::cout << " order ......................: " << order << std::endl;
430  std::cout << " number of rows .............: " << A->m_rows << std::endl;
431  std::cout << " number of columns ..........: " << A->m_cols << std::endl;
432  std::cout << " number of non-zero elements : " << A->m_elements << std::endl;
433  #endif
434 
435  // Declare
436  int dense;
437  int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
438  k2, k3, jlast, ln, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
439  ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q;
440  unsigned int h;
441 
442  // Assign input matrix attributes
443  int m = A->m_rows;
444  int n = A->m_cols;
445 
446 
447  // ======================================================================
448  // Step 1: Construct matrix C
449  // ======================================================================
450 
451  // Debug
452  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
453  std::cout << " Step 1: Construct matrix C" << std::endl;
454  #endif
455 
456  // Initialise matrix C. We do not care at this point about the actual
457  // dimension since we attribute the matrix later anyways. It's just
458  // that there is no empty constructor of a matrix, so we have to put
459  // something by default ...
460  GMatrixSparse C(m,n);
461 
462  // Get (logical) transpose of A: AT = A'
463  // NOTE: WE ONLY NEED THE LOGICAL TRANSPOSE HERE. HOWEVER, WE HAVE NO
464  // LOGICAL ADDITION SO FAR ...
465  // GMatrixSparse AT = cs_transpose(*A, 0);
466  GMatrixSparse AT = cs_transpose(*A, 1);
467 
468  // Find dense threshold
469  dense = (int)CS_MAX(16, 10 * sqrt((double)n));
470  dense = CS_MIN(n-2, dense);
471 
472  // Case A: Cholesky decomposition of a symmetric matrix: C = A + A'
473  if (order == 1 && n == m) {
474  // NOTE: HERE WE ONLY NEED LOGICAL ADDITION
475  // C = cs_add(A, AT, 0, 0);
476  C = *A + AT;
477  }
478 
479  // Case B: LU decomposition: C = A' * A with no dense rows
480  else if (order == 2) {
481 
482  // Assign A' matrix attributes
483  int* ATp = AT.m_colstart;
484  int* ATi = AT.m_rowinx;
485 
486  // Drop dense columns from AT
487  for (p2 = 0, j = 0; j < m; j++) {
488  p = ATp[j]; // column j of AT starts here
489  ATp[j] = p2; // new column j starts here
490  if (ATp[j+1] - p > dense) {
491  continue; // skip dense col j
492  }
493  for ( ; p < ATp[j+1]; p++) {
494  ATi[p2++] = ATi[p] ;
495  }
496  }
497 
498  // Finalise AT
499  ATp[m] = p2;
500 
501  // Get (logical) transpose of AT: A2 = AT'
502  // NOTE: WE ONLY NEED THE LOGICAL TRANSPOSE HERE. HOWEVER, WE HAVE NO
503  // LOGICAL MULTIPLICATION SO FAR ...
504  // GMatrixSparse A2 = cs_transpose(AT, 0);
505  GMatrixSparse A2 = cs_transpose(AT, 1);
506 
507  // C = A' * A with no dense rows
508  // NOTE: cs_multiply NOT YET IMPLEMENTED
509  // C = A2 ? cs_multiply(AT, A2) : NULL;
510  C = AT * A2;
511  }
512 
513  // Case C: QR decomposition: C = A' * A
514  else {
515  // NOTE: cs_multiply NOT YET IMPLEMENTED
516  // C = cs_multiply(AT, A);
517  C = AT * *A;
518  }
519 
520  // Debug
521  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
522  std::cout << " Matrix C " << C << std::endl;
523  #endif
524 
525  // Drop diagonal entries from C
526  cs_fkeep(&C, &cs_diag, NULL);
527 
528  // Debug
529  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
530  std::cout << " Dropped diagonal entries from matrix C " << C << std::endl;
531  #endif
532 
533  // Assign C matrix attributes
534  int* Cp = C.m_colstart;
535  int cnz = Cp[n];
536 
537  // Allocate result array
538  int* P = new int[n+1];
539 
540  // Allocate workspace
541  int wrk_size = 8*(n+1);
542  int* wrk_int = new int[wrk_size];
543 
544  // Add elbow room to C
545  int elbow_room = cnz/5 + 2*n; // Request additional # of elements
546  C.alloc_elements(cnz, elbow_room); // Appand elements
547  int nzmax = C.m_elements; // Save the maximum # of entries for garbage coll.
548  C.m_elements -= elbow_room; // Reverse element increase of alloc_elements
549 
550  // (Re-)assign C matrix attributes
551  Cp = C.m_colstart;
552  int* Ci = C.m_rowinx;
553 
554  // Debug
555  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
556  std::cout << " Added elbow room to matrix C " << C << std::endl;
557  #endif
558 
559  // Set workspace pointers
560  int* len = wrk_int;
561  int* nv = wrk_int + (n+1);
562  int* next = wrk_int + 2*(n+1);
563  int* head = wrk_int + 3*(n+1);
564  int* elen = wrk_int + 4*(n+1);
565  int* degree = wrk_int + 5*(n+1);
566  int* w = wrk_int + 6*(n+1);
567  int* hhead = wrk_int + 7*(n+1);
568  int* last = P; // use P as workspace for last
569 
570 
571  // ======================================================================
572  // Step 2: Initialize quotient graph
573  // ======================================================================
574 
575  // Debug
576  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
577  std::cout << " Step 2: Initialize quotient graph" << std::endl;
578  #endif
579 
580  // Setup array that contains for each column the number of non-zero elements
581  for (k = 0 ; k < n ; k++) {
582  len[k] = Cp[k+1] - Cp[k];
583  }
584  len[n] = 0;
585 
586  // Debug
587  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
588  std::cout << " Non-zero elements per column: ";
589  for (k = 0 ; k <= n ; k++) {
590  std::cout << " " << len[k];
591  }
592  std::cout << std::endl;
593  #endif
594 
595  // Loop over all nodes (i.e. matrix columns)
596  for (i = 0 ; i <= n ; i++) {
597  head[i] = -1; // degree list i is empty
598  last[i] = -1;
599  next[i] = -1;
600  hhead[i] = -1; // hash list i is empty
601  nv[i] = 1; // node i is just one node
602  w[i] = 1; // node i is alive
603  elen[i] = 0; // Ek of node i is empty
604  degree[i] = len[i]; // degree of node i
605  }
606 
607  // Clear w
608  mark = cs_wclear(0, 0, w, n);
609 
610  // End-point initialisations
611  elen[n] = -2; // n is a dead element
612  Cp[n] = -1; // n is a root of assembly tree
613  w[n] = 0; // n is a dead element
614 
615 
616  // ======================================================================
617  // Step 3: Initialize degree lists
618  // ======================================================================
619 
620  // Debug
621  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
622  std::cout << " Step 3: Initialize degree lists" << std::endl;
623  #endif
624 
625  // Loop over all nodes
626  for (i = 0 ; i < n ; i++) {
627 
628  // Get degree of node
629  d = degree[i];
630 
631  // Case A: Node i is empty
632  if (d == 0) {
633  elen[i] = -2; // element i is dead
634  nel++ ;
635  Cp[i] = -1; // i is a root of assembly tree
636  w[i] = 0; // i is a dead element
637  }
638 
639  // Case B: Node is dense
640  else if (d > dense) {
641  nv[i] = 0; // absorb i into element n
642  elen[i] = -1; // node i is dead
643  nel++;
644  Cp[i] = CS_FLIP(n);
645  nv[n]++;
646  }
647 
648  // Case C: Node is sparse
649  else {
650  if (head[d] != -1) last[head[d]] = i;
651  next[i] = head[d]; // put node i in degree list d
652  head[d] = i;
653  }
654 
655  } // endfor: looped over all columns
656 
657 
658  // ======================================================================
659  // Step 4: Build permutations
660  // ======================================================================
661 
662  // Debug
663  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
664  std::cout << " Step 4: Build permutations" << std::endl;
665  #endif
666 
667  // While (selecting pivots) do
668  while (nel < n) {
669 
670  // Select node of minimum approximate degree
671  for (k = -1; mindeg < n && (k = head [mindeg]) == -1; mindeg++) ;
672  if (next[k] != -1) last [next[k]] = -1;
673  head[mindeg] = next[k]; // remove k from degree list
674  elenk = elen[k]; // elenk = |Ek|
675  nvk = nv[k]; // # of nodes k represents
676  nel += nvk; // nv[k] nodes of A eliminated
677 
678  // Garbage collection. This frees the elbow rooms in matrix C
679  if (elenk > 0 && cnz + mindeg >= nzmax) {
680 
681  // Loop over all columns
682  for (j = 0; j < n; j++) {
683  // If j is a live node or element, then ...
684  if ((p = Cp[j]) >= 0) {
685  Cp[j] = Ci[p]; // save first entry of object
686  Ci[p] = CS_FLIP(j); // first entry is now CS_FLIP(j)
687  }
688  }
689 
690  // Scan all memory
691  for (q = 0, p = 0; p < cnz; ) {
692  // If we found object j, then ...
693  if ((j = CS_FLIP(Ci[p++])) >= 0) {
694  Ci[q] = Cp[j]; // restore first entry of object
695  Cp[j] = q++; // new pointer to object j
696  for (k3 = 0; k3 < len[j]-1; k3++) {
697  Ci[q++] = Ci[p++];
698  }
699  }
700  }
701  cnz = q; // Ci[cnz...nzmax-1] now free
702 
703  } // endif: garbage collection done
704 
705  // Construct new element
706  dk = 0;
707  nv[k] = -nvk; // flag k as in Lk
708  p = Cp[k];
709  pk1 = (elenk == 0) ? p : cnz; // do in place if elen[k] == 0
710  pk2 = pk1;
711  for (k1 = 1; k1 <= elenk + 1; k1++) {
712 
713  // ...
714  if (k1 > elenk) {
715  e = k; // search the nodes in k
716  pj = p; // list of nodes starts at Ci[pj]
717  ln = len[k] - elenk; // length of list of nodes in k
718  }
719  else {
720  e = Ci[p++]; // search the nodes in e
721  pj = Cp[e];
722  ln = len[e]; // length of list of nodes in e
723  }
724 
725  // ...
726  for (k2 = 1 ; k2 <= ln ; k2++) {
727  i = Ci[pj++];
728  if ((nvi = nv [i]) <= 0) { // Skip of node i is dead or seen
729  continue;
730  }
731  dk += nvi; // degree[Lk] += size of node i
732  nv[i] = -nvi; // negate nv[i] to denote i in Lk
733  Ci[pk2++] = i; // place i in Lk
734  if (next[i] != -1) {
735  last[next[i]] = last[i];
736  }
737 
738  // remove i from degree list
739  if (last[i] != -1) {
740  next[last[i]] = next[i];
741  }
742  else {
743  head[degree[i]] = next[i];
744  }
745 
746  } // endfor: k2
747 
748  // ...
749  if (e != k) {
750  Cp[e] = CS_FLIP(k); // absorb e into k
751  w[e] = 0; // e is now a dead element
752  }
753 
754  } // endfor: k1
755 
756  if (elenk != 0) {
757  cnz = pk2; // Ci [cnz...nzmax] is free
758  }
759  degree[k] = dk; // external degree of k - |Lk\i|
760  Cp[k] = pk1; // element k is in Ci[pk1..pk2-1]
761  len[k] = pk2 - pk1;
762  elen[k] = -2; // k is now an element
763 
764  // Clear w if necessary
765  mark = cs_wclear(mark, lemax, w, n);
766 
767  // Scan 1: Find set differences |Le\Lk|
768  for (pk = pk1 ; pk < pk2 ; pk++) {
769 
770  i = Ci [pk];
771  if ((eln = elen[i]) <= 0) { // skip if elen[i] empty
772  continue;
773  }
774  nvi = -nv[i]; // nv [i] was negated
775  wnvi = mark - nvi;
776 
777  // Scan Ei
778  for (p = Cp[i]; p <= Cp[i] + eln - 1; p++) {
779  e = Ci[p];
780  if (w[e] >= mark) {
781  w[e] -= nvi; // decrement |Le\Lk|
782  }
783  else if (w[e] != 0) { // ensure e is a live element
784  w[e] = degree[e] + wnvi; // 1st time e seen in scan 1
785  }
786  } // endfor: scanned Ei
787 
788  } // endfor: scan 1 finished
789 
790  // Scan 2: Degree update
791  for (pk = pk1; pk < pk2; pk++) {
792 
793  i = Ci[pk]; // consider node i in Lk
794  p1 = Cp[i];
795  p2 = p1 + elen [i] - 1;
796  pn = p1;
797 
798  // Scan Ei
799  for (h = 0, d = 0, p = p1 ; p <= p2 ; p++) {
800  e = Ci[p];
801  // If e is an unabsorbed element, then ...
802  if (w[e] != 0) {
803  dext = w[e] - mark; // dext = |Le\Lk|
804  if (dext > 0) {
805  d += dext; // sum up the set differences
806  Ci[pn++] = e; // keep e in Ei
807  h += e; // compute the hash of node i
808  }
809  else {
810  Cp[e] = CS_FLIP(k); // aggressive absorb. e->k
811  w[e] = 0; // e is a dead element
812  }
813  }
814  } // endfor: scanned Ei
815 
816  // elen[i] = |Ei|
817  elen[i] = pn - p1 + 1;
818  p3 = pn ;
819  p4 = p1 + len[i];
820 
821  // Prune edges in Ai
822  for (p = p2 + 1 ; p < p4 ; p++) {
823  j = Ci[p];
824  if ((nvj = nv[j]) <= 0) { // node j dead or in Lk
825  continue;
826  }
827  d += nvj; // degree(i) += |j|
828  Ci[pn++] = j; // place j in node list of i
829  h += j; // compute hash for node i
830  }
831 
832  // Check for mass elimination
833  if (d == 0) {
834  Cp[i] = CS_FLIP(k); // absorb i into k
835  nvi = -nv[i];
836  dk -= nvi; // |Lk| -= |i|
837  nvk += nvi; // |k| += nv[i]
838  nel += nvi;
839  nv[i] = 0;
840  elen[i] = -1; // node i is dead
841  }
842  else {
843  degree[i] = CS_MIN(degree[i], d); // update degree(i)
844  Ci[pn] = Ci[p3]; // move first node to end
845  Ci[p3] = Ci[p1]; // move 1st el. to end of Ei
846  Ci[p1] = k; // add k as 1st element in of Ei
847  len[i] = pn - p1 + 1; // new len of adj. list of node i
848  h %= n; // finalize hash of i
849  next[i] = hhead[h]; // place i in hash bucket
850  hhead[h] = i;
851  last[i] = h; // save hash of i in last[i]
852  }
853 
854  } // endfor: scan 2 is done
855 
856  // finalize |Lk|
857  degree[k] = dk;
858  lemax = CS_MAX(lemax, dk);
859 
860  // Clear w
861  mark = cs_wclear(mark+lemax, lemax, w, n);
862 
863  // Supernode detection
864  for (pk = pk1 ; pk < pk2 ; pk++) {
865 
866  i = Ci[pk];
867  if (nv[i] >= 0) {
868  continue; // skip if i is dead
869  }
870  h = last[i]; // scan hash bucket of node i
871  i = hhead[h];
872  hhead[h] = -1; // hash bucket will be empty
873 
874  // ...
875  for ( ; i != -1 && next[i] != -1 ; i = next[i], mark++) {
876  ln = len[i];
877  eln = elen[i];
878  for (p = Cp[i]+1 ; p <= Cp[i] + ln-1 ; p++) {
879  w[Ci[p]] = mark;
880  }
881  jlast = i;
882 
883  // Compare i with all j
884  for (j = next[i]; j != -1; ) {
885  ok = (len[j] == ln) && (elen[j] == eln);
886  for (p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++) {
887  if (w[Ci[p]] != mark) { // compare i and j
888  ok = 0;
889  }
890  }
891 
892  // i and j are identical
893  if (ok) {
894  Cp[j] = CS_FLIP(i); // absorb j into i
895  nv[i] += nv[j];
896  nv[j] = 0;
897  elen[j] = -1; // node j is dead
898  j = next[j]; // delete j from hash bucket
899  next[jlast] = j;
900  }
901  // j and i are different
902  else {
903  jlast = j;
904  j = next[j];
905  }
906  } // endfor: comparison of i with all j
907 
908  } // endfor: ...
909 
910  } // endfor: supernode detection
911 
912  // Finalize new element Lk
913  for (p = pk1, pk = pk1; pk < pk2; pk++) {
914  i = Ci[pk];
915  if ((nvi = -nv[i]) <= 0) { // skip if i is dead
916  continue;
917  }
918  nv[i] = nvi; // restore nv[i]
919  d = degree[i] + dk - nvi; // compute external degree(i)
920  d = CS_MIN(d, n - nel - nvi);
921  if (head[d] != -1) {
922  last[head[d]] = i;
923  }
924  next[i] = head[d]; // put i back in degree list
925  last[i] = -1;
926  head[d] = i;
927  mindeg = CS_MIN(mindeg, d); // find new minimum degree
928  degree[i] = d;
929  Ci[p++] = i; // place i in Lk
930  }
931 
932  // Number of nodes absorbed into k
933  nv[k] = nvk;
934 
935  // ...
936  if ((len [k] = p-pk1) == 0) { // length of adj list of element k
937  Cp[k] = -1; // k is a root of the tree
938  w[k] = 0; // k is now a dead element
939  }
940  if (elenk != 0) {
941  cnz = p; // free unused space in Lk
942  }
943 
944  } // endwhile
945 
946 
947  // ======================================================================
948  // Step 5: Postordering
949  // ======================================================================
950 
951  // Debug
952  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
953  std::cout << " Step 5: Postordering" << std::endl;
954  #endif
955 
956  // Fix assembly tree
957  for (i = 0; i < n; i++) {
958  Cp[i] = CS_FLIP(Cp[i]);
959  }
960  for (j = 0 ; j <= n ; j++) {
961  head[j] = -1;
962  }
963 
964  // Place unordered nodes in lists
965  for (j = n ; j >= 0 ; j--) {
966  if (nv[j] > 0) {
967  continue; // skip if j is an element
968  }
969  next[j] = head [Cp[j]]; // place j in list of its parent
970  head[Cp[j]] = j;
971  }
972 
973  // Place elements in lists
974  for (e = n ; e >= 0 ; e--) {
975  if (nv[e] <= 0) {
976  continue; // skip unless e is an element
977  }
978  if (Cp[e] != -1) {
979  next[e] = head[Cp[e]]; // place e in list of its parent
980  head[Cp[e]] = e;
981  }
982  }
983 
984 
985  // ======================================================================
986  // Step 6: Postorder the assembly tree
987  // ======================================================================
988 
989  // Debug
990  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
991  std::cout << " Step 6: Postorder the assembly tree" << std::endl;
992  #endif
993 
994  // Postorder the assembly tree
995  for (k = 0, i = 0 ; i <= n ; i++) {
996  if (Cp [i] == -1) {
997  k = cs_tdfs(i, k, head, next, P, w);
998  }
999  }
1000 
1001  // Free workspace
1002  delete [] wrk_int;
1003 
1004  // Debug
1005  #if defined(G_DEBUG_SPARSE_SYM_AMD_ORDERING)
1006  std::cout << " Permutation:" << std::endl;
1007  for (i = 0; i < n; ++i) {
1008  std::cout << " " << P[i];
1009  }
1010  std::cout << endl;
1011  std::cout << "GSparseSymbolic::cs_amd finished" << std::endl;
1012  #endif
1013 
1014  // Return result
1015  return P;
1016 }
1017 
1018 
1019 /***************************************************************************
1020  * @brief cs_counts
1021  *
1022  * @param[in] A Spare matrix.
1023  * @param[in] parent Parent array.
1024  * @param[in] post Post ordering array.
1025  * @param[in] ata 0: count LL'=A, 1: count LL'=A'A
1026  * @return Integer array of n+1 elements
1027  *
1028  * Column counts of LL'=A or LL'=A'A, given parent & post ordering
1029  ***************************************************************************/
1031  const int* parent,
1032  const int* post,
1033  int ata)
1034 {
1035  // Declare loop variables
1036  int i, j, k, J, p, q, jleaf;
1037 
1038  // Return NULL pointer if input arrays are NULL or sparse matrix is
1039  // empty
1040  if (!parent || !post) {
1041  return NULL;
1042  }
1043 
1044  // Assign input matrix attributes
1045  int m = A->m_rows;
1046  int n = A->m_cols;
1047 
1048  // Allocate result
1049  int* colcount = new int[n];
1050  int* delta = colcount;
1051 
1052  // Allocate workspace
1053  int wrk_size = 4*n + (ata ? (n+m+1) : 0);
1054  int* wrk_int = new int[wrk_size];
1055 
1056  // Get (logical) transpose of A: AT = A'
1057  GMatrixSparse AT = cs_transpose(*A, 0);
1058 
1059  // Set-up pointers to workspace
1060  int* ancestor = wrk_int;
1061  int* maxfirst = wrk_int + n;
1062  int* prevleaf = wrk_int + 2*n;
1063  int* first = wrk_int + 3*n;
1064 
1065  // Clear workspace
1066  for (k = 0; k < wrk_size; ++k) {
1067  wrk_int[k] = -1;
1068  }
1069 
1070  // Find first j
1071  for (k = 0; k < n; k++) {
1072  j = post[k];
1073  delta[j] = (first[j] == -1) ? 1 : 0; // delta[j]=1 if j is a leaf
1074  for ( ; j != -1 && first[j] == -1; j = parent[j]) {
1075  first[j] = k;
1076  }
1077  }
1078 
1079  // Assign output matrix attributes
1080  int* ATp = AT.m_colstart;
1081  int* ATi = AT.m_rowinx;
1082 
1083  // Initialise. Note that init_ata wires the working array into the *head
1084  // and *next pointers, so both arrays are just aliases of the wrk_int
1085  // array.
1086  int* head = NULL;
1087  int* next = NULL;
1088  if (ata) {
1089  init_ata(&AT, post, wrk_int, &head, &next);
1090  }
1091 
1092  // Each node in its own set
1093  for (i = 0; i < n; i++) {
1094  ancestor[i] = i;
1095  }
1096 
1097  // Loop over all columns of matrix A (or rows of matrix AT)
1098  for (k = 0; k < n; k++) {
1099  j = post[k]; // j is the kth node in postordered etree
1100  if (parent[j] != -1) {
1101  delta[parent[j]]--; // j is not a root
1102  }
1103  for (J = HEAD(k,j); J != -1; J = NEXT(J)) { // J=j for LL'=A case
1104 
1105  // Loop over all non-zero elements in column J
1106  for (p = ATp[J]; p < ATp[J+1]; p++) {
1107  i = ATi[p];
1108  q = cs_leaf(i, j, first, maxfirst, prevleaf, ancestor, &jleaf);
1109  if (jleaf >= 1) {
1110  delta[j]++; // A(i,j) is in skeleton
1111  }
1112  if (jleaf == 2) {
1113  delta[q]--; // account for overlap in q
1114  }
1115  }
1116  }
1117  if (parent[j] != -1) {
1118  ancestor[j] = parent[j];
1119  }
1120  }
1121 
1122  // Sum up delta's of each child
1123  for (j = 0; j < n; j++) {
1124  if (parent[j] != -1) {
1125  colcount[parent[j]] += colcount[j];
1126  }
1127  }
1128 
1129  // Free workspace
1130  delete [] wrk_int;
1131 
1132  // Return result
1133  return colcount;
1134 }
1135 
1136 
1137 /***************************************************************************
1138  * @brief cs_etree
1139  *
1140  * @param[in] A Sparse matrix
1141  * @param[in] ata Flag.
1142  * @return Evaluation tree of matrix @p A.
1143  *
1144  * Compute the etree of A using triu(A), or A'A without forming A'A
1145  ***************************************************************************/
1147 {
1148  // Declare loop variables
1149  int i, k, p;
1150 
1151  // Assign matrix attributes
1152  int m = A->m_rows;
1153  int n = A->m_cols;
1154 
1155  // Allocate result array
1156  int* parent = new int[n];
1157 
1158  // Allocate workspace
1159  int wrk_size = n + (ata ? m : 0);
1160  int* wrk_int = new int[wrk_size];
1161 
1162  // Set-up pointers to workspace. If 'ata=1' then we use also a prev array
1163  int* ancestor = wrk_int;
1164  int* prev = wrk_int + n;
1165 
1166  // If 'prev' is requested, initialise array with -1
1167  if (ata) {
1168  for (i = 0; i < m ; ++i) {
1169  prev[i] = -1;
1170  }
1171  }
1172 
1173  // Loop over all nodes (columns)
1174  for (k = 0; k < n; ++k) {
1175 
1176  // Initialise node k to having no parent and no ancestor
1177  parent[k] = -1; // node k has no parent yet
1178  ancestor[k] = -1; // nor does k have an ancestor
1179 
1180  // Loop over all elements in column k
1181  for (p = A->m_colstart[k]; p < A->m_colstart[k+1]; ++p) {
1182 
1183  // Get element index
1184  i = ata ? (prev[A->m_rowinx[p]]) : (A->m_rowinx[p]);
1185 
1186  // Traverse from i to k
1187  int inext;
1188  for ( ; i != -1 && i < k; i = inext) {
1189  inext = ancestor[i]; // inext = ancestor of i
1190  ancestor[i] = k; // path compression
1191  if (inext == -1) {
1192  parent[i] = k; // no ancestor, parent is k
1193  }
1194  }
1195  if (ata) {
1196  prev[A->m_rowinx[p]] = k;
1197  }
1198 
1199  } // endfor: looped over all elements in column k
1200 
1201  } // endfor: looped over all columns
1202 
1203  // Free workspace
1204  delete [] wrk_int;
1205 
1206  // Return result
1207  return parent;
1208 }
1209 
1210 
1211 /***********************************************************************//**
1212  * @brief cs_fkeep
1213  *
1214  * @param[in] A Sparse matrix.
1215  * @param[in] fkeep Screening function.
1216  * @param[in] other Other function.
1217  * @return Number of non-zero elements (-1 if not okay)
1218  *
1219  * Drop entries for which fkeep(A(i,j)) is false. Return the number of
1220  * non-zero elements if ok, otherwise return -1.
1221  ***************************************************************************/
1223  int(*fkeep)(int, int, double, void*),
1224  void* other)
1225 {
1226  // Return error if some of the input pointers is invalid
1227  if (!A || !fkeep) {
1228  return (-1);
1229  }
1230 
1231  // Initialise counter for non-zero elements
1232  int nz = 0;
1233 
1234  // Assign matrix attributes
1235  int n = A->m_cols;
1236  int* Ap = A->m_colstart;
1237  int* Ai = A->m_rowinx;
1238  double* Ax = A->m_data;
1239 
1240  // Operate only if we have some data
1241  if (Ap != NULL && Ai != NULL && Ax != NULL) {
1242 
1243  // Loop over all columns
1244  for (int j = 0; j < n; j++) {
1245  int p = Ap[j]; // get current location of col j
1246  Ap[j] = nz; // record new location of col j
1247  for ( ; p < Ap[j+1] ; p++) {
1248  if (fkeep(Ai[p], j, Ax ? Ax[p] : 1, other)) {
1249  if (Ax) {
1250  Ax[nz] = Ax[p]; // keep A(i,j)
1251  }
1252  Ai[nz++] = Ai[p];
1253  }
1254  }
1255  }
1256 
1257  // Finalise A
1258  Ap[n] = nz;
1259  A->m_elements = nz;
1260 
1261  } // endif: we had data
1262 
1263  // Remove extra space from A
1264  //cs_sprealloc(A, 0);
1265  A->free_elements(nz, (A->m_elements-nz));
1266 
1267  // Return number of non-zero elements
1268  return nz;
1269 }
1270 
1271 
1272 /***************************************************************************
1273  * @brief cs_leaf
1274  *
1275  * @param[in] i Row to consider
1276  * @param[in] j Column to consider
1277  * @param[in] first ...
1278  * @param[in] maxfirst ...
1279  * @param[in] prevleaf ...
1280  * @param[in] ancestor ...
1281  * @param[out] jleaf ...
1282  * @return q
1283  *
1284  * Consider A(i,j), node j in ith row subtree and return lca(jprev,j)
1285  ***************************************************************************/
1287  int j,
1288  const int *first,
1289  int *maxfirst,
1290  int *prevleaf,
1291  int *ancestor,
1292  int *jleaf)
1293 {
1294  // If one of the input pointers is invalid then return -1
1295  if (!first || !maxfirst || !prevleaf || !ancestor || !jleaf) {
1296  return (-1);
1297  }
1298 
1299  // Declare loop variables
1300  int q, s, sparent;
1301 
1302  // ...
1303  *jleaf = 0 ;
1304 
1305  // If j is not a leaf then return -1
1306  if (i <= j || first[j] <= maxfirst[i]) {
1307  return (-1);
1308  }
1309 
1310  // Update max first[j] seen so far
1311  maxfirst[i] = first[j];
1312 
1313  // jprev = previous leaf of ith subtree
1314  int jprev = prevleaf[i];
1315  prevleaf[i] = j;
1316 
1317  // j is first or subsequent leaf
1318  *jleaf = (jprev == -1) ? 1: 2;
1319 
1320  // If j is 1st leaf then q is the root of ith subtree
1321  if (*jleaf == 1) {
1322  return (i);
1323  }
1324 
1325  // ...
1326  for (q = jprev; q != ancestor[q]; q = ancestor[q]) ;
1327 
1328  // ...
1329  for (s = jprev; s != q; s = sparent) {
1330  sparent = ancestor[s]; // path compression
1331  ancestor[s] = q;
1332  }
1333 
1334  // Return least common ancester (jprev,j)
1335  return (q);
1336 }
1337 
1338 
1339 /***************************************************************************
1340  * @brief Inverts the permutation p[0..n-1]
1341  *
1342  * @param[in] p Integer array of n elements (NULL denotes ID)
1343  * @param[in] n Number of elements
1344  * @return Integer array of n elements pinv[0..n-1] = p[n-1..0]
1345  *
1346  * Inverts the permutation p[0..n-1]
1347  ***************************************************************************/
1348 int* GSparseSymbolic::cs_pinv(int const *p, int n)
1349 {
1350  // Return NULL pointer if input pointer is NULL or the number of
1351  // elements is zero. This denotes identity.
1352  if (!p || n < 1) {
1353  return NULL;
1354  }
1355 
1356  // Allocate result array
1357  int* pinv = new int[n];
1358 
1359  // Invert the permutation
1360  for (int k = 0; k < n; ++k) {
1361  pinv[p[k]] = k;
1362  }
1363 
1364  // Return result
1365  return pinv;
1366 }
1367 
1368 
1369 /***************************************************************************
1370  * @brief Post order a forest
1371  *
1372  * @param[in] parent Integer array of n elements
1373  * @param[in] n Number of elements
1374  * @return Integer array of n elements
1375  *
1376  * Post order a forest
1377  ***************************************************************************/
1378 int* GSparseSymbolic::cs_post(const int* parent, int n)
1379 {
1380  // Return NULL pointer if input pointer is NULL or number of elements
1381  // is zero
1382  if (!parent || n < 1) {
1383  return (NULL);
1384  }
1385 
1386  // Allocate result array
1387  int* post = new int[n];
1388 
1389  // Allocate workspace
1390  int wrk_size = 3 * n;
1391  int* wrk_int = new int[wrk_size];
1392 
1393  // Set-up pointers to workspace
1394  int* head = wrk_int;
1395  int* next = wrk_int + n;
1396  int* stack = wrk_int + 2*n; // Working array for 'cs_tdfs' function
1397 
1398  // Declare and initialise loop variables
1399  int j, k = 0;
1400 
1401  // Empty linked list
1402  for (j = 0 ; j < n ; j++) {
1403  head[j] = -1;
1404  }
1405 
1406  // Traverse nodes in reverse order to build a linked list 'head'
1407  for (j = n-1; j >= 0; j--) {
1408 
1409  // Skip if j is a root
1410  if (parent[j] == -1) {
1411  continue;
1412  }
1413 
1414  // Add j to list of its parent
1415  next[j] = head[parent[j]];
1416  head[parent[j]] = j;
1417 
1418  }
1419 
1420  // Traverse nodes
1421  for (j = 0 ; j < n ; j++) {
1422 
1423  // Skip if j is not a root
1424  if (parent[j] != -1) {
1425  continue;
1426  }
1427 
1428  // Depth-first search and postorder of a tree rooted at node j
1429  k = cs_tdfs(j, k, head, next, post, stack);
1430 
1431  }
1432 
1433  // Free workspace
1434  delete [] wrk_int;
1435 
1436  // Return result
1437  return post;
1438 }
1439 
1440 
1441 /***************************************************************************
1442  * @brief cs_tdfs
1443  *
1444  * @param[in] j Node of tree root
1445  * @param[in] k ...
1446  * @param[in,out] head Linked list
1447  * @param[in] next Linked list information
1448  * @param[in,out] post Array
1449  * @param[in,OUT] stack Working array
1450  * @return Update of k
1451  *
1452  * Depth-first search and postorder of a tree rooted at node j
1453  ***************************************************************************/
1455  int k,
1456  int* head,
1457  const int* next,
1458  int* post,
1459  int* stack)
1460 {
1461  // Return -1 if one of the array pointers is NULL
1462  if (!head || !next || !post || !stack) {
1463  return (-1);
1464  }
1465 
1466  // Decrare
1467  int i, p;
1468 
1469  // Place j on the first position of the stack
1470  stack[0] = j;
1471 
1472  // Loop while stack is not empty
1473  int top = 0;
1474  while (top >= 0) {
1475  p = stack[top]; // p = top of stack
1476  i = head[p]; // i = youngest child of p
1477  if (i == -1) {
1478  top--; // p has no unordered children left
1479  post[k++] = p; // node p is the kth postordered node
1480  }
1481  else {
1482  head[p] = next[i]; // remove i from children of p
1483  stack[++top] = i; // start depth-first search on child node i
1484  }
1485  }
1486 
1487  // Return result
1488  return k;
1489 }
1490 
1491 
1492 /*==========================================================================
1493  = =
1494  = GSparseSymbolic static member functions =
1495  = =
1496  ==========================================================================*/
1497 
1498 /***************************************************************************
1499  * @brief Initialise A'A
1500  *
1501  * @param[in] AT Sparse matrix.
1502  * @param[in] post ...
1503  * @param[in] wrk_int Workspace.
1504  * @param[out] head Pointer to linked list.
1505  * @param[out] next Pointer to linked list information.
1506  *
1507  * Initialise A'A.
1508  *
1509  * Note that the returned pointers @p head and @p info point into the
1510  * @p wrk_int workspace array.
1511  ***************************************************************************/
1513  const int* post,
1514  int* wrk_int,
1515  int** head,
1516  int** next)
1517 {
1518  // Declare loop variables
1519  int i, k, p;
1520 
1521  // Assign matrix attributes
1522  int m = AT->m_cols;
1523  int n = AT->m_rows;
1524  int* ATp = AT->m_colstart;
1525  int* ATi = AT->m_rowinx;
1526 
1527  // Set-up pointers to workspace
1528  *head = wrk_int + 4*n;
1529  *next = wrk_int + 5*n+1;
1530 
1531  // Invert post
1532  for (k = 0 ; k < n ; k++) {
1533  wrk_int[post[k]] = k;
1534  }
1535 
1536  // Loop over columns of AT matrix
1537  for (i = 0 ; i < m ; i++) {
1538 
1539  // Loop over elements in column i of AT matrix
1540  for (k = n, p = ATp[i]; p < ATp[i+1]; p++) {
1541  k = CS_MIN(k,wrk_int[ATi[p]]);
1542  }
1543 
1544  // Place row i in linked list k
1545  (*next)[i] = (*head)[k];
1546  (*head)[k] = i;
1547 
1548  } // endfor: looped over columns of matrix
1549 
1550  // Return
1551  return;
1552 }
1553 
1554 
1555 /***************************************************************************
1556  * @brief cs_diag
1557  *
1558  * @param[in] i Row index.
1559  * @param[in] j Column index.
1560  * @param[in] aij Array element (not used).
1561  * @param[in] other Not used
1562  * @return 0: off-diagonal, 1: diagonal element
1563  *
1564  * Keep off-diagonal entries and drop diagonal entries.
1565  ***************************************************************************/
1566 int GSparseSymbolic::cs_diag(int i, int j, double aij, void* other)
1567 {
1568  return (i != j);
1569 }
1570 
1571 
1572 /***************************************************************************
1573  * @brief cs_wclear
1574  *
1575  * @param[in] mark ...
1576  * @param[in] lemax ...
1577  * @param[in,out] w Working array of size @p n.
1578  * @param[in] n Size of the working array
1579  * @return Update of @p mark
1580  *
1581  * Clear w array.
1582  ***************************************************************************/
1583 int GSparseSymbolic::cs_wclear(int mark, int lemax, int* w, int n)
1584 {
1585  // ...
1586  if (mark < 2 || (mark + lemax < 0)) {
1587  for (int k = 0; k < n; k++) {
1588  if (w[k] != 0) {
1589  w[k] = 1;
1590  }
1591  }
1592  mark = 2;
1593  }
1594 
1595  // Return mark. At this point, w[0..n-1] < mark holds
1596  return mark;
1597 }
1598 
1599 
1600 /*==========================================================================
1601  = =
1602  = GSparseSymbolic friends =
1603  = =
1604  ==========================================================================*/
1605 
1606 /***************************************************************************
1607  * Output operator *
1608  ***************************************************************************/
1609 std::ostream& operator<< (std::ostream& os, const GSparseSymbolic& s)
1610 {
1611  // Put header is stream
1612  os << "=== GSparseSymbolic ===";
1613 
1614  // Show inverse permutation if it exists
1615  if (s.m_pinv != NULL) {
1616  os << std::endl << " Inverse permutation (of ";
1617  os << s.m_n_pinv << " elements)" << std::endl;
1618  for (int i = 0; i < s.m_n_pinv; ++i) {
1619  os << " " << s.m_pinv[i];
1620  }
1621  }
1622 
1623  // Show fill reducing column permutation if it exists
1624  if (s.m_q != NULL) {
1625  os << std::endl << " Fill reducing column permutation (";
1626  os << s.m_n_cp << " elements)" << std::endl;
1627  for (int i = 0; i < s.m_n_cp; ++i) {
1628  os << " " << s.m_q[i];
1629  }
1630  }
1631 
1632  // Show elimination tree for Cholesky and QR if it exists
1633  if (s.m_parent != NULL) {
1634  os << std::endl << " Elimination tree (";
1635  os << s.m_n_parent << " elements)" << std::endl;
1636  for (int i = 0; i < s.m_n_parent; ++i) {
1637  os << " " << s.m_parent[i];
1638  }
1639  }
1640 
1641  // Show column pointers for Cholesky row counts for QR if it exists
1642  if (s.m_cp != NULL) {
1643  os << std::endl << " Cholesky column pointers/QR row counts (";
1644  os << s.m_n_cp << " elements)" << std::endl;
1645  for (int i = 0; i < s.m_n_cp; ++i) {
1646  os << " " << s.m_cp[i];
1647  }
1648  }
1649 
1650  // Show leftmost[i] = min(find(A(i,:))) if it exists
1651  if (s.m_leftmost != NULL) {
1652  os << std::endl << " leftmost[i] = min(find(A(i,:))) (";
1653  os << s.m_n_leftmost << " elements)" << std::endl;
1654  for (int i = 0; i < s.m_n_leftmost; ++i) {
1655  os << " " << s.m_leftmost[i];
1656  }
1657  }
1658 
1659  // Return output stream
1660  return os;
1661 }
static void init_ata(const GMatrixSparse *AT, const int *post, int *wrk_int, int **head, int **next)
GMatrixSparse cs_symperm(const GMatrixSparse &matrix, const int *pinv)
cs_symperm
void cholesky_symbolic_analysis(int order, const GMatrixSparse &m)
int * cs_post(const int *parent, int n)
int * m_parent
elimination tree for Cholesky and QR
int * m_cp
column pointers for Cholesky, row counts for QR
double m_lnz
entries in L for LU or Cholesky; in V for QR
Sparse matrix class interface definition.
int * m_leftmost
leftmost[i] = min(find(A(i,:))), for QR
#define CS_MIN(a, b)
#define CS_MAX(a, b)
GMatrixSparse cs_transpose(const GMatrixSparse &matrix, int values)
int m_rows
Number of rows.
Sparse matrix symbolic analysis class definition.
int cs_leaf(int i, int j, const int *first, int *maxfirst, int *prevleaf, int *ancestor, int *jleaf)
int cs_tdfs(int j, int k, int *head, const int *next, int *post, int *stack)
Gammalib tools definition.
int m_n_cp
Number of elements in m_cp.
#define G_CHOLESKY
int m_n_parent
Number of elements in m_parent.
static int cs_wclear(int mark, int lemax, int *w, int n)
virtual ~GSparseSymbolic(void)
Destructor.
Sparse matrix symbolic analysis class.
GSparseSymbolic(void)
Void constructor.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
int m_n_leftmost
Number of elements in m_leftmost.
double * m_data
Matrix data.
int * m_pinv
Inverse row permutation for QR, fill reduce permutation for Cholesky.
GSparseSymbolic & operator=(const GSparseSymbolic &s)
void free_elements(const int &start, const int &num)
Free memory for obsolete matrix elements.
int m_cols
Number of columns.
std::ostream & operator<<(std::ostream &os, const GBase &base)
Output operator.
Definition: GBase.cpp:57
double cs_cumsum(int *p, int *c, int n)
cs_cumsum
int * m_q
Fill-reducing column permutation for LU and QR.
int * cs_amd(int order, const GMatrixSparse *A)
int * m_rowinx
Row-indices of all elements.
int * cs_counts(const GMatrixSparse *A, const int *parent, const int *post, int ata)
int m_n_q
Number of elements in m_q.
int * cs_etree(const GMatrixSparse *A, int ata)
#define NEXT(J)
int * cs_pinv(int const *p, int n)
int m_elements
Number of elements stored in matrix.
double m_unz
entries in U for LU; in R for QR
Sparse matrix class definition.
Exception handler interface definition.
int m_n_pinv
Number of elements in m_pinv.
#define HEAD(k, j)
int cs_fkeep(GMatrixSparse *A, int(*fkeep)(int, int, double, void *), void *other)
cs_fkeep
int m_m2
of rows for QR, after adding fictitious rows
int * m_colstart
Column start indices (m_cols+1)
static int cs_diag(int i, int j, double aij, void *other)
#define CS_FLIP(i)
#define G_CS_AMD
void alloc_elements(int start, const int &num)
Allocate memory for new matrix elements.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489