1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 2005-2021 The Octave Project Developers 4 // 5 // See the file COPYRIGHT.md in the top-level directory of this 6 // distribution or <https://octave.org/copyright/>. 7 // 8 // This file is part of Octave. 9 // 10 // Octave is free software: you can redistribute it and/or modify it 11 // under the terms of the GNU General Public License as published by 12 // the Free Software Foundation, either version 3 of the License, or 13 // (at your option) any later version. 14 // 15 // Octave is distributed in the hope that it will be useful, but 16 // WITHOUT ANY WARRANTY; without even the implied warranty of 17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 // GNU General Public License for more details. 19 // 20 // You should have received a copy of the GNU General Public License 21 // along with Octave; see the file COPYING. If not, see 22 // <https://www.gnu.org/licenses/>. 23 // 24 //////////////////////////////////////////////////////////////////////// 25 26 // This is the octave interface to ccolamd, which bore the copyright given 27 // in the help of the functions. 28 29 #if defined (HAVE_CONFIG_H) 30 # include "config.h" 31 #endif 32 33 #include <cstdlib> 34 35 #include "CSparse.h" 36 #include "Sparse.h" 37 #include "dNDArray.h" 38 #include "oct-locbuf.h" 39 #include "oct-sparse.h" 40 41 #include "defun.h" 42 #include "error.h" 43 #include "errwarn.h" 44 #include "ov.h" 45 #include "pager.h" 46 47 DEFUN (ccolamd, args, nargout, 48 doc: /* -*- texinfo -*- 49 @deftypefn {} {@var{p} =} ccolamd (@var{S}) 50 @deftypefnx {} {@var{p} =} ccolamd (@var{S}, @var{knobs}) 51 @deftypefnx {} {@var{p} =} ccolamd (@var{S}, @var{knobs}, @var{cmember}) 52 @deftypefnx {} {[@var{p}, @var{stats}] =} ccolamd (@dots{}) 53 54 Constrained column approximate minimum degree permutation. 55 56 @code{@var{p} = ccolamd (@var{S})} returns the column approximate minimum 57 degree permutation vector for the sparse matrix @var{S}. For a 58 non-symmetric matrix @var{S}, @code{@var{S}(:, @var{p})} tends to have 59 sparser LU@tie{}factors than @var{S}. 60 @code{chol (@var{S}(:, @var{p})' * @var{S}(:, @var{p}))} also tends to be 61 sparser than @code{chol (@var{S}' * @var{S})}. 62 @code{@var{p} = ccolamd (@var{S}, 1)} optimizes the ordering for 63 @code{lu (@var{S}(:, @var{p}))}. The ordering is followed by a column 64 elimination tree post-ordering. 65 66 @var{knobs} is an optional 1-element to 5-element input vector, with a 67 default value of @code{[0 10 10 1 0]} if not present or empty. Entries not 68 present are set to their defaults. 69 70 @table @code 71 @item @var{knobs}(1) 72 if nonzero, the ordering is optimized for @code{lu (S(:, p))}. It will be a 73 poor ordering for @code{chol (@var{S}(:, @var{p})' * @var{S}(:, @var{p}))}. 74 This is the most important knob for ccolamd. 75 76 @item @var{knobs}(2) 77 if @var{S} is m-by-n, rows with more than 78 @code{max (16, @var{knobs}(2) * sqrt (n))} entries are ignored. 79 80 @item @var{knobs}(3) 81 columns with more than 82 @code{max (16, @var{knobs}(3) * sqrt (min (@var{m}, @var{n})))} entries are 83 ignored and ordered last in the output permutation 84 (subject to the cmember constraints). 85 86 @item @var{knobs}(4) 87 if nonzero, aggressive absorption is performed. 88 89 @item @var{knobs}(5) 90 if nonzero, statistics and knobs are printed. 91 92 @end table 93 94 @var{cmember} is an optional vector of length @math{n}. It defines the 95 constraints on the column ordering. If @code{@var{cmember}(j) = @var{c}}, 96 then column @var{j} is in constraint set @var{c} (@var{c} must be in the 97 range 1 to n). In the output permutation @var{p}, all columns in set 1 98 appear first, followed by all columns in set 2, and so on. 99 @code{@var{cmember} = ones (1,n)} if not present or empty. 100 @code{ccolamd (@var{S}, [], 1 : n)} returns @code{1 : n} 101 102 @code{@var{p} = ccolamd (@var{S})} is about the same as 103 @code{@var{p} = colamd (@var{S})}. @var{knobs} and its default values 104 differ. @code{colamd} always does aggressive absorption, and it finds an 105 ordering suitable for both @code{lu (@var{S}(:, @var{p}))} and @code{chol 106 (@var{S}(:, @var{p})' * @var{S}(:, @var{p}))}; it cannot optimize its 107 ordering for @code{lu (@var{S}(:, @var{p}))} to the extent that 108 @code{ccolamd (@var{S}, 1)} can. 109 110 @var{stats} is an optional 20-element output vector that provides data 111 about the ordering and the validity of the input matrix @var{S}. Ordering 112 statistics are in @code{@var{stats}(1 : 3)}. @code{@var{stats}(1)} and 113 @code{@var{stats}(2)} are the number of dense or empty rows and columns 114 ignored by @sc{ccolamd} and @code{@var{stats}(3)} is the number of garbage 115 collections performed on the internal data structure used by @sc{ccolamd} 116 (roughly of size @code{2.2 * nnz (@var{S}) + 4 * @var{m} + 7 * @var{n}} 117 integers). 118 119 @code{@var{stats}(4 : 7)} provide information if CCOLAMD was able to 120 continue. The matrix is OK if @code{@var{stats}(4)} is zero, or 1 if 121 invalid. @code{@var{stats}(5)} is the rightmost column index that is 122 unsorted or contains duplicate entries, or zero if no such column exists. 123 @code{@var{stats}(6)} is the last seen duplicate or out-of-order row 124 index in the column index given by @code{@var{stats}(5)}, or zero if no 125 such row index exists. @code{@var{stats}(7)} is the number of duplicate 126 or out-of-order row indices. @code{@var{stats}(8 : 20)} is always zero in 127 the current version of @sc{ccolamd} (reserved for future use). 128 129 The authors of the code itself are @nospell{S. Larimore, T. Davis} and 130 @nospell{S. Rajamanickam} in collaboration with @nospell{J. Bilbert and E. Ng}. 131 Supported by the National Science Foundation 132 @nospell{(DMS-9504974, DMS-9803599, CCR-0203270)}, and a grant from 133 @nospell{Sandia} National Lab. 134 See @url{http://faculty.cse.tamu.edu/davis/suitesparse.html} for ccolamd, 135 csymamd, amd, colamd, symamd, and other related orderings. 136 @seealso{colamd, csymamd} 137 @end deftypefn */) 138 { 139 #if defined (HAVE_CCOLAMD) 140 141 int nargin = args.length (); 142 143 if (nargin < 1 || nargin > 3) 144 print_usage (); 145 146 octave_value_list retval (nargout == 2 ? 2 : 1); 147 int spumoni = 0; 148 149 // Get knobs 150 static_assert (CCOLAMD_KNOBS <= 40, 151 "ccolamd: # of CCOLAMD_KNOBS exceeded. Please report this to bugs.octave.org"); 152 double knob_storage[CCOLAMD_KNOBS]; 153 double *knobs = &knob_storage[0]; 154 CCOLAMD_NAME (_set_defaults) (knobs); 155 156 // Check for user-passed knobs 157 if (nargin > 1) 158 { 159 NDArray User_knobs = args(1).array_value (); 160 int nel_User_knobs = User_knobs.numel (); 161 162 if (nel_User_knobs > 0) 163 knobs[CCOLAMD_LU] = (User_knobs(0) != 0); 164 if (nel_User_knobs > 1) 165 knobs[CCOLAMD_DENSE_ROW] = User_knobs(1); 166 if (nel_User_knobs > 2) 167 knobs[CCOLAMD_DENSE_COL] = User_knobs(2); 168 if (nel_User_knobs > 3) 169 knobs[CCOLAMD_AGGRESSIVE] = (User_knobs(3) != 0); 170 if (nel_User_knobs > 4) 171 spumoni = (User_knobs(4) != 0); 172 173 // print knob settings if spumoni is set 174 if (spumoni) 175 { 176 octave_stdout << "\nccolamd version " << CCOLAMD_MAIN_VERSION << '.' 177 << CCOLAMD_SUB_VERSION << ", " << CCOLAMD_DATE 178 << ":\nknobs(1): " << User_knobs(0) << ", order for "; 179 if (knobs[CCOLAMD_LU] != 0) 180 octave_stdout << "lu (A)\n"; 181 else 182 octave_stdout << "chol (A'*A)\n"; 183 184 if (knobs[CCOLAMD_DENSE_ROW] >= 0) 185 octave_stdout << "knobs(2): " << User_knobs(1) 186 << ", rows with > max (16," 187 << knobs[CCOLAMD_DENSE_ROW] 188 << "*sqrt (size(A,2)))" 189 << " entries removed\n"; 190 else 191 octave_stdout << "knobs(2): " << User_knobs(1) 192 << ", no dense rows removed\n"; 193 194 if (knobs[CCOLAMD_DENSE_COL] >= 0) 195 octave_stdout << "knobs(3): " << User_knobs(2) 196 << ", cols with > max (16," 197 << knobs[CCOLAMD_DENSE_COL] << "*sqrt (size(A)))" 198 << " entries removed\n"; 199 else 200 octave_stdout << "knobs(3): " << User_knobs(2) 201 << ", no dense columns removed\n"; 202 203 if (knobs[CCOLAMD_AGGRESSIVE] != 0) 204 octave_stdout << "knobs(4): " << User_knobs(3) 205 << ", aggressive absorption: yes"; 206 else 207 octave_stdout << "knobs(4): " << User_knobs(3) 208 << ", aggressive absorption: no"; 209 210 octave_stdout << "knobs(5): " << User_knobs(4) 211 << ", statistics and knobs printed\n"; 212 } 213 } 214 215 octave_idx_type n_row, n_col, nnz; 216 octave_idx_type *ridx, *cidx; 217 SparseComplexMatrix scm; 218 SparseMatrix sm; 219 220 if (args(0).issparse ()) 221 { 222 if (args(0).iscomplex ()) 223 { 224 scm = args(0).sparse_complex_matrix_value (); 225 n_row = scm.rows (); 226 n_col = scm.cols (); 227 nnz = scm.nnz (); 228 ridx = scm.xridx (); 229 cidx = scm.xcidx (); 230 } 231 else 232 { 233 sm = args(0).sparse_matrix_value (); 234 235 n_row = sm.rows (); 236 n_col = sm.cols (); 237 nnz = sm.nnz (); 238 ridx = sm.xridx (); 239 cidx = sm.xcidx (); 240 } 241 } 242 else 243 { 244 if (args(0).iscomplex ()) 245 sm = SparseMatrix (real (args(0).complex_matrix_value ())); 246 else 247 sm = SparseMatrix (args(0).matrix_value ()); 248 249 n_row = sm.rows (); 250 n_col = sm.cols (); 251 nnz = sm.nnz (); 252 ridx = sm.xridx (); 253 cidx = sm.xcidx (); 254 } 255 256 // Allocate workspace for ccolamd 257 OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, p, n_col+1); 258 for (octave_idx_type i = 0; i < n_col+1; i++) 259 p[i] = cidx[i]; 260 261 octave_idx_type Alen = CCOLAMD_NAME (_recommended) (nnz, n_row, n_col); 262 OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, A, Alen); 263 for (octave_idx_type i = 0; i < nnz; i++) 264 A[i] = ridx[i]; 265 266 static_assert (CCOLAMD_STATS <= 40, 267 "ccolamd: # of CCOLAMD_STATS exceeded. Please report this to bugs.octave.org"); 268 octave::suitesparse_integer stats_storage[CCOLAMD_STATS]; 269 octave::suitesparse_integer *stats = &stats_storage[0]; 270 271 if (nargin > 2) 272 { 273 NDArray in_cmember = args(2).array_value (); 274 octave_idx_type cslen = in_cmember.numel (); 275 OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, cmember, cslen); 276 for (octave_idx_type i = 0; i < cslen; i++) 277 // convert cmember from 1-based to 0-based 278 cmember[i] = static_cast<octave::suitesparse_integer>(in_cmember(i) - 1); 279 280 if (cslen != n_col) 281 error ("ccolamd: CMEMBER must be of length equal to #cols of A"); 282 283 // Order the columns (destroys A) 284 if (! CCOLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats,cmember)) 285 { 286 CCOLAMD_NAME (_report) (stats); 287 288 error ("ccolamd: internal error!"); 289 } 290 } 291 else 292 { 293 // Order the columns (destroys A) 294 if (! CCOLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats, nullptr)) 295 { 296 CCOLAMD_NAME (_report) (stats); 297 298 error ("ccolamd: internal error!"); 299 } 300 } 301 302 // return the permutation vector 303 NDArray out_perm (dim_vector (1, n_col)); 304 for (octave_idx_type i = 0; i < n_col; i++) 305 out_perm(i) = p[i] + 1; 306 307 retval(0) = out_perm; 308 309 // print stats if spumoni > 0 310 if (spumoni > 0) 311 CCOLAMD_NAME (_report) (stats); 312 313 // Return the stats vector 314 if (nargout == 2) 315 { 316 NDArray out_stats (dim_vector (1, CCOLAMD_STATS)); 317 for (octave_idx_type i = 0 ; i < CCOLAMD_STATS ; i++) 318 out_stats(i) = stats[i]; 319 retval(1) = out_stats; 320 321 // fix stats (5) and (6), for 1-based information on 322 // jumbled matrix. note that this correction doesn't 323 // occur if symamd returns FALSE 324 out_stats(CCOLAMD_INFO1)++; 325 out_stats(CCOLAMD_INFO2)++; 326 } 327 328 return retval; 329 330 #else 331 332 octave_unused_parameter (args); 333 octave_unused_parameter (nargout); 334 335 err_disabled_feature ("ccolamd", "CCOLAMD"); 336 337 #endif 338 } 339 340 DEFUN (csymamd, args, nargout, 341 doc: /* -*- texinfo -*- 342 @deftypefn {} {@var{p} =} csymamd (@var{S}) 343 @deftypefnx {} {@var{p} =} csymamd (@var{S}, @var{knobs}) 344 @deftypefnx {} {@var{p} =} csymamd (@var{S}, @var{knobs}, @var{cmember}) 345 @deftypefnx {} {[@var{p}, @var{stats}] =} csymamd (@dots{}) 346 347 For a symmetric positive definite matrix @var{S}, return the permutation 348 vector @var{p} such that @code{@var{S}(@var{p},@var{p})} tends to have a 349 sparser Cholesky@tie{}factor than @var{S}. 350 351 Sometimes @code{csymamd} works well for symmetric indefinite matrices too. 352 The matrix @var{S} is assumed to be symmetric; only the strictly lower 353 triangular part is referenced. @var{S} must be square. The ordering is 354 followed by an elimination tree post-ordering. 355 356 @var{knobs} is an optional 1-element to 3-element input vector, with a 357 default value of @code{[10 1 0]}. Entries not present are set to their 358 defaults. 359 360 @table @code 361 @item @var{knobs}(1) 362 If @var{S} is n-by-n, then rows and columns with more than 363 @code{max(16,@var{knobs}(1)*sqrt(n))} entries are ignored, and ordered 364 last in the output permutation (subject to the cmember constraints). 365 366 @item @var{knobs}(2) 367 If nonzero, aggressive absorption is performed. 368 369 @item @var{knobs}(3) 370 If nonzero, statistics and knobs are printed. 371 372 @end table 373 374 @var{cmember} is an optional vector of length n. It defines the constraints 375 on the ordering. If @code{@var{cmember}(j) = @var{S}}, then row/column j is 376 in constraint set @var{c} (@var{c} must be in the range 1 to n). In the 377 output permutation @var{p}, rows/columns in set 1 appear first, followed 378 by all rows/columns in set 2, and so on. @code{@var{cmember} = ones (1,n)} 379 if not present or empty. @code{csymamd (@var{S},[],1:n)} returns 380 @code{1:n}. 381 382 @code{@var{p} = csymamd (@var{S})} is about the same as 383 @code{@var{p} = symamd (@var{S})}. @var{knobs} and its default values 384 differ. 385 386 @code{@var{stats}(4:7)} provide information if CCOLAMD was able to 387 continue. The matrix is OK if @code{@var{stats}(4)} is zero, or 1 if 388 invalid. @code{@var{stats}(5)} is the rightmost column index that is 389 unsorted or contains duplicate entries, or zero if no such column exists. 390 @code{@var{stats}(6)} is the last seen duplicate or out-of-order row 391 index in the column index given by @code{@var{stats}(5)}, or zero if no 392 such row index exists. @code{@var{stats}(7)} is the number of duplicate 393 or out-of-order row indices. @code{@var{stats}(8:20)} is always zero in 394 the current version of @sc{ccolamd} (reserved for future use). 395 396 The authors of the code itself are @nospell{S. Larimore, T. Davis} and 397 @nospell{S. Rajamanickam} in collaboration with @nospell{J. Bilbert and E. Ng}. 398 Supported by the National Science Foundation 399 @nospell{(DMS-9504974, DMS-9803599, CCR-0203270)}, and a grant from 400 @nospell{Sandia} National Lab. 401 See @url{http://faculty.cse.tamu.edu/davis/suitesparse.html} for ccolamd, 402 colamd, csymamd, amd, colamd, symamd, and other related orderings. 403 @seealso{symamd, ccolamd} 404 @end deftypefn */) 405 { 406 #if defined (HAVE_CCOLAMD) 407 408 int nargin = args.length (); 409 410 if (nargin < 1 || nargin > 3) 411 print_usage (); 412 413 octave_value_list retval (nargout == 2 ? 2 : 1); 414 int spumoni = 0; 415 416 // Get knobs 417 static_assert (CCOLAMD_KNOBS <= 40, 418 "csymamd: # of CCOLAMD_KNOBS exceeded. Please report this to bugs.octave.org"); 419 double knob_storage[CCOLAMD_KNOBS]; 420 double *knobs = &knob_storage[0]; 421 CCOLAMD_NAME (_set_defaults) (knobs); 422 423 // Check for user-passed knobs 424 if (nargin > 1) 425 { 426 NDArray User_knobs = args(1).array_value (); 427 int nel_User_knobs = User_knobs.numel (); 428 429 if (nel_User_knobs > 0) 430 knobs[CCOLAMD_DENSE_ROW] = User_knobs(0); 431 if (nel_User_knobs > 1) 432 knobs[CCOLAMD_AGGRESSIVE] = User_knobs(1); 433 if (nel_User_knobs > 2) 434 spumoni = static_cast<int> (User_knobs(2)); 435 436 // print knob settings if spumoni is set 437 if (spumoni) 438 { 439 octave_stdout << "\ncsymamd version " << CCOLAMD_MAIN_VERSION 440 << '.' << CCOLAMD_SUB_VERSION 441 << ", " << CCOLAMD_DATE << "\n"; 442 443 if (knobs[CCOLAMD_DENSE_ROW] >= 0) 444 octave_stdout << "knobs(1): " << User_knobs(0) 445 << ", rows/cols with > max (16," 446 << knobs[CCOLAMD_DENSE_ROW] 447 << "*sqrt (size(A,2)))" 448 << " entries removed\n"; 449 else 450 octave_stdout << "knobs(1): " << User_knobs(0) 451 << ", no dense rows/cols removed\n"; 452 453 if (knobs[CCOLAMD_AGGRESSIVE] != 0) 454 octave_stdout << "knobs(2): " << User_knobs(1) 455 << ", aggressive absorption: yes"; 456 else 457 octave_stdout << "knobs(2): " << User_knobs(1) 458 << ", aggressive absorption: no"; 459 460 octave_stdout << "knobs(3): " << User_knobs(2) 461 << ", statistics and knobs printed\n"; 462 } 463 } 464 465 octave_idx_type n_row, n_col; 466 octave_idx_type *ridx, *cidx; 467 SparseMatrix sm; 468 SparseComplexMatrix scm; 469 470 if (args(0).issparse ()) 471 { 472 if (args(0).iscomplex ()) 473 { 474 scm = args(0).sparse_complex_matrix_value (); 475 n_row = scm.rows (); 476 n_col = scm.cols (); 477 ridx = scm.xridx (); 478 cidx = scm.xcidx (); 479 } 480 else 481 { 482 sm = args(0).sparse_matrix_value (); 483 n_row = sm.rows (); 484 n_col = sm.cols (); 485 ridx = sm.xridx (); 486 cidx = sm.xcidx (); 487 } 488 } 489 else 490 { 491 if (args(0).iscomplex ()) 492 sm = SparseMatrix (real (args(0).complex_matrix_value ())); 493 else 494 sm = SparseMatrix (args(0).matrix_value ()); 495 496 n_row = sm.rows (); 497 n_col = sm.cols (); 498 ridx = sm.xridx (); 499 cidx = sm.xcidx (); 500 } 501 502 if (n_row != n_col) 503 err_square_matrix_required ("csymamd", "S"); 504 505 // Allocate workspace for symamd 506 OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, perm, n_col+1); 507 static_assert (CCOLAMD_STATS <= 40, 508 "csymamd: # of CCOLAMD_STATS exceeded. Please report this to bugs.octave.org"); 509 octave::suitesparse_integer stats_storage[CCOLAMD_STATS]; 510 octave::suitesparse_integer *stats = &stats_storage[0]; 511 512 if (nargin > 2) 513 { 514 NDArray in_cmember = args(2).array_value (); 515 octave_idx_type cslen = in_cmember.numel (); 516 OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, cmember, cslen); 517 for (octave_idx_type i = 0; i < cslen; i++) 518 // convert cmember from 1-based to 0-based 519 cmember[i] = static_cast<octave_idx_type> (in_cmember(i) - 1); 520 521 if (cslen != n_col) 522 error ("csymamd: CMEMBER must be of length equal to #cols of A"); 523 524 if (! CSYMAMD_NAME () (n_col, 525 octave::to_suitesparse_intptr (ridx), 526 octave::to_suitesparse_intptr (cidx), 527 perm, knobs, stats, &calloc, &free, cmember, -1)) 528 { 529 CSYMAMD_NAME (_report)(stats); 530 531 error ("csymamd: internal error!"); 532 } 533 } 534 else 535 { 536 if (! CSYMAMD_NAME () (n_col, 537 octave::to_suitesparse_intptr (ridx), 538 octave::to_suitesparse_intptr (cidx), 539 perm, knobs, stats, &calloc, &free, nullptr, -1)) 540 { 541 CSYMAMD_NAME (_report)(stats); 542 543 error ("csymamd: internal error!"); 544 } 545 } 546 547 // return the permutation vector 548 NDArray out_perm (dim_vector (1, n_col)); 549 for (octave_idx_type i = 0; i < n_col; i++) 550 out_perm(i) = perm[i] + 1; 551 552 retval(0) = out_perm; 553 554 // print stats if spumoni > 0 555 if (spumoni > 0) 556 CSYMAMD_NAME (_report)(stats); 557 558 // Return the stats vector 559 if (nargout == 2) 560 { 561 NDArray out_stats (dim_vector (1, CCOLAMD_STATS)); 562 for (octave_idx_type i = 0 ; i < CCOLAMD_STATS ; i++) 563 out_stats(i) = stats[i]; 564 retval(1) = out_stats; 565 566 // fix stats (5) and (6), for 1-based information on 567 // jumbled matrix. note that this correction doesn't 568 // occur if symamd returns FALSE 569 out_stats(CCOLAMD_INFO1)++; 570 out_stats(CCOLAMD_INFO2)++; 571 } 572 573 return retval; 574 575 #else 576 577 octave_unused_parameter (args); 578 octave_unused_parameter (nargout); 579 580 err_disabled_feature ("csymamd", "CCOLAMD"); 581 582 #endif 583 } 584