diff --git a/README.md b/README.md index 2e4f77a..b89277a 100644 --- a/README.md +++ b/README.md @@ -10,6 +10,7 @@ Features -------- The core of the library is a multithreaded function that computes the row echelon form of a sparse matrix modulo a word-sized prime. + This enables several of useful computations on sparse matrices modulo _p_: * Basis of the row space * Basis of the kernel @@ -18,7 +19,8 @@ This enables several of useful computations on sparse matrices modulo _p_: In addition, SpaSM is capable of computing a full PLUQ factorization, but this is slower than simple echelonization. This enables: * Solution of linear systems - * Rank certificates + * Extraction of a square submatrix of maximal rank + * Production of rank certificates SpaSM works with all odd 32-bit prime moduli. While it would be possible to work with _p = 2_, this would require non-trivial tweaks to the code. diff --git a/src/spasm.h b/src/spasm.h index af5f86d..36841b4 100644 --- a/src/spasm.h +++ b/src/spasm.h @@ -15,6 +15,10 @@ typedef int32_t i32; #include #endif +#define SPASM_VERSION "1.3" +#define SPASM_BUG_ADDRESS "" + + /* --- primary SpaSM routines and data structures --- */ // unfortunately we use "n" for #rows and "m" for #columns whereas the rest of the world (BLAS...) diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 4c5094e..87fddee 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -1,5 +1,7 @@ include_directories (../src) +add_library(spasmtools common.c) + ############### tools add_executable(bitmap bitmap.c) @@ -24,10 +26,10 @@ add_executable(dm dm.c) target_link_libraries(dm PUBLIC spasm) add_executable(rank rank.c) -target_link_libraries(rank PUBLIC spasm) +target_link_libraries(rank PUBLIC spasm spasmtools) add_executable(echelonize echelonize.c) -target_link_libraries(echelonize PUBLIC spasm) +target_link_libraries(echelonize PUBLIC spasm spasmtools) add_executable(kernel kernel.c) -target_link_libraries(kernel PUBLIC spasm) +target_link_libraries(kernel PUBLIC spasm spasmtools) diff --git a/tools/common.c b/tools/common.c new file mode 100644 index 0000000..696cd17 --- /dev/null +++ b/tools/common.c @@ -0,0 +1,134 @@ +#include +#include +#include +#include + +#include "spasm.h" +#include "common.h" + +const char *argp_program_version = SPASM_VERSION; +const char *argp_program_bug_address = SPASM_BUG_ADDRESS; + +struct argp_option input_options[] = { + {0, 0, 0, 0, "Input problem", 1 }, + {"matrix", 'm', "FILE", 0, "Read the input matrix from FILE", 1 }, + {"modulus", 'p', "P", 0, "Perform arithmetic modulo P", 1 }, + {0} +}; + +/* Parse a single option --- rank program. */ +error_t parse_input_opt(int key, char *arg, struct argp_state *state) +{ + struct input_matrix *arguments = state->input; + switch (key) { + case ARGP_KEY_INIT: /* set defaults */ + arguments->filename = NULL; + arguments->prime = 42013; + break; + case 'm': + arguments->filename = arg; + break; + case 'p': + arguments->prime = atoll(arg); + break; + default: + return ARGP_ERR_UNKNOWN; + } + return 0; +} + +struct argp input_argp = { input_options, parse_input_opt, NULL, NULL, NULL, NULL, NULL}; + + +/* The options of the echelonization code */ +enum ech_opt_key { + NO_LOW_RANK, NO_DENSE, NO_GPLU, + MAX_ITER, DENSE_THR, MIN_PIV_RATIO, + DENSE_BLKSZ, MIN_RANK_RATIO, MAX_ASPECT_RATIO +}; + +struct argp_option echelonize_options[] = { + {0, 0, 0, 0, "Echelonization sub-algorithms", -2}, + {"no-low-rank-mode", NO_LOW_RANK, 0, 0, "Disable the (dense) low-rank mode", -2 }, + {"no-dense-mode", NO_DENSE, 0, 0, "Don't use FFPACK", -2 }, + {"no-GPLU", NO_GPLU, 0, 0, "Don't use GPLU", -2 }, + + {0, 0, 0, 0, "Main echelonization options", -3 }, + {"max-iterations", MAX_ITER, "N", 0, "Compute at most N sparse Schur complements ", -3}, + {"dense-threshold", DENSE_THR, "D", 0, "Use FFPACK when the density is greater than D", -3}, + {"min-pivot-proportion", MIN_PIV_RATIO, "P", 0, "Stop if the proportion of pivots found is less than P", -3}, + + {0, 0, 0, 0, "Dense code options", -4}, + {"dense-block-size", DENSE_BLKSZ, "N", 0, "Use dense matrices of at most N rows", -4}, + {"min-rank-ratio", MIN_RANK_RATIO, "R", 0, "Use low-rank mode if k rows have rank <= k * X", -4}, + {"max-aspect-ratio", MAX_ASPECT_RATIO, "R", 0, "Use low-rank mode if #rows / #columns >= R", -4}, + + { 0 } +}; + +/* Parse a single option --- rank program. */ +error_t parse_ech_opt(int key, char *arg, struct argp_state *state) +{ + struct echelonize_opts *opts = state->input; + switch (key) { + case ARGP_KEY_INIT: /* set defaults */ + spasm_echelonize_init_opts(opts); + break; + case NO_LOW_RANK: + opts->enable_tall_and_skinny = 0; + break; + case NO_DENSE: + opts->enable_dense = 0; + break; + case NO_GPLU: + opts->enable_GPLU = 0; + break; + case MAX_ITER: + opts->max_round = atoi(arg); + break; + case DENSE_THR: + opts->sparsity_threshold = atof(arg); + break; + case MIN_PIV_RATIO: + opts->min_pivot_proportion = atof(arg); + break; + case DENSE_BLKSZ: + opts->dense_block_size = atoi(arg); + break; + case MIN_RANK_RATIO: + opts->low_rank_ratio = atof(arg); + break; + case MAX_ASPECT_RATIO: + opts->tall_and_skinny_ratio = atof(arg); + break; + default: + return ARGP_ERR_UNKNOWN; + } + return 0; +} + +struct argp echelonize_argp = { echelonize_options, parse_ech_opt, NULL, NULL, NULL, NULL, NULL}; + + +spasm_triplet * load_input_matrix(struct input_matrix *in) +{ + if (in->filename == NULL) + return spasm_load_sms(stdin, in->prime); /* load from stdin */ + FILE *f = fopen(in->filename, "r"); + if (f == NULL) + err(1, "Cannot open %s", in->filename); + spasm_triplet *T = spasm_load_sms(f, in->prime); + fclose(f); + return T; +} + + +FILE * open_output(const char *filename) +{ + if (filename == NULL) + return stdout; + FILE *f = fopen(filename, "w"); + if (f == NULL) + err(1, "Cannot open %s", filename); + return f; +} diff --git a/tools/common.h b/tools/common.h new file mode 100644 index 0000000..b93538b --- /dev/null +++ b/tools/common.h @@ -0,0 +1,20 @@ +#ifndef _SPASM_COMMON_H +#define _SPASM_COMMON_H + +#include +#include + +struct input_matrix { + char *filename; + i64 prime; +}; + +extern struct argp echelonize_argp; +extern struct argp input_argp; +extern const char *argp_program_version; +extern const char *argp_program_bug_address; + +spasm_triplet * load_input_matrix(struct input_matrix *in); +FILE * open_output(const char *filename); + +#endif \ No newline at end of file diff --git a/tools/kernel.c b/tools/kernel.c index 8b7ca7d..9cd8cbb 100644 --- a/tools/kernel.c +++ b/tools/kernel.c @@ -5,57 +5,73 @@ #include #include "spasm.h" +#include "common.h" -i64 prime = 42013; -struct echelonize_opts opts; -bool left = 0; +/* Program documentation. */ +char doc[] = "Compute a kernel basis of a sparse matrix"; -void parse_command_line_options(int argc, char **argv) +struct cmdline_args { + /* input problem */ + struct input_matrix input; + + /* common options to all programs that call spasm_echelonize() */ + struct echelonize_opts opts; + + /* options specific to the kernel program */ + bool left; + char *output_filename; +}; + +/* The options we understand. */ +struct argp_option options[] = { + {0, 0, 0, 0, "Kernel options", 2 }, + {"left", 'l', 0, 0, "Compute the left-kernel", 2}, + {"output", 'o', "FILE", 0, "Write the kernel basis in FILE", 2 }, + { 0 } +}; + +/* Parse a single option --- rank program. */ +error_t parse_ker_opt(int key, char *arg, struct argp_state *state) { - struct option longopts[] = { - {"modulus", required_argument, NULL, 'p'}, - {"left", no_argument, NULL, 'L'}, - {"no-greedy-pivot-search", no_argument, NULL, 'g'}, - {"no-low-rank-mode", no_argument, NULL, 'l'}, - {"dense-block-size", required_argument, NULL, 'd'}, - {"low-rank-start-weight", required_argument, NULL, 'w'}, - {NULL, 0, NULL, 0} - }; - char ch; - while ((ch = getopt_long(argc, argv, "", longopts, NULL)) != -1) { - switch (ch) { - case 'p': - prime = atoll(optarg); - break; - case 'L': - left = 1; - break; - case 'g': - opts.enable_greedy_pivot_search = 0; - break; - case 'd': - opts.dense_block_size = atoi(optarg); - fprintf(stderr, "Using dense block size %d\n", opts.dense_block_size); - break; - case 'l': - opts.enable_tall_and_skinny = 0; - break; - case 'w': - opts.low_rank_start_weight = atoi(optarg); - break; - default: - errx(1, "Unknown option\n"); - } + struct cmdline_args *arguments = state->input; + switch (key) { + case 'l': + arguments->left = 1; + break; + case 'o': + arguments->output_filename = arg; + break; + case ARGP_KEY_INIT: + arguments->left = 0; + arguments->output_filename = NULL; + state->child_inputs[0] = &arguments->input; + state->child_inputs[1] = &arguments->opts; + break; + default: + return ARGP_ERR_UNKNOWN; } + return 0; } +struct argp_child children_parsers[] = { + { &input_argp, 0, 0, 0 }, + { &echelonize_argp, 0, 0, 0 }, + { 0 } +}; + +/* argp parser for the rank program */ +struct argp argp = { options, parse_ker_opt, NULL, doc, children_parsers, NULL, NULL }; + + int main(int argc, char **argv) { - spasm_echelonize_init_opts(&opts); - parse_command_line_options(argc, argv); + /* process command-line options */ + struct cmdline_args args; + argp_parse(&argp, argc, argv, 0, 0, &args); - spasm_triplet *T = spasm_load_sms(stdin, prime); - if (left) { + /* load input matrix */ + spasm_triplet *T = load_input_matrix(&args.input); + if (args.left) { fprintf(stderr, "Left-kernel, transposing\n"); spasm_triplet_transpose(T); } @@ -63,14 +79,13 @@ int main(int argc, char **argv) spasm_triplet_free(T); /* echelonize A */ - spasm_lu *fact = spasm_echelonize(A, &opts); + spasm_lu *fact = spasm_echelonize(A, &args.opts); spasm_csr_free(A); /* kernel basis */ spasm *K = spasm_kernel(fact); - spasm_save_csr(stdout, K); fprintf(stderr, "Kernel basis matrix is %d x %d with %" PRId64 " nz\n", K->n, K->m, spasm_nnz(K)); - spasm_lu_free(fact); - spasm_csr_free(K); - exit(EXIT_SUCCESS); -} + + FILE *f = open_output(args.output_filename); + spasm_save_csr(f, K); +} \ No newline at end of file diff --git a/tools/rank.c b/tools/rank.c index c3192f1..8e58030 100644 --- a/tools/rank.c +++ b/tools/rank.c @@ -1,52 +1,83 @@ #include #include #include -#include +#include #include "spasm.h" +#include "common.h" -i64 prime = 42013; -bool compute_cert = 0; -bool allow_transpose = 1; /* transpose ON by default */ -struct echelonize_opts opts; +/* Program documentation. */ +char doc[] = "Compute the rank of a sparse matrix"; -void parse_command_line_options(int argc, char **argv) +struct cmdline_args { + /* input problem */ + struct input_matrix input; + + /* common options to all programs that call spasm_echelonize() */ + struct echelonize_opts opts; + + /* options specific to the rank program */ + bool certificate; + char *cert_file; /* for the certificate */ + bool allow_transpose; /* transpose ON by default */ +}; + +/* The options we understand. */ +struct argp_option options[] = { + {0, 0, 0, 0, "Rank options", 2 }, + {"no-transpose", 't', 0, 0, "Do not transpose the input matrix", 2}, + {"certificate", 'c', 0, 0, "Output a rank certificate", 2 }, + {"output", 'o', "FILE", 0, "Write the rank certificate in FILE", 2 }, + { 0 } +}; + +/* Parse a single option --- rank program. */ +error_t parse_rank_opt(int key, char *arg, struct argp_state *state) { - struct option longopts[7] = { - {"no-transpose", no_argument, NULL, 'a'}, - {"modulus", required_argument, NULL, 'p'}, - {"certificate", no_argument, NULL, 'c'}, - {NULL, 0, NULL, 0} - }; - char ch; - while ((ch = getopt_long(argc, argv, "", longopts, NULL)) != -1) { - switch (ch) { - case 'a': - allow_transpose = 0; - break; - case 'p': - prime = atoll(optarg); - break; - case 'c': - compute_cert = 1; - break; - default: - printf("Unknown option\n"); - exit(1); - } + struct cmdline_args *arguments = state->input; + switch (key) { + case 't': + arguments->allow_transpose = 0; + break; + case 'c': + arguments->certificate = 1; + break; + case 'o': + arguments->cert_file = arg; + break; + case ARGP_KEY_INIT: + arguments->allow_transpose = 1; + arguments->certificate = 0; + arguments->cert_file = NULL; + state->child_inputs[0] = &arguments->input; + state->child_inputs[1] = &arguments->opts; + break; + default: + return ARGP_ERR_UNKNOWN; } + return 0; } +struct argp_child children_parsers[] = { + { &input_argp, 0, 0, 0 }, + { &echelonize_argp, 0, 0, 0 }, + {0} +}; + +/* argp parser for the rank program */ +struct argp argp = { options, parse_rank_opt, NULL, doc, children_parsers, NULL, NULL }; + /** Computes the rank of the input matrix using the hybrid strategy described in [PASCO'17] */ int main(int argc, char **argv) { - spasm_echelonize_init_opts(&opts); - parse_command_line_options(argc, argv); + /* process command-line options */ + struct cmdline_args args; + argp_parse(&argp, argc, argv, 0, 0, &args); - /* load matrix from standard input */ - spasm_triplet *T = spasm_load_sms(stdin, prime); - if (allow_transpose && (T->n < T->m)) { + /* load input matrix */ + spasm_triplet *T = load_input_matrix(&args.input); + if (args.allow_transpose && (T->n < T->m)) { fprintf(stderr, "[rank] transposing matrix : "); fflush(stderr); spasm_triplet_transpose(T); @@ -59,17 +90,17 @@ int main(int argc, char **argv) spasm_human_format(spasm_nnz(A), hnnz); fprintf(stderr, "start. A is %d x %d (%s nnz)\n", n, m, hnnz); - if (compute_cert) - opts.L = 1; + if (args.certificate) + args.opts.L = 1; double start_time = spasm_wtime(); - spasm_lu *fact = spasm_echelonize(A, &opts); /* NULL = default options */ + spasm_lu *fact = spasm_echelonize(A, &args.opts); /* NULL = default options */ double end_time = spasm_wtime(); fprintf(stderr, "done in %.3f s rank = %d\n", end_time - start_time, fact->U->n); spasm_rank_certificate *proof = NULL; - if (compute_cert) { + if (args.certificate) { assert(spasm_factorization_verify(A, fact, 42)); assert(spasm_factorization_verify(A, fact, 1337)); assert(spasm_factorization_verify(A, fact, 21011984));