1 /* 2 HMat-OSS (HMatrix library, open source software) 3 4 Copyright (C) 2014-2015 Airbus Group SAS 5 6 This program is free software; you can redistribute it and/or 7 modify it under the terms of the GNU General Public License 8 as published by the Free Software Foundation; either version 2 9 of the License, or (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with this program; if not, write to the Free Software 18 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 20 http://github.com/jeromerobert/hmat-oss 21 */ 22 23 /*! \file 24 \ingroup HMatrix 25 \brief C interface to the HMatrix library. 26 */ 27 #ifndef _HMAT_H 28 #define _HMAT_H 29 30 #include <stdlib.h> 31 #include "hmat/config.h" 32 33 #if defined _WIN32 34 # define HMAT_HELPER_DLL_IMPORT __declspec(dllimport) 35 # define HMAT_HELPER_DLL_EXPORT __declspec(dllexport) 36 #else 37 # define HMAT_HELPER_DLL_IMPORT 38 # define HMAT_HELPER_DLL_EXPORT 39 #endif 40 41 #ifdef HMAT_STATIC /* We are building hmat as a static library */ 42 # define HMAT_API 43 #elif defined HMAT_DLL_EXPORTS /* We are building hmat as a shared library */ 44 # define HMAT_API HMAT_HELPER_DLL_EXPORT 45 #else /* We are using hmat library */ 46 # define HMAT_API HMAT_HELPER_DLL_IMPORT 47 #endif /* !HMAT_STATIC */ 48 49 #ifdef __cplusplus 50 extern "C" { 51 #endif 52 53 /*! \brief All the scalar types */ 54 typedef enum { 55 HMAT_SIMPLE_PRECISION=0, 56 HMAT_DOUBLE_PRECISION=1, 57 HMAT_SIMPLE_COMPLEX=2, 58 HMAT_DOUBLE_COMPLEX=3 59 } hmat_value_t; 60 61 /** Choice of the compression method */ 62 typedef enum { 63 hmat_compress_svd, 64 hmat_compress_aca_full, 65 hmat_compress_aca_partial, 66 hmat_compress_aca_plus, 67 hmat_compress_aca_random 68 } hmat_compress_t; 69 70 typedef enum { 71 hmat_block_full, 72 hmat_block_null, 73 hmat_block_sparse 74 } hmat_block_t; 75 76 typedef enum { 77 hmat_factorization_none = -1, 78 hmat_factorization_lu, 79 hmat_factorization_ldlt, 80 hmat_factorization_llt 81 } hmat_factorization_t; 82 83 typedef struct hmat_block_info_struct { 84 hmat_block_t block_type; 85 /** 86 * user data to pass from prepare function to compute function. 87 * Will also contains user data required to execute is_guaranteed_null_row and 88 * is_guaranteed_null_col 89 */ 90 void * user_data; 91 /*! \brief Function provided by the user in hmat_prepare_func_t and used to release user_data */ 92 void (*release_user_data)(void* user_data); 93 /*! \brief Function provided by the user in hmat_prepare_func_t and used to quickly detect null rows 94 95 It should be able to detect null rows very quickly, possibly without to compute them, 96 at the cost of possibly missing some null rows. It should return '\1' for a null row, 97 '\0' means that the row can be anything (null or non-null). 98 The block_row_offset argument is the row index within this block. 99 100 Note 1: We have to use hmat_block_info_struct in this definition, but 101 user prototype should be 102 char is_guaranteed_null_row(const hmat_block_info_t * block_info, int i, int stratum); 103 104 Note 2: It is convenient to store informations about null rows and columns for each stratum 105 inside the prepare function, so that this function only performs a lookup to quickly 106 determine whether this row is null or not. 107 */ 108 char (*is_guaranteed_null_row)(const struct hmat_block_info_struct * block_info, int block_row_offset, int stratum); 109 /*! \brief Function provided by the user in hmat_prepare_func_t and used to quickly detect null columns (equivalent to is_guaranteed_null_row) */ 110 char (*is_guaranteed_null_col)(const struct hmat_block_info_struct * block_info, int block_col_offset, int stratum); 111 /** 112 * The memory needed to assemble the block. 113 * When set to 0, the hmat_prepare_func_t should reset it 114 * to the expected value and return. The hmat_prepare_func_t will then 115 * be called a second time to run the actual preparation. 116 */ 117 size_t needed_memory; 118 /** the number of strata in the block */ 119 int number_of_strata; 120 } hmat_block_info_t; 121 122 /*! \brief Prepare block assembly. 123 \param row_start starting row 124 \param row_count number of rows 125 \param col_start starting column 126 \param col_count number of columns 127 \param row_hmat2client renumbered rows -> global row indices mapping 128 \param row_client2hmat global row indices -> renumbered rows mapping 129 \param col_hmat2client renumbered cols -> global col indices mapping 130 \param col_client2hmat global col indices -> renumbered cols mapping 131 \param context user provided data 132 */ 133 typedef void (*hmat_prepare_func_t)(int row_start, 134 int row_count, 135 int col_start, 136 int col_count, 137 int *row_hmat2client, 138 int *row_client2hmat, 139 int *col_hmat2client, 140 int *col_client2hmat, 141 void *context, 142 hmat_block_info_t * block_info); 143 144 /*! \brief Compute a sub-block. 145 146 \warning The computation has to be prepared with \a prepare_func() first. 147 148 The supported cases are: 149 - Compute the full block 150 - Compute exactly one row 151 - Compute exactly one column 152 153 Regarding the indexing: For all indices in this function, the "C" 154 conventions are followed, ie indices start at 0, the lower bound is 155 included and the upper bound is excluded. In contrast with 156 hmat_prepare_func_t, the indexing is relative 157 to the block, that is the row i is the i-th row within the block. 158 159 \param v_data opaque pointer, as set by \a prepare_func() in field user_data of hmat_block_info_t 160 \param block_row_start starting row in the block 161 \param block_row_count number of rows to be computed 162 \param block_col_start starting column in the block 163 \param block_col_count number of columns to be computed 164 \param block pointer to the output buffer. No padding is allowed, 165 that is the leading dimension of the buffer must be its real leading 166 dimension (1 for a row). Column-major order is assumed. 167 */ 168 typedef void (*hmat_compute_func_t)(void* v_data, int block_row_start, int block_row_count, 169 int block_col_start, int block_col_count, void* block); 170 171 /*! \brief Compute a single matrix term 172 173 \param user_context pointer to user data, see \a assemble_hmatrix_simple_interaction function 174 \param row row index 175 \param col column index 176 \param result address where result is stored; result is a pointer to a double for real matrices, 177 and a pointer to a double complex for complex matrices. 178 */ 179 typedef void (*hmat_interaction_func_t)(void* user_context, int row, int col, void* result); 180 181 typedef struct hmat_clustering_algorithm hmat_clustering_algorithm_t; 182 183 /* Opaque pointer */ 184 typedef struct hmat_cluster_tree_struct hmat_cluster_tree_t; 185 186 /* Median clustering */ 187 HMAT_API hmat_clustering_algorithm_t* hmat_create_clustering_median(void); 188 /* NTiles recursive clustering */ 189 HMAT_API hmat_clustering_algorithm_t* hmat_create_clustering_ntilesrecursive(int nTiles); 190 /* Geometric clustering */ 191 HMAT_API hmat_clustering_algorithm_t* hmat_create_clustering_geometric(void); 192 /* Hybrid clustering */ 193 HMAT_API hmat_clustering_algorithm_t* hmat_create_clustering_hybrid(void); 194 /* Create a new clustering algorithm by setting the maximum number of degrees of freedom in a leaf */ 195 HMAT_API hmat_clustering_algorithm_t* hmat_create_clustering_max_dof(const hmat_clustering_algorithm_t* algo, int max_dof); 196 197 /** 198 * Create a new clustering algorithm which put large span apart and 199 * delegate to an other algorithm. 200 * @param algo the delegate algorithm 201 * @param ratio ratio between the number of DOF in a cluster and 202 * a span size so a DOF is concidered as large 203 */ 204 HMAT_API hmat_clustering_algorithm_t* hmat_create_clustering_span( 205 const hmat_clustering_algorithm_t* algo, double ratio); 206 207 /* Create a new clustering algorithm (for tests purpose only) */ 208 HMAT_API hmat_clustering_algorithm_t* hmat_create_void_clustering(const hmat_clustering_algorithm_t* algo); 209 /* Create a new clustering algorithm (for tests purpose only) */ 210 HMAT_API hmat_clustering_algorithm_t* hmat_create_shuffle_clustering(const hmat_clustering_algorithm_t* algo, int from_divider, int to_divider); 211 /* Delete clustering algorithm */ 212 HMAT_API void hmat_delete_clustering(hmat_clustering_algorithm_t *algo); 213 /* Set the clustering divider parameter */ 214 HMAT_API void hmat_set_clustering_divider(hmat_clustering_algorithm_t* algo, int divider); 215 216 /*! \brief Create a ClusterTree from the DoFs coordinates. 217 218 \param coord DoFs coordinates 219 \param dimension spatial dimension 220 \param size number of DoFs 221 \param algo pointer to clustering algorithm 222 \return an opaque pointer to a ClusterTree, or NULL in case of error. 223 */ 224 HMAT_API hmat_cluster_tree_t * hmat_create_cluster_tree(double* coord, int dimension, int size, hmat_clustering_algorithm_t* algo); 225 226 /* Opaque pointer */ 227 typedef struct hmat_cluster_tree_builder hmat_cluster_tree_builder_t; 228 229 HMAT_API hmat_cluster_tree_builder_t* hmat_create_cluster_tree_builder(const hmat_clustering_algorithm_t* algo); 230 231 /* Specify an algorithm for nodes at given depth and below */ 232 HMAT_API void hmat_cluster_tree_builder_add_algorithm(hmat_cluster_tree_builder_t* ctb, int level, const hmat_clustering_algorithm_t* algo); 233 234 HMAT_API void hmat_delete_cluster_tree_builder(hmat_cluster_tree_builder_t* ctb); 235 236 /*! \brief Create a ClusterTree from the DoFs coordinates. 237 238 \param coord DoFs coordinates 239 \param dimension spatial dimension 240 \param size number of DoFs 241 \param ctb pointer to an opaque ClusterTreeBuilder 242 \return an opaque pointer to a ClusterTree, or NULL in case of error. 243 */ 244 HMAT_API hmat_cluster_tree_t * hmat_create_cluster_tree_from_builder(double* coord, int dimension, int size, const hmat_cluster_tree_builder_t* ctb); 245 246 HMAT_API void hmat_delete_cluster_tree(const hmat_cluster_tree_t * tree); 247 248 HMAT_API hmat_cluster_tree_t * hmat_copy_cluster_tree(const hmat_cluster_tree_t * tree); 249 250 struct hmat_cluster_tree_create_context_t { 251 /** Spatial dimension */ 252 unsigned dimension; 253 unsigned number_of_points; 254 double * coordinates; 255 unsigned number_of_dof; 256 /** 257 * Offset of each span in the spans array. 258 * if span_offset is NULL each span is considered as a single point. 259 * if not NULL the size of this array must be number_of_dof. 260 * span_offsets[0] is the offset of the dof 1 (not 0). 261 * span_offset[number_of_dof-1] is the length of the spans array. 262 */ 263 unsigned * span_offsets; 264 /** The id of the points in each span */ 265 unsigned * spans; 266 /** Group index of dofs (may be NULL) */ 267 int * group_index; 268 /** pointer to an opaque ClusterTreeBuilder */ 269 const hmat_cluster_tree_builder_t* builder; 270 }; 271 272 HMAT_API hmat_cluster_tree_t * hmat_create_cluster_tree_generic(struct hmat_cluster_tree_create_context_t*); 273 274 /*! 275 * Return the number of nodes in a cluster tree 276 */ 277 HMAT_API int hmat_tree_nodes_count(const hmat_cluster_tree_t * tree); 278 /*! 279 * Returns the index-th son of a cluster tree root 280 */ 281 HMAT_API hmat_cluster_tree_t *hmat_cluster_get_son( hmat_cluster_tree_t * tree, int index ); 282 283 /** Information on a cluster tree */ 284 typedef struct 285 { 286 /* ! Spatial dimension of degrees of freedom */ 287 int spatial_dimension; 288 289 /* ! Number of degrees of freedom */ 290 int dimension; 291 292 /* ! Number of tree nodes */ 293 size_t nr_tree_nodes; 294 } hmat_cluster_info_t; 295 296 HMAT_API int hmat_cluster_get_info(hmat_cluster_tree_t *tree, hmat_cluster_info_t* info); 297 298 /** 299 * @brief Return the renumbering of the cluster tree 300 * @param tree a cluster tree 301 * @return The original position (before renumbering) of the ith value in the cluster 302 * tree (after renumbering). This array must not be freed nor modified by the caller. 303 */ 304 HMAT_API const int * hmat_cluster_get_indices(const hmat_cluster_tree_t *tree); 305 306 typedef struct { 307 /** eta for Hackbusch condition */ 308 double eta; 309 /** ratio (in [0, 0.5]) to prevent tall and skinny blocks */ 310 double ratio; 311 /** maximum block width */ 312 size_t max_width; 313 } hmat_admissibility_param_t; 314 315 /** Init an hmat_admissibility_param structure with default values */ 316 HMAT_API void hmat_init_admissibility_param(hmat_admissibility_param_t *); 317 318 /* Opaque pointer */ 319 typedef struct hmat_admissibility_condition hmat_admissibility_t; 320 321 /** Create an admissibility condition from parameters */ 322 HMAT_API hmat_admissibility_t* hmat_create_admissibility(hmat_admissibility_param_t *); 323 324 /** Update an admissibility condition with parameters */ 325 HMAT_API void hmat_update_admissibility(hmat_admissibility_t*, hmat_admissibility_param_t *); 326 327 /* Create a standard (Hackbusch) admissibility condition, with a given eta */ 328 HMAT_API hmat_admissibility_t* hmat_create_admissibility_standard(double eta); 329 330 /** 331 * @brief Create an admissibility condiction which set all blocks as admissible 332 * @param max_block_size The maximum acceptable block size in number of values (rows * cols) 333 * @param min_nr_block The minimum acceptable number of blocks created with this condition 334 * @param split_rows Tel whether or not to split rows 335 * @param split_cols Tel whether or not to split cols 336 */ 337 HMAT_API hmat_admissibility_t* hmat_create_admissibility_always( 338 size_t max_size, unsigned int min_block, int split_rows, int split_cols); 339 340 /** 341 * @brief Create an admissibility condiction which set all blocks as full 342 * @param max_block_size The maximum acceptable block size in number of values (rows * cols). 343 * Use 0 to switch to autodetect. 344 * @param min_nr_block The minimum acceptable number of blocks created with this condition 345 * Use 0 to switch to autodetect. 346 * @param split_rows Tel whether or not to split rows 347 * @param split_cols Tel whether or not to split cols 348 */ 349 HMAT_API hmat_admissibility_t* hmat_create_admissibility_never( 350 size_t max_size, unsigned int min_block, int split_rows, int split_cols); 351 352 /* Delete admissibility condition */ 353 HMAT_API void hmat_delete_admissibility(hmat_admissibility_t * cond); 354 355 /* Opaque pointer */ 356 typedef struct hmat_compression_algorithm_struct hmat_compression_algorithm_t; 357 358 /* Create a procedure to compress matrix blocks */ 359 HMAT_API hmat_compression_algorithm_t* hmat_create_compression_svd(double epsilon); 360 HMAT_API hmat_compression_algorithm_t* hmat_create_compression_aca_full(double epsilon); 361 HMAT_API hmat_compression_algorithm_t* hmat_create_compression_aca_partial(double epsilon); 362 HMAT_API hmat_compression_algorithm_t* hmat_create_compression_aca_plus(double epsilon); 363 HMAT_API hmat_compression_algorithm_t* hmat_create_compression_aca_random(double epsilon); 364 365 /* Delete a compression algorithm */ 366 HMAT_API void hmat_delete_compression(const hmat_compression_algorithm_t* algo); 367 368 /** Information on the HMatrix */ 369 typedef struct 370 { 371 /*! Number of allocated zeros */ 372 int full_zeros; 373 /** Total number of terms stored in the full leaves of the HMatrix */ 374 size_t full_size; 375 376 /** @deprecated Use compressed_size, uncompressed_size or full_size */ 377 size_t rk_size; 378 379 /** Total number of full leaves of the HMatrix */ 380 size_t full_count; 381 382 /** Total number of rk leaves of the HMatrix */ 383 size_t rk_count; 384 385 /** Total number of terms stored in the HMatrix */ 386 size_t compressed_size; 387 388 /*! Total number of terms that would be stored if the matrix was not compressed */ 389 size_t uncompressed_size; 390 391 /*! Total number of block cluster tree nodes in the HMatrix */ 392 int nr_block_clusters; 393 394 /*! Number of rows in the largest Rk matrice with rows + cols criteria */ 395 int largest_rk_dim_rows; 396 /*! Number of cols in the largest Rk matrice with rows + cols criteria */ 397 int largest_rk_dim_cols; 398 /*! Number of rows in the largest Rk matrice with memory criteria */ 399 int largest_rk_mem_rows; 400 /*! Number of cols in the largest Rk matrice with memory criteria */ 401 int largest_rk_mem_cols; 402 /*! Rank of the largest Rk matrice with memory criteria */ 403 int largest_rk_mem_rank; 404 } hmat_info_t; 405 406 typedef struct hmat_matrix_struct hmat_matrix_t; 407 408 /** Allow to implement a progress bar associated to assemble or factorize */ 409 typedef struct hmat_progress_struct { 410 int max; 411 int current; 412 /** Called each time assembling or factorization progress */ 413 void (*update)(struct hmat_progress_struct* context); 414 void * user_data; 415 } hmat_progress_t; 416 417 /** 418 * Return the default progress bar. 419 * This is a singleton which must/can not be freed. 420 */ 421 hmat_progress_t * hmat_default_progress(void); 422 423 /** 424 * Function representing a generic stream. 425 * It could be implemented as FILE*, unix fd, C++ iostream while using a C API. 426 * @param buffer the buffer to read or write 427 * @param n the size of the buffer 428 */ 429 typedef void (*hmat_iostream)(void * buffer, size_t n, void *user_data); 430 431 /** */ 432 struct hmat_block_compute_context_t { 433 /** 434 * opaque pointer as set by \a prepare_func() in 435 * hmat_block_info_t.user_data 436 */ 437 void* user_data; 438 /** Location and size of the block */ 439 int row_start, row_count, col_start, col_count; 440 /** stratum id */ 441 int stratum; 442 /** 443 * block pointer to the output buffer. No padding is allowed, 444 * that is the leading dimension of the buffer must be its real 445 * leading. 446 */ 447 void* block; 448 }; 449 /** 450 * Argument of the assemble_generic function. 451 * Only one of block_compute/advanced_compute, simple_compute or assembly can be non NULL. 452 * If both block_compute & advanced_compute are set, advanced is used. 453 */ 454 typedef struct { 455 /** 456 * The hmat_matrix_t matrix has previouly been initialized, it is a 457 * hierarchical structure containing inner nodes and leaves, which 458 * contain either a FullMatrix or an RkMatrix (compressed form). 459 * 460 * There are 4 different ways to perform a matrix assembly; 461 * the current recommended way is via advanced_compute, others 462 * are kept for legacy code. 463 * 464 * 1. assembly: this member is a pointer to an Assembly instance; 465 * it provides an assemble method: 466 * 467 * void assemble(const LocalSettings & settings, 468 * const ClusterTree & rows, const ClusterTree & cols, 469 * bool admissible, 470 * FullMatrix<T> * & fullMatrix, RkMatrix<T> * & rkMatrix, 471 * const AllocationObserver & ao) 472 * 473 * If admissible argument is true, this method must compute the rkMatrix 474 * argument, otherwise fullMatrix is computed. 475 * 476 * 2. simple_compute: this member is a pointer to a function 477 * 478 * void interaction(void* user_context, int row, int col, void* result) 479 * 480 * If block is full, interaction is called for all (row,col) tuple of 481 * this block, the first argument user_context contains user data to 482 * perform these computations. 483 * 484 * If block is compressed, the compression algorithm will either retrieve 485 * the full block (with SVD or ACA full algorithm) and compress it, or 486 * retrieve only some rows and columns (ACA partial or ACA+ algorithms). 487 * In all cases, it uses the interaction function. 488 * Note that row and col are numbered here according to user numbering. 489 * 490 * 3. block_compute: this member is a pointer to a function 491 * 492 * void compute(void* user_context, int block_row_start, int block_row_count, 493 * int block_col_start, int block_col_count, void* block) 494 * 495 * Moreover, user must also provide a prepare function. 496 * It is called to fill up an hmat_block_info_t structure, and potentially 497 * allocate hmat_block_info_t.user_data member (in which case 498 * hmat_block_info_t.release_user_data function pointer must be provided). 499 * If block_type is set to hmat_block_null in prepare function, nothing is 500 * computed. 501 * If block is full, compute is called on the whole block. 502 * 503 * If block is compressed, the compression algorithm will either retrieve 504 * the full block (with SVD or ACA full algorithm) and compress it, or 505 * retrieve only some rows and columns (ACA partial or ACA+ algorithms). 506 * 507 * 4. advanced_compute: this member is a pointer to a function 508 * 509 * void compute(struct hmat_block_compute_context_t*) 510 * 511 * This case is similar to the previous one when block is full or compressed 512 * with SVD or ACA full algorithms. But with ACA partial or ACA+ algorithms, 513 * assembly is done by material (called stratum), and interactions are summed 514 * up. The prepare function must set hmat_block_info_t.number_of_strata, and 515 * a loop on strata is performed. By convention, if hmat_block_info_t.stratum 516 * is -1, this callback must sum up interaction for all strata. Otherwise, it 517 * must compute only the interactions of the given stratum. 518 * 519 * Only one of advanced_compute, block_compute, simple_compute and 520 * assembly function pointers must be non null. 521 * 522 */ 523 /** First scenario */ 524 void * assembly; 525 /** Second scenario */ 526 hmat_interaction_func_t simple_compute; 527 /** Third scenario */ 528 hmat_compute_func_t block_compute; 529 /** Fourth scenario */ 530 void (*advanced_compute)(struct hmat_block_compute_context_t*); 531 532 /** 533 * The user context used in all scenarii but the first one. The default is NULL. 534 */ 535 void* user_context; 536 537 /** Auxiliary method used by third and fourth scenarii 538 * This member is a pointer to function 539 * void prepare(int row_start, int row_count, int col_start, int col_count, 540 * int *row_hmat2client, int *row_client2hmat, 541 * int *col_hmat2client, int *col_client2hmat, 542 * void *context, 543 * hmat_block_info_t * block_info) 544 * Eight first arguments are input arguments 545 * Block_info is an output argument, it is initialized by hmat, and this 546 * callback must fill it up; see comments in hmat_block_info_t. 547 * Context should be passed to hmat_block_info_t.user_data. 548 */ 549 hmat_prepare_func_t prepare; 550 551 /** Compression algorithm (created by calling an hmat_create_compression_* function, or NULL for no compression) **/ 552 const hmat_compression_algorithm_t * compression; 553 554 /** Copy left lower values to the upper right of the matrix */ 555 int lower_symmetric; 556 /** The type of factorization to do after this assembling. The default is hmat_factorization_none. */ 557 hmat_factorization_t factorization; 558 /** NULL disable progress display. The default is to use the hmat progress internal implementation. */ 559 hmat_progress_t * progress; 560 } hmat_assemble_context_t; 561 562 /** Init a hmat_assemble_context_t with default values */ 563 HMAT_API void hmat_assemble_context_init(hmat_assemble_context_t * context); 564 565 typedef struct { 566 /** The type of factorization to do after this assembling. The default is hmat_factorization_lu. */ 567 hmat_factorization_t factorization; 568 /** NULL disable progress display. The default is to use the hmat progress internal implementation. */ 569 hmat_progress_t * progress; 570 } hmat_factorization_context_t; 571 572 /** Init a hmat_factorization_context_t with default values */ 573 HMAT_API void hmat_factorization_context_init(hmat_factorization_context_t * context); 574 575 /** Context for the get_values and get_block function */ 576 struct hmat_get_values_context_t { 577 /** The matrix from witch to get values */ 578 hmat_matrix_t* matrix; 579 /** Where output values are stored. Must be allocated by the caller */ 580 void * values; 581 /** 582 * @brief min row and col indices to get 583 */ 584 int row_offset, col_offset; 585 /** number of rows and cols to get */ 586 int row_size, col_size; 587 /** 588 * @brief Indirection array for numbering. 589 * Those are pointer to internal data and must not be freed by the caller. 590 */ 591 int * row_indices, * col_indices; 592 /** 593 * @brief If true renumber rows when getting blocks. 594 * This is only supported if row_offset = 0 and row_size = matrix row size. 595 */ 596 int renumber_rows:1; 597 }; 598 599 /* Opaque pointer */ 600 typedef struct 601 { 602 hmat_value_t value_type; 603 /** For internal use only */ 604 void * internal; 605 } hmat_procedure_t; 606 607 /* Delete a procedure */ 608 HMAT_API void hmat_delete_procedure(hmat_procedure_t* proc); 609 610 /* Opaque pointer */ 611 typedef struct 612 { 613 hmat_value_t value_type; 614 /** For internal use only */ 615 const void * internal; 616 } hmat_leaf_procedure_t; 617 618 /* Delete a leaf procedure */ 619 HMAT_API void hmat_delete_leaf_procedure(hmat_leaf_procedure_t* proc); 620 621 typedef struct 622 { 623 /*! Create an empty (not assembled) HMatrix from 2 \a ClusterTree instances 624 and an admissibility condition. 625 626 The HMatrix is built on a row and a column tree, that are provided as 627 arguments. 628 629 \param stype the scalar type 630 \param rows_tree a ClusterTree as returned by \a hmat_create_cluster_tree(). 631 \param cols_tree a ClusterTree as returned by \a hmat_create_cluster_tree(). 632 \param cond an admissibility condition, as returned by \a hmat_create_admissibility_standard(). 633 \param lower_symmetric 1 if the matrix is lower symmetric, 0 otherwise 634 \return an opaque pointer to an HMatrix, or NULL in case of error. 635 */ 636 hmat_matrix_t* (*create_empty_hmatrix_admissibility)(const hmat_cluster_tree_t* rows_tree, const hmat_cluster_tree_t* cols_tree, 637 int lower_symmetric, hmat_admissibility_t* cond); 638 639 /*! Declare that the given HMatrix owns its cluster trees, which means that they will be 640 automatically deallocated when HMatrix is destroyed. 641 642 \param hmatrix HMatrix 643 \param owns_row if 1, declare that this HMatrix owns its row cluster tree 644 \param owns_col if 1, declare that this HMatrix owns its column cluster tree 645 */ 646 void (*own_cluster_trees)(hmat_matrix_t* hmatrix, int owns_row, int owns_col); 647 648 /*! Set threshold for low-rank recompressions. 649 650 \param hmatrix HMatrix 651 \param epsilon Threshold for low-rank recompressions 652 */ 653 void (*set_low_rank_epsilon)(hmat_matrix_t* hmatrix, double epsilon); 654 655 /*! Assemble a HMatrix */ 656 int (*assemble_generic)(hmat_matrix_t* matrix, hmat_assemble_context_t * context); 657 658 /*! \brief Return a copy of a HMatrix. 659 660 \param from_hmat the source matrix 661 \return a copy of the matrix, or NULL 662 */ 663 hmat_matrix_t* (*copy)(hmat_matrix_t* hmatrix); 664 665 /** Return a null matrix with the same structure */ 666 hmat_matrix_t* (*copy_struct)(hmat_matrix_t* hmatrix); 667 668 /*! \brief Compute the norm of a HMatrix. 669 670 \param hmatrix the matrix of which to compute the norm 671 \return the norm 672 */ 673 double (*norm)(hmat_matrix_t* hmatrix); 674 /*! \brief Inverse a HMatrix in place. 675 676 \param hmatrix the matrix to inverse 677 \return 0 678 */ 679 int (*inverse)(hmat_matrix_t* hmatrix); 680 681 /*! Factorize a HMatrix */ 682 int (*factorize_generic)(hmat_matrix_t* matrix, hmat_factorization_context_t * context); 683 684 /*! \brief Destroy a HMatrix. 685 686 \param hmatrix the matrix to destroy 687 \return 0 688 */ 689 int (*destroy)(hmat_matrix_t* hmatrix); 690 691 /*! \brief Get one of the children of an HMatrix. 692 693 \param hmatrix the main matrix 694 \param i the row coordinate of the child 695 \param j the column coordinate of the child 696 \return 0 697 */ 698 hmat_matrix_t * (*get_child)( hmat_matrix_t *hmatrix, int i, int j ); 699 700 /*! \brief Destroy a child HMatrix. 701 702 \param hmatrix the child matrix to destroy 703 \return 0 704 */ 705 int (*destroy_child)(hmat_matrix_t* hmatrix); 706 707 /*! \brief Solve A X = B, with X overwriting B. 708 709 In this function, B is a H-matrix 710 \param hmatrix 711 \param hmatrixb 712 \return 0 for success 713 */ 714 int (*solve_mat)(hmat_matrix_t* hmatrix, hmat_matrix_t* hmatrixB); 715 /*! \brief Solve A x = b, with x overwriting b. 716 717 In this function, b is a multi-column vector, with nrhs RHS. 718 719 \param hmatrix 720 \param b 721 \param nrhs 722 \return 0 for success 723 */ 724 int (*solve_systems)(hmat_matrix_t* hmatrix, void* b, int nrhs); 725 /*! \brief Solve A x = b, with x overwriting b. 726 727 In this function, b is a multi-column vector, with nrhs RHS. 728 729 Functions vector_reorder and vector_restore must be called before 730 and after this function to transform original numbering of b 731 into/from hmat internal numbering. 732 733 \param hmatrix 734 \param b 735 \param nrhs 736 \return 0 for success 737 */ 738 int (*solve_dense)(hmat_matrix_t* hmatrix, void* b, int nrhs); 739 /*! \brief Transpose an HMatrix in place. 740 741 \return 0 for success. 742 */ 743 int (*transpose)(hmat_matrix_t *hmatrix); 744 /*! \brief A <- alpha * A 745 746 \param alpha 747 \param hmatrix 748 */ 749 int (*scale)(void * alpha, hmat_matrix_t *hmatrix); 750 /*! \brief Recompress all Rk matrices to their respective epsilon_ values. 751 752 This function is useful only after epsilon_ values have been modified. 753 754 \param hmatrix 755 */ 756 int (*truncate)(hmat_matrix_t *hmatrix); 757 /*! \brief Permute values in a dense array 758 759 \param vec_b values 760 \param rows_ct cluster tree for rows; may be NULL, in which case rows argument must be set to define 761 \param rows when rows_ct is NULL, this argument contains the number of rows, it is required to know data shape. 762 When rows_ct is not NULL, this argument is unused. 763 \param cols_ct cluster tree for cols; may be NULL, in which case cols argument must be set to define 764 \param cols when cols_ct is NULL, this argument contains the number of columns, it is required to know data shape. 765 When cols_ct is not NULL, this argument is unused. 766 */ 767 int (*vector_reorder)(void* vec_b, const hmat_cluster_tree_t *rols_ct, int rows, const hmat_cluster_tree_t *cols_ct, int cols); 768 /*! \brief Renumber values back to original numbering 769 770 \param vec_b values 771 \param rows_ct cluster tree for rows; may be NULL, in which case rows argument must be set to define 772 \param rows when rows_ct is NULL, this argument contains the number of rows, it is required to know data shape. 773 When rows_ct is not NULL, this argument is unused. 774 \param cols_ct cluster tree for cols; may be NULL, in which case cols argument must be set to define 775 \param cols when cols_ct is NULL, this argument contains the number of columns, it is required to know data shape. 776 When cols_ct is not NULL, this argument is unused. 777 */ 778 int (*vector_restore)(void* vec_b, const hmat_cluster_tree_t *rols_ct, int rows, const hmat_cluster_tree_t *cols_ct, int cols); 779 /*! \brief C <- alpha * A * B + beta * C 780 781 \param trans_a 'N' or 'T' 782 \param trans_b 'N' or 'T' 783 \param alpha 784 \param hmatrix 785 \param hmatrix_b 786 \param beta 787 \param hmatrix_c 788 \return 0 for success 789 */ 790 int (*gemm)(char trans_a, char trans_b, void* alpha, hmat_matrix_t* hmatrix, 791 hmat_matrix_t* hmatrix_b, void* beta, hmat_matrix_t* hmatrix_c); 792 793 /*! \brief y := y + a * x */ 794 int (*axpy)(void* a, hmat_matrix_t* x, hmat_matrix_t* hmatrix_y); 795 796 /*! \brief solve \alpha A * x */ 797 int (*trsm)( char side, char uplo, char transa, char diag, int m, int n, 798 void *alpha, hmat_matrix_t *A, int is_b_hmat, void *B ); 799 800 /*! @deprecated \brief c <- alpha * A * b + beta * c 801 802 \param trans_a 'N' or 'T' 803 \param alpha 804 \param hmatrix 805 \param vec_b 806 \param beta 807 \param vec_c 808 \param nrhs 809 \return 0 for success 810 */ 811 int (*gemv)(char trans_a, void* alpha, hmat_matrix_t* hmatrix, void* vec_b, 812 void* beta, void* vec_c, int nrhs); 813 /*! @deprecated \brief Same as gemv, but without renumbering on vec_b and vec_c 814 */ 815 int (*gemm_scalar)(char trans_a, void* alpha, hmat_matrix_t* hmatrix, void* vec_b, 816 void* beta, void* vec_c, int nrhs); 817 /*! @deprecated \brief C <- alpha * A * B + beta * C 818 819 820 In this version, a, c: FullMatrix, b: HMatrix. 821 822 \param trans_a 823 \param trans_b 824 \param mc number of rows of C 825 \param nc number of columns of C 826 \param c 827 \param alpha 828 \param a 829 \param hmat_b 830 \param beta 831 \return 0 for success 832 */ 833 int (*full_gemm)(char trans_a, char trans_b, int mc, int nc, void* c, 834 void* alpha, void* a, hmat_matrix_t* hmat_b, void* beta); 835 /*! \brief Y <- alpha * op(B) * op(X) + beta * Y if side is 'L' 836 or Y <- alpha * op(X) * op(B) + beta * Y if side is 'R' 837 838 Functions vector_reorder and vector_restore must be called before 839 and after this function to transform original numbering of X and Y 840 into/from hmat internal numbering. 841 842 \param trans_b 'N', 'T' or 'C' 843 \param trans_x 'N', 'T' or 'C' 844 \param side 'L' or 'R' 845 \param alpha 846 \param hmatrix 847 \param vec_x 848 \param beta 849 \param vec_y 850 \param nrhs 851 \return 0 for success 852 */ 853 int (*gemm_dense)(char trans_b, char trans_x, char side, void* alpha, hmat_matrix_t* holder, 854 void* vec_x, void* beta, void* vec_y, int nrhs); 855 /*! \brief A <- A + alpha Id 856 857 \param hmatrix 858 \param alpha 859 \return 0 for success 860 */ 861 int (*add_identity)(hmat_matrix_t* hmatrix, void *alpha); 862 /*! \brief Initialize library 863 */ 864 int (*init)(void); 865 866 /*! \brief Do the cleanup 867 */ 868 int (*finalize)(void); 869 870 /*! \brief Get current informations 871 \param hmatrix A hmatrix 872 \param info A structure to fill with current informations 873 */ 874 int (*get_info)(hmat_matrix_t *hmatrix, hmat_info_t* info); 875 876 /*! \brief Dump json & postscript informations about matrix 877 \param hmatrix A hmatrix 878 \param prefix A string to prefix files output */ 879 int (*dump_info)(hmat_matrix_t *hmatrix, char* prefix); 880 881 /** 882 * @brief Get cluster trees 883 * \param hmatrix A hmatrix 884 * \param rows if not NULL, will contain a pointer to rows cluster tree 885 * \param cols if not NULL, will contain a pointer to cols cluster tree 886 */ 887 int (*get_cluster_trees)(hmat_matrix_t* hmatrix, const hmat_cluster_tree_t ** rows, const hmat_cluster_tree_t ** cols); 888 889 /** 890 * @brief Replace the cluster tree in a hmatrix 891 * The provided cluster trees must be compatible with the structure of 892 * the matrix. 893 */ 894 int (*set_cluster_trees)(hmat_matrix_t* hmatrix, const hmat_cluster_tree_t * rows, const hmat_cluster_tree_t * cols); 895 896 /** 897 * @brief Extract matrix diagonal 898 * \param hmatrix A hmatrix 899 * \param diag allocated memory area in which diagonal values are written 900 * \param size memory size 901 */ 902 int (*extract_diagonal)(hmat_matrix_t* holder, void* diag, int size); 903 904 /** 905 * @brief Solve system op(L)*X=B 906 \warning There is no check to ensure that matrix has been factorized. 907 * \param hmatrix A hmatrix 908 * \param transpose if different from 0, transposed matrix is used 909 * \param b right-hand sides, overwritten by solution X at exit 910 * \param nrhs number of right-hand sides 911 */ 912 int (*solve_lower_triangular)(hmat_matrix_t* hmatrix, int transpose, void* b, int nrhs); 913 914 /** 915 * @brief Solve system op(L)*X=B 916 \warning There is no check to ensure that matrix has been factorized. 917 * \param hmatrix A hmatrix 918 * \param transpose if different from 0, transposed matrix is used 919 * \param b right-hand sides, overwritten by solution X at exit 920 * \param nrhs number of right-hand sides 921 */ 922 int (*solve_lower_triangular_dense)(hmat_matrix_t* hmatrix, int transpose, void* b, int nrhs); 923 924 /** 925 * @brief Extract and uncompress a block of the matrix. 926 * After calling this function row_indices and col_indices will contains the 927 * indices of the extract block. 928 * @see struct hmat_get_values_context_t 929 */ 930 int (*get_block)(struct hmat_get_values_context_t * ctx); 931 932 /** 933 * @brief Extract a set of values from the matrix 934 * row_indices and col_indices must contains the rows ans columns to get. 935 * row_offset, col_offset and renumber_rows are ignored 936 * @see struct hmat_get_values_context_t 937 */ 938 int (*get_values)(struct hmat_get_values_context_t * ctx); 939 940 /** 941 * @brief Apply a procedure to all nodes of a matrix 942 * \param hmatrix A hmatrix 943 * \param proc 944 */ 945 int (*walk)(hmat_matrix_t* hmatrix, hmat_procedure_t* proc); 946 947 /** 948 * @brief Apply a procedure to all leaves of a matrix 949 * \param hmatrix A hmatrix 950 * \param proc 951 */ 952 int (*apply_on_leaf)(hmat_matrix_t* hmatrix, const hmat_leaf_procedure_t* proc); 953 954 hmat_value_t value_type; 955 956 /** For internal use only */ 957 void * internal; 958 959 hmat_matrix_t * (*read_struct)(hmat_iostream readfunc, void * user_data); 960 void (*write_struct)(hmat_matrix_t* matrix, hmat_iostream writefunc, void * user_data); 961 void (*read_data)(hmat_matrix_t* matrix, hmat_iostream readfunc, void * user_data); 962 void (*write_data)(hmat_matrix_t* matrix, hmat_iostream writefunc, void * user_data); 963 964 /** 965 * @brief Set the progress bar associated to a matrix. 966 * 967 * Note that progress bar is also set by assemble_generic and factorize_generic. 968 * @param matrix the matrix whose one want to change the progressbar 969 * @param progress the new progress bar implementation. NULL disable progress 970 * reporting. 971 */ 972 void (*set_progressbar)(hmat_matrix_t * matrix, hmat_progress_t * progress); 973 } hmat_interface_t; 974 975 HMAT_API void hmat_init_default_interface(hmat_interface_t * i, hmat_value_t type); 976 977 typedef struct 978 { 979 /*! \brief svd compression if max(rows->n, cols->n) < compressionMinLeafSize.*/ 980 int compressionMinLeafSize; 981 /*! \brief Tolerance for coarsening */ 982 double coarseningEpsilon; 983 /*! \brief Maximum size of a leaf in a ClusterTree (and of a non-admissible block in an HMatrix) */ 984 int maxLeafSize; 985 /*! \brief Coarsen the matrix structure after assembly. */ 986 int coarsening; 987 /*! \brief Validate the detection of null rows and columns */ 988 int validateNullRowCol; 989 /*! \brief Validate the rk-matrices after compression */ 990 int validateCompression; 991 /*! \brief Dump trace at the end of the algorithms (depends on the runtime) */ 992 int dumpTrace; 993 /*! \brief For blocks above error threshold, re-run the compression algorithm */ 994 int validationReRun; 995 /*! \brief For blocks above error threshold, dump the faulty block to disk */ 996 int validationDump; 997 /*! \brief Error threshold for the compression validation */ 998 double validationErrorThreshold; 999 } hmat_settings_t; 1000 1001 /*! \brief Get current settings 1002 \param settings A structure to fill with current settings 1003 */ 1004 HMAT_API void hmat_get_parameters(hmat_settings_t * settings); 1005 /*! \brief Set current settings 1006 1007 \param structure containing parameters 1008 \return 1 on failure, 0 otherwise. 1009 */ 1010 HMAT_API int hmat_set_parameters(hmat_settings_t*); 1011 1012 /*! 1013 * \brief hmat_get_version 1014 * \return The version of this library 1015 */ 1016 HMAT_API const char * hmat_get_version(void); 1017 /*! 1018 * \brief hmat_get_build_date 1019 * \return The build date and time of this library 1020 */ 1021 HMAT_API void hmat_get_build_date(const char**, const char**); 1022 1023 /*! 1024 \brief hmat_tracing_dump Dumps the trace info in the given filename 1025 1026 The file is in json format. Hmat library must be compiled with -DHAVE_CONTEXT for this to work. 1027 \param filename the name of the output json file 1028 */ 1029 HMAT_API void hmat_tracing_dump(char *filename) ; 1030 1031 /** \brief Set the function used to get the worker index. 1032 1033 The function f() must return the worker Id (between 0 and nbWorkers-1) or -1 in a sequential section. 1034 This function is used within hmat-oss for timers and traces. It must be set if hmat-oss is called by 1035 multiple threads simultaneously. 1036 */ 1037 HMAT_API void hmat_set_worker_index_function(int (*f)(void)); 1038 1039 #ifdef __cplusplus 1040 } 1041 #endif 1042 #endif /* _HMAT_H */ 1043