forked from RcppCore/RcppEigen
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathRcppEigenCholmod.h
More file actions
1090 lines (922 loc) · 49.2 KB
/
Copy pathRcppEigenCholmod.h
File metadata and controls
1090 lines (922 loc) · 49.2 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
//
// cholmod.h: selected headers from Tim Davis's CHOLMOD package
//
// Copyright (C) 2011 Douglas Bates, Martin Maechler, Dirk Eddelbuettel and Romain Francois
//
// This file is part of RcppEigen.
//
// RcppEigen is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 2 of the License, or
// (at your option) any later version.
//
// RcppEigen is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with RcppEigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef RcppEigen_CHOLMOD_H
#define RcppEigen_CHOLMOD_H
/* Original of this is Matrix/inst/include/cholmod.h --- sync with it (e.g for lme4) !
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <stddef.h>
#include <limits.h>
// Rather use C99 -- which we require in R anyway
#include <inttypes.h>
#ifdef __cplusplus
extern "C" {
#endif
// from Matrix/src/SuiteSparse_config/SuiteSparse_config.h - line 51 :
#ifndef SuiteSparse_long
/* #ifdef _WIN64 */
/* #define SuiteSparse_long __int64 */
/* #define SuiteSparse_long_max _I64_MAX */
/* #define SuiteSparse_long_idd "I64d" */
/* #else */
/* #define SuiteSparse_long long */
/* #define SuiteSparse_long_max LONG_MAX */
/* #define SuiteSparse_long_idd "ld" */
/* #endif */
#define SuiteSparse_long int64_t
// typically long (but on WIN64)
#define SuiteSparse_ulong uint64_t
// only needed for ../COLAMD/Source/colamd.c (original has 'unsigned Int' which fails!!)
#define SuiteSparse_long_max 9223372036854775801
// typically LONG_MAX (but ..)
#define SuiteSparse_long_idd PRId64
// typically "ld"
#define SuiteSparse_long_id "%" SuiteSparse_long_idd
#endif
/* For backward compatibility with prior versions of SuiteSparse. The UF_*
* macros are deprecated and will be removed in a future version. */
#ifndef UF_long
#define UF_long SuiteSparse_long
#define UF_long_max SuiteSparse_long_max
#define UF_long_idd SuiteSparse_long_idd
#define UF_long_id SuiteSparse_long_id
#endif
#define CHOLMOD_HAS_VERSION_FUNCTION
#define CHOLMOD_DATE "April 25, 2013"
#define CHOLMOD_VER_CODE(main,sub) ((main) * 1000 + (sub))
#define CHOLMOD_MAIN_VERSION 2
#define CHOLMOD_SUB_VERSION 1
#define CHOLMOD_SUBSUB_VERSION 2
#define CHOLMOD_VERSION \
CHOLMOD_VER_CODE(CHOLMOD_MAIN_VERSION,CHOLMOD_SUB_VERSION)
// from ../../src/CHOLMOD/Include/cholmod_core.h - line 275 :
/* Each CHOLMOD object has its own type code. */
#define CHOLMOD_COMMON 0
#define CHOLMOD_SPARSE 1
#define CHOLMOD_FACTOR 2
#define CHOLMOD_DENSE 3
#define CHOLMOD_TRIPLET 4
/* ========================================================================== */
/* === CHOLMOD Common ======================================================= */
/* ========================================================================== */
/* itype defines the types of integer used: */
#define CHOLMOD_INT 0 /* all integer arrays are int */
#define CHOLMOD_INTLONG 1 /* most are int, some are SuiteSparse_long */
#define CHOLMOD_LONG 2 /* all integer arrays are SuiteSparse_long */
/* The itype of all parameters for all CHOLMOD routines must match.
* FUTURE WORK: CHOLMOD_INTLONG is not yet supported.
*/
/* dtype defines what the numerical type is (double or float): */
#define CHOLMOD_DOUBLE 0 /* all numerical values are double */
#define CHOLMOD_SINGLE 1 /* all numerical values are float */
/* The dtype of all parameters for all CHOLMOD routines must match.
*
* Scalar floating-point values are always passed as double arrays of size 2
* (for the real and imaginary parts). They are typecast to float as needed.
* FUTURE WORK: the float case is not supported yet.
*/
/* xtype defines the kind of numerical values used: */
#define CHOLMOD_PATTERN 0 /* pattern only, no numerical values */
#define CHOLMOD_REAL 1 /* a real matrix */
#define CHOLMOD_COMPLEX 2 /* a complex matrix (ANSI C99 compatible) */
#define CHOLMOD_ZOMPLEX 3 /* a complex matrix (MATLAB compatible) */
/* Definitions for cholmod_common: */
#define CHOLMOD_MAXMETHODS 9 /* maximum number of different methods that */
/* cholmod_analyze can try. Must be >= 9. */
/* Common->status values. zero means success, negative means a fatal error,
* positive is a warning. */
#define CHOLMOD_OK 0 /* success */
#define CHOLMOD_NOT_INSTALLED (-1) /* failure: method not installed */
#define CHOLMOD_OUT_OF_MEMORY (-2) /* failure: out of memory */
#define CHOLMOD_TOO_LARGE (-3) /* failure: integer overflow occured */
#define CHOLMOD_INVALID (-4) /* failure: invalid input */
#define CHOLMOD_GPU_PROBLEM (-5) /* failure: GPU fatal error */
#define CHOLMOD_NOT_POSDEF (1) /* warning: matrix not pos. def. */
#define CHOLMOD_DSMALL (2) /* warning: D for LDL' or diag(L) or */
/* LL' has tiny absolute value */
/* ordering method (also used for L->ordering) */
#define CHOLMOD_NATURAL 0 /* use natural ordering */
#define CHOLMOD_GIVEN 1 /* use given permutation */
#define CHOLMOD_AMD 2 /* use minimum degree (AMD) */
#define CHOLMOD_METIS 3 /* use METIS' nested dissection */
#define CHOLMOD_NESDIS 4 /* use CHOLMOD's version of nested dissection:*/
/* node bisector applied recursively, followed
* by constrained minimum degree (CSYMAMD or
* CCOLAMD) */
#define CHOLMOD_COLAMD 5 /* use AMD for A, COLAMD for A*A' */
/* POSTORDERED is not a method, but a result of natural ordering followed by a
* weighted postorder. It is used for L->ordering, not method [ ].ordering. */
#define CHOLMOD_POSTORDERED 6 /* natural ordering, postordered. */
/* supernodal strategy (for Common->supernodal) */
#define CHOLMOD_SIMPLICIAL 0 /* always do simplicial */
#define CHOLMOD_AUTO 1 /* select simpl/super depending on matrix */
#define CHOLMOD_SUPERNODAL 2 /* always do supernodal */
typedef struct cholmod_common_struct
{
/* ---------------------------------------------------------------------- */
/* parameters for symbolic/numeric factorization and update/downdate */
/* ---------------------------------------------------------------------- */
double dbound ; /* Smallest absolute value of diagonal entries of D
* for LDL' factorization and update/downdate/rowadd/
* rowdel, or the diagonal of L for an LL' factorization.
* Entries in the range 0 to dbound are replaced with dbound.
* Entries in the range -dbound to 0 are replaced with -dbound. No
* changes are made to the diagonal if dbound <= 0. Default: zero */
double grow0 ; /* For a simplicial factorization, L->i and L->x can
* grow if necessary. grow0 is the factor by which
* it grows. For the initial space, L is of size MAX (1,grow0) times
* the required space. If L runs out of space, the new size of L is
* MAX(1.2,grow0) times the new required space. If you do not plan on
* modifying the LDL' factorization in the Modify module, set grow0 to
* zero (or set grow2 to 0, see below). Default: 1.2 */
double grow1 ;
size_t grow2 ; /* For a simplicial factorization, each column j of L
* is initialized with space equal to
* grow1*L->ColCount[j] + grow2. If grow0 < 1, grow1 < 1, or grow2 == 0,
* then the space allocated is exactly equal to L->ColCount[j]. If the
* column j runs out of space, it increases to grow1*need + grow2 in
* size, where need is the total # of nonzeros in that column. If you do
* not plan on modifying the factorization in the Modify module, set
* grow2 to zero. Default: grow1 = 1.2, grow2 = 5. */
size_t maxrank ; /* rank of maximum update/downdate. Valid values:
* 2, 4, or 8. A value < 2 is set to 2, and a
* value > 8 is set to 8. It is then rounded up to the next highest
* power of 2, if not already a power of 2. Workspace (Xwork, below) of
* size nrow-by-maxrank double's is allocated for the update/downdate.
* If an update/downdate of rank-k is requested, with k > maxrank,
* it is done in steps of maxrank. Default: 8, which is fastest.
* Memory usage can be reduced by setting maxrank to 2 or 4.
*/
double supernodal_switch ; /* supernodal vs simplicial factorization */
int supernodal ; /* If Common->supernodal <= CHOLMOD_SIMPLICIAL
* (0) then cholmod_analyze performs a
* simplicial analysis. If >= CHOLMOD_SUPERNODAL (2), then a supernodal
* analysis is performed. If == CHOLMOD_AUTO (1) and
* flop/nnz(L) < Common->supernodal_switch, then a simplicial analysis
* is done. A supernodal analysis done otherwise.
* Default: CHOLMOD_AUTO. Default supernodal_switch = 40 */
int final_asis ; /* If TRUE, then ignore the other final_* parameters
* (except for final_pack).
* The factor is left as-is when done. Default: TRUE.*/
int final_super ; /* If TRUE, leave a factor in supernodal form when
* supernodal factorization is finished. If FALSE,
* then convert to a simplicial factor when done.
* Default: TRUE */
int final_ll ; /* If TRUE, leave factor in LL' form when done.
* Otherwise, leave in LDL' form. Default: FALSE */
int final_pack ; /* If TRUE, pack the columns when done. If TRUE, and
* cholmod_factorize is called with a symbolic L, L is
* allocated with exactly the space required, using L->ColCount. If you
* plan on modifying the factorization, set Common->final_pack to FALSE,
* and each column will be given a little extra slack space for future
* growth in fill-in due to updates. Default: TRUE */
int final_monotonic ; /* If TRUE, ensure columns are monotonic when done.
* Default: TRUE */
int final_resymbol ;/* if cholmod_factorize performed a supernodal
* factorization, final_resymbol is true, and
* final_super is FALSE (convert a simplicial numeric factorization),
* then numerically zero entries that resulted from relaxed supernodal
* amalgamation are removed. This does not remove entries that are zero
* due to exact numeric cancellation, since doing so would break the
* update/downdate rowadd/rowdel routines. Default: FALSE. */
/* supernodal relaxed amalgamation parameters: */
double zrelax [3] ;
size_t nrelax [3] ;
/* Let ns be the total number of columns in two adjacent supernodes.
* Let z be the fraction of zero entries in the two supernodes if they
* are merged (z includes zero entries from prior amalgamations). The
* two supernodes are merged if:
* (ns <= nrelax [0]) || (no new zero entries added) ||
* (ns <= nrelax [1] && z < zrelax [0]) ||
* (ns <= nrelax [2] && z < zrelax [1]) || (z < zrelax [2])
*
* Default parameters result in the following rule:
* (ns <= 4) || (no new zero entries added) ||
* (ns <= 16 && z < 0.8) || (ns <= 48 && z < 0.1) || (z < 0.05)
*/
int prefer_zomplex ; /* X = cholmod_solve (sys, L, B, Common) computes
* x=A\b or solves a related system. If L and B are
* both real, then X is real. Otherwise, X is returned as
* CHOLMOD_COMPLEX if Common->prefer_zomplex is FALSE, or
* CHOLMOD_ZOMPLEX if Common->prefer_zomplex is TRUE. This parameter
* is needed because there is no supernodal zomplex L. Suppose the
* caller wants all complex matrices to be stored in zomplex form
* (MATLAB, for example). A supernodal L is returned in complex form
* if A is zomplex. B can be real, and thus X = cholmod_solve (L,B)
* should return X as zomplex. This cannot be inferred from the input
* arguments L and B. Default: FALSE, since all data types are
* supported in CHOLMOD_COMPLEX form and since this is the native type
* of LAPACK and the BLAS. Note that the MATLAB/cholmod.c mexFunction
* sets this parameter to TRUE, since MATLAB matrices are in
* CHOLMOD_ZOMPLEX form.
*/
int prefer_upper ; /* cholmod_analyze and cholmod_factorize work
* fastest when a symmetric matrix is stored in
* upper triangular form when a fill-reducing ordering is used. In
* MATLAB, this corresponds to how x=A\b works. When the matrix is
* ordered as-is, they work fastest when a symmetric matrix is in lower
* triangular form. In MATLAB, R=chol(A) does the opposite. This
* parameter affects only how cholmod_read returns a symmetric matrix.
* If TRUE (the default case), a symmetric matrix is always returned in
* upper-triangular form (A->stype = 1). */
int quick_return_if_not_posdef ; /* if TRUE, the supernodal numeric
* factorization will return quickly if
* the matrix is not positive definite. Default: FALSE. */
/* ---------------------------------------------------------------------- */
/* printing and error handling options */
/* ---------------------------------------------------------------------- */
int print ; /* print level. Default: 3 */
int precise ; /* if TRUE, print 16 digits. Otherwise print 5 */
int (*print_function) (const char *, ...) ; /* pointer to printf */
int try_catch ; /* if TRUE, then ignore errors; CHOLMOD is in the middle
* of a try/catch block. No error message is printed
* and the Common->error_handler function is not called. */
void (*error_handler) (int status, const char *file,
int line, const char *message) ;
/* Common->error_handler is the user's error handling routine. If not
* NULL, this routine is called if an error occurs in CHOLMOD. status
* can be CHOLMOD_OK (0), negative for a fatal error, and positive for
* a warning. file is a string containing the name of the source code
* file where the error occured, and line is the line number in that
* file. message is a string describing the error in more detail. */
/* ---------------------------------------------------------------------- */
/* ordering options */
/* ---------------------------------------------------------------------- */
/* The cholmod_analyze routine can try many different orderings and select
* the best one. It can also try one ordering method multiple times, with
* different parameter settings. The default is to use three orderings,
* the user's permutation (if provided), AMD which is the fastest ordering
* and generally gives good fill-in, and METIS. CHOLMOD's nested dissection
* (METIS with a constrained AMD) usually gives a better ordering than METIS
* alone (by about 5% to 10%) but it takes more time.
*
* If you know the method that is best for your matrix, set Common->nmethods
* to 1 and set Common->method [0] to the set of parameters for that method.
* If you set it to 1 and do not provide a permutation, then only AMD will
* be called.
*
* If METIS is not available, the default # of methods tried is 2 (the user
* permutation, if any, and AMD).
*
* To try other methods, set Common->nmethods to the number of methods you
* want to try. The suite of default methods and their parameters is
* described in the cholmod_defaults routine, and summarized here:
*
* Common->method [i]:
* i = 0: user-provided ordering (cholmod_analyze_p only)
* i = 1: AMD (for both A and A*A')
* i = 2: METIS
* i = 3: CHOLMOD's nested dissection (NESDIS), default parameters
* i = 4: natural
* i = 5: NESDIS with nd_small = 20000
* i = 6: NESDIS with nd_small = 4, no constrained minimum degree
* i = 7: NESDIS with no dense node removal
* i = 8: AMD for A, COLAMD for A*A'
*
* You can modify the suite of methods you wish to try by modifying
* Common.method [...] after calling cholmod_start or cholmod_defaults.
*
* For example, to use AMD, followed by a weighted postordering:
*
* Common->nmethods = 1 ;
* Common->method [0].ordering = CHOLMOD_AMD ;
* Common->postorder = TRUE ;
*
* To use the natural ordering (with no postordering):
*
* Common->nmethods = 1 ;
* Common->method [0].ordering = CHOLMOD_NATURAL ;
* Common->postorder = FALSE ;
*
* If you are going to factorize hundreds or more matrices with the same
* nonzero pattern, you may wish to spend a great deal of time finding a
* good permutation. In this case, try setting Common->nmethods to 9.
* The time spent in cholmod_analysis will be very high, but you need to
* call it only once.
*
* cholmod_analyze sets Common->current to a value between 0 and nmethods-1.
* Each ordering method uses the set of options defined by this parameter.
*/
int nmethods ; /* The number of ordering methods to try. Default: 0.
* nmethods = 0 is a special case. cholmod_analyze
* will try the user-provided ordering (if given) and AMD. Let fl and
* lnz be the flop count and nonzeros in L from AMD's ordering. Let
* anz be the number of nonzeros in the upper or lower triangular part
* of the symmetric matrix A. If fl/lnz < 500 or lnz/anz < 5, then this
* is a good ordering, and METIS is not attempted. Otherwise, METIS is
* tried. The best ordering found is used. If nmethods > 0, the
* methods used are given in the method[ ] array, below. The first
* three methods in the default suite of orderings is (1) use the given
* permutation (if provided), (2) use AMD, and (3) use METIS. Maximum
* allowed value is CHOLMOD_MAXMETHODS. */
int current ; /* The current method being tried. Default: 0. Valid
* range is 0 to nmethods-1. */
int selected ; /* The best method found. */
/* The suite of ordering methods and parameters: */
struct cholmod_method_struct
{
/* statistics for this method */
double lnz ; /* nnz(L) excl. zeros from supernodal amalgamation,
* for a "pure" L */
double fl ; /* flop count for a "pure", real simplicial LL'
* factorization, with no extra work due to
* amalgamation. Subtract n to get the LDL' flop count. Multiply
* by about 4 if the matrix is complex or zomplex. */
/* ordering method parameters */
double prune_dense ;/* dense row/col control for AMD, SYMAMD, CSYMAMD,
* and NESDIS (cholmod_nested_dissection). For a
* symmetric n-by-n matrix, rows/columns with more than
* MAX (16, prune_dense * sqrt (n)) entries are removed prior to
* ordering. They appear at the end of the re-ordered matrix.
*
* If prune_dense < 0, only completely dense rows/cols are removed.
*
* This paramater is also the dense column control for COLAMD and
* CCOLAMD. For an m-by-n matrix, columns with more than
* MAX (16, prune_dense * sqrt (MIN (m,n))) entries are removed prior
* to ordering. They appear at the end of the re-ordered matrix.
* CHOLMOD factorizes A*A', so it calls COLAMD and CCOLAMD with A',
* not A. Thus, this parameter affects the dense *row* control for
* CHOLMOD's matrix, and the dense *column* control for COLAMD and
* CCOLAMD.
*
* Removing dense rows and columns improves the run-time of the
* ordering methods. It has some impact on ordering quality
* (usually minimal, sometimes good, sometimes bad).
*
* Default: 10. */
double prune_dense2 ;/* dense row control for COLAMD and CCOLAMD.
* Rows with more than MAX (16, dense2 * sqrt (n))
* for an m-by-n matrix are removed prior to ordering. CHOLMOD's
* matrix is transposed before ordering it with COLAMD or CCOLAMD,
* so this controls the dense *columns* of CHOLMOD's matrix, and
* the dense *rows* of COLAMD's or CCOLAMD's matrix.
*
* If prune_dense2 < 0, only completely dense rows/cols are removed.
*
* Default: -1. Note that this is not the default for COLAMD and
* CCOLAMD. -1 is best for Cholesky. 10 is best for LU. */
double nd_oksep ; /* in NESDIS, when a node separator is computed, it
* discarded if nsep >= nd_oksep*n, where nsep is
* the number of nodes in the separator, and n is the size of the
* graph being cut. Valid range is 0 to 1. If 1 or greater, the
* separator is discarded if it consists of the entire graph.
* Default: 1 */
double other1 [4] ; /* future expansion */
size_t nd_small ; /* do not partition graphs with fewer nodes than
* nd_small, in NESDIS. Default: 200 (same as
* METIS) */
size_t other2 [4] ; /* future expansion */
int aggressive ; /* Aggresive absorption in AMD, COLAMD, SYMAMD,
* CCOLAMD, and CSYMAMD. Default: TRUE */
int order_for_lu ; /* CCOLAMD can be optimized to produce an ordering
* for LU or Cholesky factorization. CHOLMOD only
* performs a Cholesky factorization. However, you may wish to use
* CHOLMOD as an interface for CCOLAMD but use it for your own LU
* factorization. In this case, order_for_lu should be set to FALSE.
* When factorizing in CHOLMOD itself, you should *** NEVER *** set
* this parameter FALSE. Default: TRUE. */
int nd_compress ; /* If TRUE, compress the graph and subgraphs before
* partitioning them in NESDIS. Default: TRUE */
int nd_camd ; /* If 1, follow the nested dissection ordering
* with a constrained minimum degree ordering that
* respects the partitioning just found (using CAMD). If 2, use
* CSYMAMD instead. If you set nd_small very small, you may not need
* this ordering, and can save time by setting it to zero (no
* constrained minimum degree ordering). Default: 1. */
int nd_components ; /* The nested dissection ordering finds a node
* separator that splits the graph into two parts,
* which may be unconnected. If nd_components is TRUE, each of
* these connected components is split independently. If FALSE,
* each part is split as a whole, even if it consists of more than
* one connected component. Default: FALSE */
/* fill-reducing ordering to use */
int ordering ;
size_t other3 [4] ; /* future expansion */
} method [CHOLMOD_MAXMETHODS + 1] ;
int postorder ; /* If TRUE, cholmod_analyze follows the ordering with a
* weighted postorder of the elimination tree. Improves
* supernode amalgamation. Does not affect fundamental nnz(L) and
* flop count. Default: TRUE. */
/* ---------------------------------------------------------------------- */
/* memory management routines */
/* ---------------------------------------------------------------------- */
void *(*malloc_memory) (size_t) ; /* pointer to malloc */
void *(*realloc_memory) (void *, size_t) ; /* pointer to realloc */
void (*free_memory) (void *) ; /* pointer to free */
void *(*calloc_memory) (size_t, size_t) ; /* pointer to calloc */
/* ---------------------------------------------------------------------- */
/* routines for complex arithmetic */
/* ---------------------------------------------------------------------- */
int (*complex_divide) (double ax, double az, double bx, double bz,
double *cx, double *cz) ;
/* flag = complex_divide (ax, az, bx, bz, &cx, &cz) computes the complex
* division c = a/b, where ax and az hold the real and imaginary part
* of a, and b and c are stored similarly. flag is returned as 1 if
* a divide-by-zero occurs, or 0 otherwise. By default, the function
* pointer Common->complex_divide is set equal to cholmod_divcomplex.
*/
double (*hypotenuse) (double x, double y) ;
/* s = hypotenuse (x,y) computes s = sqrt (x*x + y*y), but does so more
* accurately. By default, the function pointer Common->hypotenuse is
* set equal to cholmod_hypot. See also the hypot function in the C99
* standard, which has an identical syntax and function. If you have
* a C99-compliant compiler, you can set Common->hypotenuse = hypot. */
/* ---------------------------------------------------------------------- */
/* METIS workarounds */
/* ---------------------------------------------------------------------- */
double metis_memory ; /* This is a parameter for CHOLMOD's interface to
* METIS, not a parameter to METIS itself. METIS
* uses an amount of memory that is difficult to estimate precisely
* beforehand. If it runs out of memory, it terminates your program.
* All routines in CHOLMOD except for CHOLMOD's interface to METIS
* return an error status and safely return to your program if they run
* out of memory. To mitigate this problem, the CHOLMOD interface
* can allocate a single block of memory equal in size to an empirical
* upper bound of METIS's memory usage times the Common->metis_memory
* parameter, and then immediately free it. It then calls METIS. If
* this pre-allocation fails, it is possible that METIS will fail as
* well, and so CHOLMOD returns with an out-of-memory condition without
* calling METIS.
*
* METIS_NodeND (used in the CHOLMOD_METIS ordering option) with its
* default parameter settings typically uses about (4*nz+40n+4096)
* times sizeof(int) memory, where nz is equal to the number of entries
* in A for the symmetric case or AA' if an unsymmetric matrix is
* being ordered (where nz includes both the upper and lower parts
* of A or AA'). The observed "upper bound" (with 2 exceptions),
* measured in an instrumented copy of METIS 4.0.1 on thousands of
* matrices, is (10*nz+50*n+4096) * sizeof(int). Two large matrices
* exceeded this bound, one by almost a factor of 2 (Gupta/gupta2).
*
* If your program is terminated by METIS, try setting metis_memory to
* 2.0, or even higher if needed. By default, CHOLMOD assumes that METIS
* does not have this problem (so that CHOLMOD will work correctly when
* this issue is fixed in METIS). Thus, the default value is zero.
* This work-around is not guaranteed anyway.
*
* If a matrix exceeds this predicted memory usage, AMD is attempted
* instead. It, too, may run out of memory, but if it does so it will
* not terminate your program.
*/
double metis_dswitch ; /* METIS_NodeND in METIS 4.0.1 gives a seg */
size_t metis_nswitch ; /* fault with one matrix of order n = 3005 and
* nz = 6,036,025. This is a very dense graph.
* The workaround is to use AMD instead of METIS for matrices of dimension
* greater than Common->metis_nswitch (default 3000) or more and with
* density of Common->metis_dswitch (default 0.66) or more.
* cholmod_nested_dissection has no problems with the same matrix, even
* though it uses METIS_NodeComputeSeparator on this matrix. If this
* seg fault does not affect you, set metis_nswitch to zero or less,
* and CHOLMOD will not switch to AMD based just on the density of the
* matrix (it will still switch to AMD if the metis_memory parameter
* causes the switch).
*/
/* ---------------------------------------------------------------------- */
/* workspace */
/* ---------------------------------------------------------------------- */
/* CHOLMOD has several routines that take less time than the size of
* workspace they require. Allocating and initializing the workspace would
* dominate the run time, unless workspace is allocated and initialized
* just once. CHOLMOD allocates this space when needed, and holds it here
* between calls to CHOLMOD. cholmod_start sets these pointers to NULL
* (which is why it must be the first routine called in CHOLMOD).
* cholmod_finish frees the workspace (which is why it must be the last
* call to CHOLMOD).
*/
size_t nrow ; /* size of Flag and Head */
SuiteSparse_long mark ; /* mark value for Flag array */
size_t iworksize ; /* size of Iwork. Upper bound: 6*nrow+ncol */
size_t xworksize ; /* size of Xwork, in bytes.
* maxrank*nrow*sizeof(double) for update/downdate.
* 2*nrow*sizeof(double) otherwise */
/* initialized workspace: contents needed between calls to CHOLMOD */
void *Flag ; /* size nrow, an integer array. Kept cleared between
* calls to cholmod rouines (Flag [i] < mark) */
void *Head ; /* size nrow+1, an integer array. Kept cleared between
* calls to cholmod routines (Head [i] = EMPTY) */
void *Xwork ; /* a double array. Its size varies. It is nrow for
* most routines (cholmod_rowfac, cholmod_add,
* cholmod_aat, cholmod_norm, cholmod_ssmult) for the real case, twice
* that when the input matrices are complex or zomplex. It is of size
* 2*nrow for cholmod_rowadd and cholmod_rowdel. For cholmod_updown,
* its size is maxrank*nrow where maxrank is 2, 4, or 8. Kept cleared
* between calls to cholmod (set to zero). */
/* uninitialized workspace, contents not needed between calls to CHOLMOD */
void *Iwork ; /* size iworksize, 2*nrow+ncol for most routines,
* up to 6*nrow+ncol for cholmod_analyze. */
int itype ; /* If CHOLMOD_LONG, Flag, Head, and Iwork are SuiteSparse_long.
* Otherwise all three arrays are int. */
int dtype ; /* double or float */
/* Common->itype and Common->dtype are used to define the types of all
* sparse matrices, triplet matrices, dense matrices, and factors
* created using this Common struct. The itypes and dtypes of all
* parameters to all CHOLMOD routines must match. */
int no_workspace_reallocate ; /* this is an internal flag, used as a
* precaution by cholmod_analyze. It is normally false. If true,
* cholmod_allocate_work is not allowed to reallocate any workspace;
* they must use the existing workspace in Common (Iwork, Flag, Head,
* and Xwork). Added for CHOLMOD v1.1 */
/* ---------------------------------------------------------------------- */
/* statistics */
/* ---------------------------------------------------------------------- */
/* fl and lnz are set only in cholmod_analyze and cholmod_rowcolcounts,
* in the Cholesky modudle. modfl is set only in the Modify module. */
int status ; /* error code */
double fl ; /* LL' flop count from most recent analysis */
double lnz ; /* fundamental nz in L */
double anz ; /* nonzeros in tril(A) if A is symmetric/lower,
* triu(A) if symmetric/upper, or tril(A*A') if
* unsymmetric, in last call to cholmod_analyze. */
double modfl ; /* flop count from most recent update/downdate/
* rowadd/rowdel (excluding flops to modify the
* solution to Lx=b, if computed) */
size_t malloc_count ; /* # of objects malloc'ed minus the # free'd*/
size_t memory_usage ; /* peak memory usage in bytes */
size_t memory_inuse ; /* current memory usage in bytes */
double nrealloc_col ; /* # of column reallocations */
double nrealloc_factor ;/* # of factor reallocations due to col. reallocs */
double ndbounds_hit ; /* # of times diagonal modified by dbound */
double rowfacfl ; /* # of flops in last call to cholmod_rowfac */
double aatfl ; /* # of flops to compute A(:,f)*A(:,f)' */
/* ---------------------------------------------------------------------- */
/* future expansion */
/* ---------------------------------------------------------------------- */
/* To allow CHOLMOD to be updated without recompiling the user application,
* additional space is set aside here for future statistics, parameters,
* and workspace. Note: additional entries were added in v1.1 to the
* method array, above, and thus v1.0 and v1.1 are not binary compatible.
*
* v1.1 to the current version are binary compatible.
*/
/* ---------------------------------------------------------------------- */
double other1 [10] ;
double SPQR_xstat [4] ; /* for SuiteSparseQR statistics */
/* SuiteSparseQR control parameters: */
double SPQR_grain ; /* task size is >= max (total flops / grain) */
double SPQR_small ; /* task size is >= small */
/* ---------------------------------------------------------------------- */
SuiteSparse_long SPQR_istat [10] ; /* for SuiteSparseQR statistics */
SuiteSparse_long other2 [6] ; /* reduced from size 16 in v1.6 */
/* ---------------------------------------------------------------------- */
int other3 [10] ; /* reduced from size 16 in v1.1. */
int prefer_binary ; /* cholmod_read_triplet converts a symmetric
* pattern-only matrix into a real matrix. If
* prefer_binary is FALSE, the diagonal entries are set to 1 + the degree
* of the row/column, and off-diagonal entries are set to -1 (resulting
* in a positive definite matrix if the diagonal is zero-free). Most
* symmetric patterns are the pattern a positive definite matrix. If
* this parameter is TRUE, then the matrix is returned with a 1 in each
* entry, instead. Default: FALSE. Added in v1.3. */
/* control parameter (added for v1.2): */
int default_nesdis ; /* Default: FALSE. If FALSE, then the default
* ordering strategy (when Common->nmethods == 0)
* is to try the given ordering (if present), AMD, and then METIS if AMD
* reports high fill-in. If Common->default_nesdis is TRUE then NESDIS
* is used instead in the default strategy. */
/* statistic (added for v1.2): */
int called_nd ; /* TRUE if the last call to
* cholmod_analyze called NESDIS or METIS. */
int blas_ok ; /* FALSE if BLAS int overflow; TRUE otherwise */
/* SuiteSparseQR control parameters: */
int SPQR_shrink ; /* controls stack realloc method */
int SPQR_nthreads ; /* number of TBB threads, 0 = auto */
/* ---------------------------------------------------------------------- */
size_t other4 [16] ;
/* ---------------------------------------------------------------------- */
void *other5 [16] ;
} cholmod_common ;
// in ../../src/CHOLMOD/Include/cholmod_core.h skip forward to - line 1170 :
/* A sparse matrix stored in compressed-column form. */
typedef struct cholmod_sparse_struct
{
size_t nrow ; /* the matrix is nrow-by-ncol */
size_t ncol ;
size_t nzmax ; /* maximum number of entries in the matrix */
/* pointers to int or SuiteSparse_long: */
void *p ; /* p [0..ncol], the column pointers */
void *i ; /* i [0..nzmax-1], the row indices */
/* for unpacked matrices only: */
void *nz ; /* nz [0..ncol-1], the # of nonzeros in each col. In
* packed form, the nonzero pattern of column j is in
* A->i [A->p [j] ... A->p [j+1]-1]. In unpacked form, column j is in
* A->i [A->p [j] ... A->p [j]+A->nz[j]-1] instead. In both cases, the
* numerical values (if present) are in the corresponding locations in
* the array x (or z if A->xtype is CHOLMOD_ZOMPLEX). */
/* pointers to double or float: */
void *x ; /* size nzmax or 2*nzmax, if present */
void *z ; /* size nzmax, if present */
int stype ; /* Describes what parts of the matrix are considered:
*
* 0: matrix is "unsymmetric": use both upper and lower triangular parts
* (the matrix may actually be symmetric in pattern and value, but
* both parts are explicitly stored and used). May be square or
* rectangular.
* >0: matrix is square and symmetric, use upper triangular part.
* Entries in the lower triangular part are ignored.
* <0: matrix is square and symmetric, use lower triangular part.
* Entries in the upper triangular part are ignored.
*
* Note that stype>0 and stype<0 are different for cholmod_sparse and
* cholmod_triplet. See the cholmod_triplet data structure for more
* details.
*/
int itype ; /* CHOLMOD_INT: p, i, and nz are int.
* CHOLMOD_INTLONG: p is SuiteSparse_long,
* i and nz are int.
* CHOLMOD_LONG: p, i, and nz are SuiteSparse_long */
int xtype ; /* pattern, real, complex, or zomplex */
int dtype ; /* x and z are double or float */
int sorted ; /* TRUE if columns are sorted, FALSE otherwise */
int packed ; /* TRUE if packed (nz ignored), FALSE if unpacked
* (nz is required) */
} cholmod_sparse ;
// in ../../src/CHOLMOD/Include/cholmod_core.h skip forward to - line 1552 :
/* A symbolic and numeric factorization, either simplicial or supernodal.
* In all cases, the row indices in the columns of L are kept sorted. */
typedef struct cholmod_factor_struct
{
/* ---------------------------------------------------------------------- */
/* for both simplicial and supernodal factorizations */
/* ---------------------------------------------------------------------- */
size_t n ; /* L is n-by-n */
size_t minor ; /* If the factorization failed, L->minor is the column
* at which it failed (in the range 0 to n-1). A value
* of n means the factorization was successful or
* the matrix has not yet been factorized. */
/* ---------------------------------------------------------------------- */
/* symbolic ordering and analysis */
/* ---------------------------------------------------------------------- */
void *Perm ; /* size n, permutation used */
void *ColCount ; /* size n, column counts for simplicial L */
void *IPerm ; /* size n, inverse permutation. Only created by
* cholmod_solve2 if Bset is used. */
/* ---------------------------------------------------------------------- */
/* simplicial factorization */
/* ---------------------------------------------------------------------- */
size_t nzmax ; /* size of i and x */
void *p ; /* p [0..ncol], the column pointers */
void *i ; /* i [0..nzmax-1], the row indices */
void *x ; /* x [0..nzmax-1], the numerical values */
void *z ;
void *nz ; /* nz [0..ncol-1], the # of nonzeros in each column.
* i [p [j] ... p [j]+nz[j]-1] contains the row indices,
* and the numerical values are in the same locatins
* in x. The value of i [p [k]] is always k. */
void *next ; /* size ncol+2. next [j] is the next column in i/x */
void *prev ; /* size ncol+2. prev [j] is the prior column in i/x.
* head of the list is ncol+1, and the tail is ncol. */
/* ---------------------------------------------------------------------- */
/* supernodal factorization */
/* ---------------------------------------------------------------------- */
/* Note that L->x is shared with the simplicial data structure. L->x has
* size L->nzmax for a simplicial factor, and size L->xsize for a supernodal
* factor. */
size_t nsuper ; /* number of supernodes */
size_t ssize ; /* size of s, integer part of supernodes */
size_t xsize ; /* size of x, real part of supernodes */
size_t maxcsize ; /* size of largest update matrix */
size_t maxesize ; /* max # of rows in supernodes, excl. triangular part */
void *super ; /* size nsuper+1, first col in each supernode */
void *pi ; /* size nsuper+1, pointers to integer patterns */
void *px ; /* size nsuper+1, pointers to real parts */
void *s ; /* size ssize, integer part of supernodes */
/* ---------------------------------------------------------------------- */
/* factorization type */
/* ---------------------------------------------------------------------- */
int ordering ; /* ordering method used */
int is_ll ; /* TRUE if LL', FALSE if LDL' */
int is_super ; /* TRUE if supernodal, FALSE if simplicial */
int is_monotonic ; /* TRUE if columns of L appear in order 0..n-1.
* Only applicable to simplicial numeric types. */
/* There are 8 types of factor objects that cholmod_factor can represent
* (only 6 are used):
*
* Numeric types (xtype is not CHOLMOD_PATTERN)
* --------------------------------------------
*
* simplicial LDL': (is_ll FALSE, is_super FALSE). Stored in compressed
* column form, using the simplicial components above (nzmax, p, i,
* x, z, nz, next, and prev). The unit diagonal of L is not stored,
* and D is stored in its place. There are no supernodes.
*
* simplicial LL': (is_ll TRUE, is_super FALSE). Uses the same storage
* scheme as the simplicial LDL', except that D does not appear.
* The first entry of each column of L is the diagonal entry of
* that column of L.
*
* supernodal LDL': (is_ll FALSE, is_super TRUE). Not used.
* FUTURE WORK: add support for supernodal LDL'
*
* supernodal LL': (is_ll TRUE, is_super TRUE). A supernodal factor,
* using the supernodal components described above (nsuper, ssize,
* xsize, maxcsize, maxesize, super, pi, px, s, x, and z).
*
*
* Symbolic types (xtype is CHOLMOD_PATTERN)
* -----------------------------------------
*
* simplicial LDL': (is_ll FALSE, is_super FALSE). Nothing is present
* except Perm and ColCount.
*
* simplicial LL': (is_ll TRUE, is_super FALSE). Identical to the
* simplicial LDL', except for the is_ll flag.
*
* supernodal LDL': (is_ll FALSE, is_super TRUE). Not used.
* FUTURE WORK: add support for supernodal LDL'
*
* supernodal LL': (is_ll TRUE, is_super TRUE). A supernodal symbolic
* factorization. The simplicial symbolic information is present
* (Perm and ColCount), as is all of the supernodal factorization
* except for the numerical values (x and z).
*/
int itype ; /* The integer arrays are Perm, ColCount, p, i, nz,
* next, prev, super, pi, px, and s. If itype is
* CHOLMOD_INT, all of these are int arrays.
* CHOLMOD_INTLONG: p, pi, px are SuiteSparse_long, others int.
* CHOLMOD_LONG: all integer arrays are SuiteSparse_long. */
int xtype ; /* pattern, real, complex, or zomplex */
int dtype ; /* x and z double or float */
} cholmod_factor ;
// in ../../src/CHOLMOD/Include/cholmod_core.h skip forward to - line 1836 :
/* A dense matrix in column-oriented form. It has no itype since it contains
* no integers. Entry in row i and column j is located in x [i+j*d].
*/
typedef struct cholmod_dense_struct
{
size_t nrow ; /* the matrix is nrow-by-ncol */
size_t ncol ;
size_t nzmax ; /* maximum number of entries in the matrix */
size_t d ; /* leading dimension (d >= nrow must hold) */
void *x ; /* size nzmax or 2*nzmax, if present */
void *z ; /* size nzmax, if present */
int xtype ; /* pattern, real, complex, or zomplex */
int dtype ; /* x and z double or float */
} cholmod_dense ;
// in ../../src/CHOLMOD/Include/cholmod_core.h skip forward to - line 2033 :
/* A sparse matrix stored in triplet form. */
typedef struct cholmod_triplet_struct
{
size_t nrow ; /* the matrix is nrow-by-ncol */
size_t ncol ;
size_t nzmax ; /* maximum number of entries in the matrix */
size_t nnz ; /* number of nonzeros in the matrix */
void *i ; /* i [0..nzmax-1], the row indices */
void *j ; /* j [0..nzmax-1], the column indices */
void *x ; /* size nzmax or 2*nzmax, if present */
void *z ; /* size nzmax, if present */
int stype ; /* Describes what parts of the matrix are considered:
*
* 0: matrix is "unsymmetric": use both upper and lower triangular parts
* (the matrix may actually be symmetric in pattern and value, but
* both parts are explicitly stored and used). May be square or
* rectangular.
* >0: matrix is square and symmetric. Entries in the lower triangular
* part are transposed and added to the upper triangular part when
* the matrix is converted to cholmod_sparse form.
* <0: matrix is square and symmetric. Entries in the upper triangular
* part are transposed and added to the lower triangular part when
* the matrix is converted to cholmod_sparse form.
*
* Note that stype>0 and stype<0 are different for cholmod_sparse and
* cholmod_triplet. The reason is simple. You can permute a symmetric
* triplet matrix by simply replacing a row and column index with their
* new row and column indices, via an inverse permutation. Suppose
* P = L->Perm is your permutation, and Pinv is an array of size n.
* Suppose a symmetric matrix A is represent by a triplet matrix T, with
* entries only in the upper triangular part. Then the following code:
*
* Ti = T->i ;
* Tj = T->j ;
* for (k = 0 ; k < n ; k++) Pinv [P [k]] = k ;
* for (k = 0 ; k < nz ; k++) Ti [k] = Pinv [Ti [k]] ;
* for (k = 0 ; k < nz ; k++) Tj [k] = Pinv [Tj [k]] ;
*
* creates the triplet form of C=P*A*P'. However, if T initially
* contains just the upper triangular entries (T->stype = 1), after
* permutation it has entries in both the upper and lower triangular
* parts. These entries should be transposed when constructing the
* cholmod_sparse form of A, which is what cholmod_triplet_to_sparse
* does. Thus:
*
* C = cholmod_triplet_to_sparse (T, 0, &Common) ;
*
* will return the matrix C = P*A*P'.
*
* Since the triplet matrix T is so simple to generate, it's quite easy
* to remove entries that you do not want, prior to converting T to the
* cholmod_sparse form. So if you include these entries in T, CHOLMOD
* assumes that there must be a reason (such as the one above). Thus,
* no entry in a triplet matrix is ever ignored.
*/
int itype ; /* CHOLMOD_LONG: i and j are SuiteSparse_long. Otherwise int */
int xtype ; /* pattern, real, complex, or zomplex */
int dtype ; /* x and z are double or float */
} cholmod_triplet ;
// -------- our (Matrix) short and const_ forms of of the pointers :
typedef cholmod_common* CHM_CM;
typedef cholmod_dense* CHM_DN;
typedef const cholmod_dense* const_CHM_DN;
typedef cholmod_factor* CHM_FR;
typedef const cholmod_factor* const_CHM_FR;
typedef cholmod_sparse* CHM_SP;
typedef const cholmod_sparse* const_CHM_SP;
typedef cholmod_triplet* CHM_TR;
typedef const cholmod_triplet* const_CHM_TR;