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