Skip to content

Commit

Permalink
type flexibility
Browse files Browse the repository at this point in the history
  • Loading branch information
@cbouilla
cbouilla committed Sep 21, 2023
1 parent db4ea4c commit c7d0951
Show file tree
Hide file tree
Showing 16 changed files with 265 additions and 99 deletions.
9 changes: 7 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,10 @@ This enables several of useful computations on sparse matrices modulo _p_:
* Basis of the kernel
* Reduced Row Echelon form
* Rank
* Solution of linear systems (WIP --- performance not as good as the rest)

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

SpaSM works with all prime moduli in the range [3; 189,812,507]. While it would be possible to work with _p = 2_ or with any 32-bit prime _p_, this would require tweaks to the code (and is not tested).

Expand All @@ -24,7 +27,9 @@ In addition, SpaSM contains code to compute the Dulmage-Mendelson decomposition
The following algorithms algorithms are used in SpaSM:
* [Gilbert-Peierls sparse triangular solving with sparse right-hand side](https://doi.org/10.1137/0909058)
* [Faugère-Lachartre pivot selection](http://www-almasty.lip6.fr/~bouillaguet/pub/CASC16.pdf)
* [Improved greedy pivot selection](http://www-almasty.lip6.fr/~bouillaguet/pub/PASCO17.pdf).
* [Improved greedy pivot selection](http://www-almasty.lip6.fr/~bouillaguet/pub/PASCO17.pdf)
* [Dense linear algebra mod _p_](https://hal.science/hal-00018223/)
* [Linear-size rank certificates](https://prism.ucalgary.ca/server/api/core/bitstreams/b00bb76d-12bf-41c2-9fed-88cd774d3b29/content)

Initial versions of SpaSM were heavily influence by
[Tim Davis](http://faculty.cse.tamu.edu/davis/)'s [CSparse](http://faculty.cse.tamu.edu/davis/publications_files/CSparse.zip).
Expand Down
20 changes: 14 additions & 6 deletions src/spasm.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,8 @@ typedef struct {
spasm_ZZp *y; /* size r */
} spasm_rank_certificate;

typedef enum {SPASM_DOUBLE, SPASM_FLOAT, SPASM_I64} spasm_datatype;

#define SPASM_IDENTITY_PERMUTATION NULL
#define SPASM_IGNORE NULL
#define SPASM_IGNORE_VALUES 0
Expand Down Expand Up @@ -156,7 +158,7 @@ int spasm_get_thread_num();
static inline i64 spasm_get_prime(const spasm *A) { return A->field->p; }

/* spasm_triplet.c */
void spasm_add_entry(spasm_triplet * T, int i, int j, spasm_ZZp x);
void spasm_add_entry(spasm_triplet *T, int i, int j, i64 x);
void spasm_triplet_transpose(spasm_triplet * T);
spasm *spasm_compress(const spasm_triplet * T);

Expand Down Expand Up @@ -201,8 +203,10 @@ int spasm_sparse_triangular_solve(const spasm *U, const spasm *B, int k, int *xj
spasm *spasm_schur(const spasm *A, const int *p, int n, const spasm_lu *fact,
double est_density, spasm_triplet *L, const int *p_in, int *p_out);
double spasm_schur_estimate_density(const spasm * A, const int *p, int n, const spasm *U, const int *qinv, int R);
int spasm_schur_dense(const spasm *A, const int *p, int k, const spasm *U, const int *qinv, double *S, int *q);
void spasm_schur_dense_randomized(const spasm *A, const int *p, int n, const spasm *U, const int *qinv, double *S, int *q, int N, int w);
int spasm_schur_dense(const spasm *A, const int *p, int n, const spasm *U, const int *qinv,
void *S, spasm_datatype datatype, int *q);
void spasm_schur_dense_randomized(const spasm *A, const int *p, int n, const spasm *U, const int *qinv,
void *S, spasm_datatype datatype, int *q, int N, int w);

/* spasm_pivots.c */
int spasm_pivots_extract_structural(const spasm *A, const int *p_in, spasm_lu *fact, int *p, struct echelonize_opts *opts);
Expand All @@ -221,9 +225,13 @@ spasm_dm *spasm_dulmage_mendelsohn(const spasm *A);
spasm_dm *spasm_strongly_connected_components(const spasm *A);

/* spasm_ffpack.cpp */
int spasm_ffpack_rref_double(i64 prime, int n, int m, double *A, int ldA, size_t *qinv);
int spasm_ffpack_rref_float(i64 prime, int n, int m, float *A, int ldA, size_t *qinv);
int spasm_ffpack_rref_i64(i64 prime, int n, int m, i64 *A, int ldA, size_t *qinv);
int spasm_ffpack_rref(i64 prime, int n, int m, void *A, int ldA, spasm_datatype datatype, size_t *qinv);
spasm_ZZp spasm_datatype_read(const void *A, size_t i, spasm_datatype datatype);
void spasm_datatype_write(void *A, size_t i, spasm_datatype datatype, spasm_ZZp value);
size_t spasm_datatype_size(spasm_datatype datatype);
spasm_datatype spasm_datatype_choose(i64 prime);
const char * spasm_datatype_name(spasm_datatype datatype);


/* spasm_echelonize */
void spasm_echelonize_init_opts(struct echelonize_opts *opts);
Expand Down
53 changes: 27 additions & 26 deletions src/spasm_echelonize.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,6 @@
* echelonize_XXX(spasm *A, const int *p, int n, spasm *U, int *Uqinv, struct echelonize_opts *opts);
*
* passing NULL as the options leads to default choices.
* U must be preallocated with sufficiently many rows (e.g. with A->n rows to be ready for the worst case).
*
* presently, this does not provide a way to store L; we will do something about this later on.
*/


Expand Down Expand Up @@ -48,16 +45,17 @@ bool spasm_echelonize_test_completion(const spasm *A, const int *p, int n, spasm
if (n == 0 || spasm_nnz(A) == 0)
return 1;
int m = A->m;
int Sm = m - U->n;
i64 Sm = m - U->n;
i64 prime = spasm_get_prime(A);
int Sn = ceil(128 / log2(prime));
double *S = spasm_malloc(Sn * Sm * sizeof(*S));
spasm_datatype datatype = spasm_datatype_choose(prime);
i64 Sn = ceil(128 / log2(prime));
void *S = spasm_malloc(Sn * Sm * spasm_datatype_size(datatype));
int *q = spasm_malloc(Sm * sizeof(*q));
size_t *Sp = spasm_malloc(Sm * sizeof(*Sp)); /* for FFPACK */
fprintf(stderr, "[echelonize/completion] Testing completion with %d random linear combinations (rank %d)\n", Sn, U->n);
fprintf(stderr, "[echelonize/completion] Testing completion with %" PRId64" random linear combinations (rank %d)\n", Sn, U->n);
fflush(stderr);
spasm_schur_dense_randomized(A, p, n, U, Uqinv, S, q, Sn, 0);
int rr = spasm_ffpack_rref_double(prime, Sn, Sm, S, Sm, Sp);
spasm_schur_dense_randomized(A, p, n, U, Uqinv, S, datatype, q, Sn, 0);
int rr = spasm_ffpack_rref(prime, Sn, Sm, S, Sm, datatype, Sp);
free(S);
free(Sp);
free(q);
Expand Down Expand Up @@ -201,8 +199,8 @@ void spasm_echelonize_GPLU(const spasm *A, const int *p, int n, spasm_lu *fact,
free(xj);
}


static void dense_update_U(spasm *U, int rr, int Sm, const double *S, const size_t *Sqinv, const int *q, int *Uqinv)
static void dense_update_U(spasm *U, int rr, int Sm, const void *S, spasm_datatype datatype,
const size_t *Sqinv, const int *q, int *Uqinv)
{
i64 extra_nnz = ((i64) (1 + Sm - rr)) * rr; /* maximum size increase */
i64 unz = spasm_nnz(U);
Expand All @@ -219,10 +217,11 @@ static void dense_update_U(spasm *U, int rr, int Sm, const double *S, const size
Uqinv[q[j]] = U->n;
for (i64 k = rr; k < Sm; k++) {
i64 j = Sqinv[k];
if (S[i * Sm + k] == 0)
spasm_ZZp x = spasm_datatype_read(S, i * Sm + k, datatype);
if (x == 0)
continue; /* don't store zero */
Uj[unz] = q[j];
Ux[unz] = S[i * Sm + k]; // reduce?
Ux[unz] = x; // reduce?
unz += 1;
}
U->n += 1;
Expand All @@ -237,16 +236,17 @@ static void spasm_echelonize_dense_lowrank(const spasm *A, const int *p, int n,
int m = A->m;
int Sm = m - U->n;
i64 prime = spasm_get_prime(A);
spasm_datatype datatype = spasm_datatype_choose(prime);

i64 size_S = (i64) opts->dense_block_size * (i64) Sm * sizeof(double);
double *S = spasm_malloc(size_S);
i64 size_S = (i64) opts->dense_block_size * (i64) Sm * spasm_datatype_size(datatype);
void *S = spasm_malloc(size_S);
int *q = spasm_malloc(Sm * sizeof(*q));
size_t *Sp = spasm_malloc(Sm * sizeof(*Sp)); /* for FFPACK */
double start = spasm_wtime();
int old_un = U->n;
int round = 0;
fprintf(stderr, "[echelonize/dense/low-rank] processing dense schur complement of dimension %d x %d; block size=%d\n",
n, Sm, opts->dense_block_size);
fprintf(stderr, "[echelonize/dense/low-rank] processing dense schur complement of dimension %d x %d; block size=%d, type %s\n",
n, Sm, opts->dense_block_size, spasm_datatype_name(datatype));

/*
* stupid algorithm to decide a starting weight:
Expand All @@ -268,8 +268,8 @@ static void spasm_echelonize_dense_lowrank(const spasm *A, const int *p, int n,
break;
fprintf(stderr, "[echelonize/dense/low-rank] Round %d. Weight %d. Processing chunk (%d x %d), |U| = %"PRId64"\n",
round, w, Sn, Sm, spasm_nnz(U));
spasm_schur_dense_randomized(A, p, n, U, Uqinv, S, q, Sn, w);
int rr = spasm_ffpack_rref_double(prime, Sn, Sm, S, Sm, Sp);
spasm_schur_dense_randomized(A, p, n, U, Uqinv, S, datatype, q, Sn, w);
int rr = spasm_ffpack_rref(prime, Sn, Sm, S, Sm, datatype, Sp);

if (rr == 0) {
if (spasm_echelonize_test_completion(A, p, n, U, Uqinv))
Expand All @@ -282,7 +282,7 @@ static void spasm_echelonize_dense_lowrank(const spasm *A, const int *p, int n,
w *= 2;
fprintf(stderr, "[echelonize/dense/low-rank] Not enough pivots, increasing weight to %d\n", w);
}
dense_update_U(U, rr, Sm, S, Sp, q, Uqinv);
dense_update_U(U, rr, Sm, S, datatype, Sp, q, Uqinv);
n -= rr;
Sm -= rr;
round += 1;
Expand All @@ -304,16 +304,17 @@ static void spasm_echelonize_dense(const spasm *A, const int *p, int n, spasm *U
int m = A->m;
int Sm = m - U->n;
i64 prime = spasm_get_prime(A);
spasm_datatype datatype = spasm_datatype_choose(prime);

double *S = spasm_malloc(opts->dense_block_size * Sm * sizeof(*S));
void *S = spasm_malloc(opts->dense_block_size * Sm * spasm_datatype_size(datatype));
int *q = spasm_malloc(Sm * sizeof(*q));
size_t *Sp = spasm_malloc(Sm * sizeof(*Sp)); /* for FFPACK */
int processed = 0;
double start = spasm_wtime();
int old_un = U->n;
int round = 0;
fprintf(stderr, "[echelonize/dense] processing dense schur complement of dimension %d x %d; block size=%d\n",
n, Sm, opts->dense_block_size);
fprintf(stderr, "[echelonize/dense] processing dense schur complement of dimension %d x %d; block size=%d, type %s\n",
n, Sm, opts->dense_block_size, spasm_datatype_name(datatype));
bool lowrank_mode = 0;
int rank_ub = spasm_min(A->n - U->n, A->m - U->n);

Expand All @@ -324,11 +325,11 @@ static void spasm_echelonize_dense(const spasm *A, const int *p, int n, spasm *U
break;

fprintf(stderr, "[echelonize/dense] Round %d. processing S[%d:%d] (%d x %d)\n", round, processed, processed + Sn, Sn, Sm);
int r = spasm_schur_dense(A, p, Sn, U, Uqinv, S, q);
int rr = spasm_ffpack_rref_double(prime, r, Sm, S, Sm, Sp);
int r = spasm_schur_dense(A, p, Sn, U, Uqinv, S, datatype, q);
int rr = spasm_ffpack_rref(prime, r, Sm, S, Sm, datatype, Sp);

/* update U */
dense_update_U(U, rr, Sm, S, Sp, q, Uqinv);
dense_update_U(U, rr, Sm, S, datatype, Sp, q, Uqinv);

/* move on to the next chunk */
round += 1;
Expand Down
60 changes: 53 additions & 7 deletions src/spasm_ffpack.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,19 +43,65 @@ static inline int spasm_ffpack_rref(i64 prime, int n, int m, T *A, int ldA, size
return rank;
}

/* force instatiations */
int spasm_ffpack_rref(i64 prime, int n, int m, void *A, int ldA, spasm_datatype datatype, size_t *qinv)
{
switch (datatype) {
case SPASM_DOUBLE: return spasm_ffpack_rref<double>(prime, n, m, (double *) A, ldA, qinv);
case SPASM_FLOAT: return spasm_ffpack_rref<float>(prime, n, m, (float *) A, ldA, qinv);
case SPASM_I64: return spasm_ffpack_rref<i64>(prime, n, m, (i64 *) A, ldA, qinv);
}
assert(false);
}

/* code to make the "runtime type selection" mechanisme work */

spasm_ZZp spasm_datatype_read(const void *A, size_t i, spasm_datatype datatype)
{
switch (datatype) {
case SPASM_DOUBLE: return ((double *) A)[i];
case SPASM_FLOAT: return ((float *) A)[i];
case SPASM_I64: return ((i64 *) A)[i];
}
assert(false);
}

int spasm_ffpack_rref_double(i64 prime, int n, int m, double *A, int ldA, size_t *qinv)
void spasm_datatype_write(void *A, size_t i, spasm_datatype datatype, spasm_ZZp value)
{
return spasm_ffpack_rref(prime, n, m, A, ldA, qinv);
switch (datatype) {
case SPASM_DOUBLE: ((double *) A)[i] = value; return;
case SPASM_FLOAT: ((float *) A)[i] = value; return;
case SPASM_I64: ((i64 *) A)[i] = value; return;
}
assert(false);
}

int spasm_ffpack_rref_float(i64 prime, int n, int m, float *A, int ldA, size_t *qinv)
size_t spasm_datatype_size(spasm_datatype datatype)
{
return spasm_ffpack_rref(prime, n, m, A, ldA, qinv);
switch (datatype) {
case SPASM_DOUBLE: return sizeof(double);
case SPASM_FLOAT: return sizeof(float);
case SPASM_I64: return sizeof(i64);
}
assert(false);
}

int spasm_ffpack_rref_i64(i64 prime, int n, int m, i64 *A, int ldA, size_t *qinv)
spasm_datatype spasm_datatype_choose(i64 prime)
{
return spasm_ffpack_rref(prime, n, m, A, ldA, qinv);
// return SPASM_DOUBLE;
if (prime <= 8191)
return SPASM_FLOAT;
else if (prime <= 189812531)
return SPASM_DOUBLE;
else
return SPASM_I64;
}

const char * spasm_datatype_name(spasm_datatype datatype)
{
switch (datatype) {
case SPASM_DOUBLE: return "double";
case SPASM_FLOAT: return "float";
case SPASM_I64: return "i64";
}
assert(false);
}
2 changes: 1 addition & 1 deletion src/spasm_kernel.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ spasm * spasm_kernel(const spasm_lu *fact)

/* write the new row in K */
Kj[local_nnz] = j;
Kx[local_nnz] = prime - 1;
Kx[local_nnz] = -1;
local_nnz += 1;
for (int px = top; px < n; px++) {
int jj = xj[px];
Expand Down
Loading

0 comments on commit c7d0951

Please sign in to comment.