From e432b3743d3d32f6060600af580e49d7033dcae9 Mon Sep 17 00:00:00 2001 From: Vitaly Magerya Date: Thu, 31 May 2018 17:43:56 +0200 Subject: [PATCH] Add Markowitz-ordered Gaussian elimination algorithm. This algorithm avoids the 'fill-in' problem of Gaussian elimination and significantly improves the times for solving large sparse systems. See: . --- AUTHORS | 2 +- check/check_lsolve.cpp | 55 ++++++++--------- ginac/flags.h | 8 ++- ginac/inifcns.cpp | 26 +++++++- ginac/matrix.cpp | 134 +++++++++++++++++++++++++++++++++++++++-- ginac/matrix.h | 1 + 6 files changed, 189 insertions(+), 37 deletions(-) diff --git a/AUTHORS b/AUTHORS index d6eaa874..7fb3466b 100644 --- a/AUTHORS +++ b/AUTHORS @@ -12,7 +12,7 @@ Contributors of patches ----------------------- Roberto Bagnara, Do Hoang Son, Markus Nullmeier, Pearu Peterson, Benedikt -Pluemper, Ben Sapp, Stefan Weinzierl. +Pluemper, Ben Sapp, Stefan Weinzierl, Vitaly Magerya. (Please send email if you think you were forgotten.) diff --git a/check/check_lsolve.cpp b/check/check_lsolve.cpp index 8c4f438b..247776ba 100644 --- a/check/check_lsolve.cpp +++ b/check/check_lsolve.cpp @@ -118,7 +118,6 @@ static unsigned check_inifcns_lsolve(unsigned n) } lst eqns; // equation list lst vars; // variable list - ex sol; // solution // Create a random linear system... for (unsigned i=0; i(e)) { + es.insert(e); + } else { + for (const ex &sube : e) { + insert_symbols(es, sube); + } + } +} + +static exset symbolset(const ex &e) +{ + exset s; + insert_symbols(s, e); + return s; +} + ex lsolve(const ex &eqns, const ex &symbols, unsigned options) { // solve a system of linear equations @@ -1077,8 +1095,10 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options) for (size_t r=0; r(symbols.op(c)),1); linpart -= co*symbols.op(c); sys(r,c) = co; @@ -1088,11 +1108,13 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options) } // test if system is linear and fill vars matrix + const exset sys_syms = symbolset(sys); + const exset rhs_syms = symbolset(rhs); for (size_t i=0; i colid(aug.cols()); + for (unsigned c = 0; c < aug.cols(); c++) { + colid[c] = c; + } switch(algo) { case solve_algo::gauss: aug.gauss_elimination(); @@ -1050,8 +1054,13 @@ matrix matrix::solve(const matrix & vars, aug.division_free_elimination(); break; case solve_algo::bareiss: - default: aug.fraction_free_elimination(); + break; + case solve_algo::markowitz: + colid = aug.markowitz_elimination(n); + break; + default: + throw std::invalid_argument("matrix::solve(): 'algo' is not one of the solve_algo enum"); } // assemble the solution matrix: @@ -1071,18 +1080,18 @@ matrix matrix::solve(const matrix & vars, // assign solutions for vars between fnz+1 and // last_assigned_sol-1: free parameters for (unsigned c=fnz; c= n). */ +std::vector +matrix::markowitz_elimination(unsigned n) +{ + GINAC_ASSERT(n <= col); + std::vector rowcnt(row, 0); + std::vector colcnt(col, 0); + // Normalize everything before start. We'll keep all the + // cells normalized throughout the algorithm to properly + // handle unnormal zeros. + for (unsigned r = 0; r < row; r++) { + for (unsigned c = 0; c < col; c++) { + if (!m[r*col + c].is_zero()) { + m[r*col + c] = m[r*col + c].normal(); + rowcnt[r]++; + colcnt[c]++; + } + } + } + std::vector colid(col); + for (unsigned c = 0; c < col; c++) { + colid[c] = c; + } + exvector ab(row); + for (unsigned k = 0; (k < col) && (k < row - 1); k++) { + // Find the pivot that minimizes (rowcnt[r]-1)*(colcnt[c]-1). + unsigned pivot_r = row + 1; + unsigned pivot_c = col + 1; + int pivot_m = row*col; + for (unsigned r = k; r < row; r++) { + for (unsigned c = k; c < n; c++) { + const ex &mrc = m[r*col + c]; + if (mrc.is_zero()) + continue; + GINAC_ASSERT(rowcnt[r] > 0); + GINAC_ASSERT(colcnt[c] > 0); + int measure = (rowcnt[r] - 1)*(colcnt[c] - 1); + if (measure < pivot_m) { + pivot_m = measure; + pivot_r = r; + pivot_c = c; + } + } + } + if (pivot_m == row*col) { + // The rest of the matrix is zero. + break; + } + GINAC_ASSERT(k <= pivot_r && pivot_r < row); + GINAC_ASSERT(k <= pivot_c && pivot_c < col); + // Swap the pivot into (k, k). + if (pivot_c != k) { + for (unsigned r = 0; r < row; r++) { + m[r*col + pivot_c].swap(m[r*col + k]); + } + std::swap(colid[pivot_c], colid[k]); + std::swap(colcnt[pivot_c], colcnt[k]); + } + if (pivot_r != k) { + for (unsigned c = k; c < col; c++) { + m[pivot_r*col + c].swap(m[k*col + c]); + } + std::swap(rowcnt[pivot_r], rowcnt[k]); + } + // No normalization before is_zero() here, because + // we maintain the matrix normalized throughout the + // algorithm. + ex a = m[k*col + k]; + GINAC_ASSERT(!a.is_zero()); + // Subtract the pivot row KJI-style (so: loop by pivot, then + // column, then row) to maximally exploit pivot row zeros (at + // the expense of the pivot column zeros). The speedup compared + // to the usual KIJ order is not really significant though... + for (unsigned r = k + 1; r < row; r++) { + const ex &b = m[r*col + k]; + if (!b.is_zero()) { + ab[r] = b/a; + rowcnt[r]--; + } + } + colcnt[k] = rowcnt[k] = 0; + for (unsigned c = k + 1; c < col; c++) { + const ex &mr0c = m[k*col + c]; + if (mr0c.is_zero()) + continue; + colcnt[c]--; + for (unsigned r = k + 1; r < row; r++) { + if (ab[r].is_zero()) + continue; + bool waszero = m[r*col + c].is_zero(); + m[r*col + c] = (m[r*col + c] - ab[r]*mr0c).normal(); + bool iszero = m[r*col + c].is_zero(); + if (waszero && !iszero) { + rowcnt[r]++; + colcnt[c]++; + } + if (!waszero && iszero) { + rowcnt[r]--; + colcnt[c]--; + } + } + } + for (unsigned r = k + 1; r < row; r++) { + ab[r] = m[r*col + k] = _ex0; + } + } + return colid; +} /** Perform the steps of division free elimination to bring the m x n matrix * into an upper echelon form. diff --git a/ginac/matrix.h b/ginac/matrix.h index 1a95d0aa..cf7dd039 100644 --- a/ginac/matrix.h +++ b/ginac/matrix.h @@ -159,6 +159,7 @@ protected: int gauss_elimination(const bool det = false); int division_free_elimination(const bool det = false); int fraction_free_elimination(const bool det = false); + std::vector markowitz_elimination(unsigned n); int pivot(unsigned ro, unsigned co, bool symbolic = true); void print_elements(const print_context & c, const char *row_start, const char *row_end, const char *row_sep, const char *col_sep) const; -- 2.45.0