isoron/onerow
Archived
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Create repository for single-row wedge cuts

Browse Source
This commit is contained in:
Alinson S. Xavier
2017-11-18 17:01:42 -06:00
commit 1494e2f759
771 changed files with 741801 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
qxx/CMakeLists.txt Normal file
View File

@@ -0,0 +1,21 @@
set(COMMON_SOURCES
src/dlu.cpp
src/dmat.cpp
src/dvec.cpp
src/rational.cpp
src/slu.cpp
src/smat.cpp
src/strprintf.cpp
src/svec.cpp
include/qxx/dlu.hpp
include/qxx/dmat.hpp
include/qxx/dvec.hpp
include/qxx/rational.hpp
include/qxx/slu.hpp
include/qxx/smat.hpp
include/qxx/strprintf.hpp
include/qxx/svec.hpp)
add_library(qxx_static ${COMMON_SOURCES})
set_target_properties(qxx_static PROPERTIES OUTPUT_NAME qxx)
target_include_directories (qxx_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})

54
qxx/include/qxx/debug.hpp Normal file
View File

@@ -0,0 +1,54 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef QXX_DEBUG_H
#define QXX_DEBUG_H
#include <cassert>
#include <cstdio>
// **************************************************************************
//
// **************************************************************************
#define Q_RANGE_CHECK(i, lb, ub) \
do { \
if (((i) < (lb)) || ((i) > (ub))) { \
fprintf(stderr, "%s:%d: in %s(): range check failed " \
"(%s not in %s..%s | %d not in %d..%d)\n", \
__FILE__, __LINE__, __func__, \
#i, #lb, #ub, (i), (lb), (ub)); \
assert(!"range check"); \
} \
} while(0)
#define Q_MATCH(a, b) \
do { \
if ((a) != (b)) { \
fprintf(stderr, "%s:%d: in %s(): match failed " \
"(%s != %s | %d != %d)\n", \
__FILE__, __LINE__, __func__, \
#a, #b, (a), (b)); \
assert(!"match"); \
} \
} while(0)
// **************************************************************************
//
// **************************************************************************
#endif

60
qxx/include/qxx/dlu.hpp Normal file
View File

@@ -0,0 +1,60 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef QXX_DENSE_LU_HPP
#define QXX_DENSE_LU_HPP
#include "rational.hpp"
#include "dmat.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
class dlu {
public:
dlu(const dmat &a);
~dlu();
dmat L() const;
dmat U() const;
dvec solve_Ax(const dvec &b) const;
dvec solve_xA(const dvec &b) const;
mpq det() const;
private:
void pivot(int k);
void factorize(const dmat &a);
int n;
std::vector<int> row_bwd;
dmat g;
};
// **************************************************************************
//
// **************************************************************************
};
#endif

104
qxx/include/qxx/dmat.hpp Normal file
View File

@@ -0,0 +1,104 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef QXX_DMAT_HPP
#define QXX_DMAT_HPP
#include <cstdio>
#include "rational.hpp"
#include "dvec.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
class dmat {
public:
dmat();
dmat(const dmat &b);
dmat(int m, int n);
int rows() const;
int cols() const;
void resize(int m, int n);
void clear();
const mpq &get(int i, int j) const;
void set(int i, int j, const mpq &v);
void set(const mpq &v);
void set_identity();
dvec get_col(int j) const;
void set_col(int j, const dvec &v);
const dvec &get_row(int i) const;
void set_row(int i, const dvec &v);
dvec &operator[] (int i);
const dvec &operator[] (int i) const;
mpq &operator() (int i, int j);
const mpq &operator() (int i, int j) const;
const dmat &operator+() const;
dmat operator-() const;
dmat operator+(const dmat &b) const;
dmat operator-(const dmat &b) const;
dmat operator*(const dmat &b) const;
dvec operator*(const dvec &v) const;
dmat operator*(const mpq &v) const;
dmat operator/(const mpq &v) const;
dmat inv() const;
dmat t() const;
mpq det() const;
bool operator==(const dmat &b) const;
bool operator!=(const dmat &b) const;
dmat &operator=(const dmat &b);
dmat &operator=(const dvec &v);
dmat &operator+=(const dmat &b);
dmat &operator-=(const dmat &b);
dmat &operator*=(const dmat &b);
dmat &operator*=(const mpq &v);
dmat &operator/=(const mpq &v);
void fdump(FILE *f, const std::string &name = "") const;
void dump(const std::string &name = "") const;
private:
int m, n;
std::vector<dvec> d;
};
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const dmat &a);
// **************************************************************************
//
// **************************************************************************
}
#endif

97
qxx/include/qxx/dvec.hpp Normal file
View File

@@ -0,0 +1,97 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef QXX_DENSE_VECTOR_HPP
#define QXX_DENSE_VECTOR_HPP
#include <cstdio>
#include <vector>
#include "rational.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
class svec;
// **************************************************************************
//
// **************************************************************************
class dvec {
public:
dvec();
dvec(const dvec &b);
dvec(int n);
dvec(const svec &b);
int size() const;
void resize(int n);
void clear();
mpq &operator[] (int i);
const mpq &operator[] (int i) const;
void set(const mpq &v);
void set(int i0, int i1, const mpq &v);
void assign(int n, const mpq &v);
dvec operator+(const dvec &b) const;
dvec operator-(const dvec &b) const;
mpq operator*(const dvec &b) const;
mpq operator*(const svec &b) const;
const dvec &operator+() const;
dvec operator-() const;
dvec operator+(const mpq &v) const;
dvec operator-(const mpq &v) const;
dvec operator*(const mpq &v) const;
dvec operator/(const mpq &v) const;
bool operator==(const dvec &b) const;
bool operator!=(const dvec &b) const;
dvec &operator=(const dvec &b);
dvec &operator=(const svec &b);
dvec &operator+=(const dvec &b);
dvec &operator-=(const dvec &b);
dvec &operator+=(const mpq &v);
dvec &operator-=(const mpq &v);
dvec &operator*=(const mpq &v);
dvec &operator/=(const mpq &v);
void fdump(FILE *f, const std::string &name = "") const;
void dump(const std::string &name = "") const;
private:
std::vector<mpq> d;
};
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const dvec &v);
// **************************************************************************
//
// **************************************************************************
}
#endif

114
qxx/include/qxx/rational.hpp Normal file
View File

@@ -0,0 +1,114 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef QXX_TYPES_H
#define QXX_TYPES_H
#include <cstdio>
#include <iostream>
#include <gmp.h>
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
class mpq {
public:
mpq();
mpq(const mpq &b);
mpq(long num, long den);
mpq(int num, int den);
mpq(long i);
mpq(int i);
mpq(double d);
mpq(const char *s);
~mpq();
const mpq &operator+() const;
mpq operator-() const;
mpq operator+(const mpq &b) const;
mpq operator-(const mpq &b) const;
mpq operator*(const mpq &b) const;
mpq operator/(const mpq &b) const;
mpq &operator=(const mpq &b);
mpq &operator=(int i);
mpq &operator+=(const mpq &b);
mpq &operator-=(const mpq &b);
mpq &operator*=(const mpq &b);
mpq &operator/=(const mpq &b);
bool operator<(const mpq &b) const;
bool operator>(const mpq &b) const;
bool operator<=(const mpq &b) const;
bool operator>=(const mpq &b) const;
bool operator==(const mpq &b) const;
bool operator!=(const mpq &b) const;
void set_neg();
void set_neg(const mpq &b);
void set_inv();
void set_inv(const mpq &b);
void set_abs();
void set_abs(const mpq &b);
mpq neg() const;
mpq inv() const;
mpq abs() const;
int sign() const;
mpq num() const;
mpq den() const;
int bits() const;
mpq frac() const;
mpq floor() const;
mpq trunc() const;
mpq ceil() const;
long get_long_num() const;
long get_long_den() const;
double get_double() const;
mpq reduce(mpq max_den) const;
void fdump(FILE *f, const std::string &name = "") const;
void dump(const std::string &name = "") const;
public:
mpq_t v;
};
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const mpq &v);
// **************************************************************************
//
// **************************************************************************
}
#endif

86
qxx/include/qxx/slu.hpp Normal file
View File

@@ -0,0 +1,86 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef QXX_SLU_H
#define QXX_SLU_H
#include <vector>
#include "qxx/smat.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
class perm {
public:
perm();
~perm();
int size() const;
void resize(int n);
void clear();
void id();
void id(int n);
void pivot(int k, int i);
public:
std::vector<int> fwd;
std::vector<int> bwd;
};
// **************************************************************************
//
// **************************************************************************
class slu {
public:
slu(const smat &a);
~slu();
dvec solve_Ax(const dvec &b) const;
dvec solve_xA(const dvec &b) const;
mpq det() const;
public:
void pivot(int k);
void factorize(const smat &a);
dvec solve_gen(const smat &l, const smat &u,
const perm &row, const perm &col, const dvec &b) const;
int n;
std::vector<int> cnz;
perm row, col;
smat w, lt, u;
smat l, ut;
};
// **************************************************************************
//
// **************************************************************************
};
#endif

97
qxx/include/qxx/smat.hpp Normal file
View File

@@ -0,0 +1,97 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef QXX_SMAT_HPP
#define QXX_SMAT_HPP
#include <cstdio>
#include "rational.hpp"
#include "dvec.hpp"
#include "dmat.hpp"
#include "svec.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
class smat {
public:
smat();
smat(const smat &b);
smat(int m, int n);
int rows() const;
int cols() const;
void resize(int m, int n);
void clear();
const mpq &get(int i, int j) const;
void set(int i, int j, const mpq &v);
void set_zero();
void set_identity();
void gather(const dmat &src);
void spread(dmat &r) const;
dmat dense() const;
svec get_col(int j) const;
void set_col(int j, const svec &v);
const svec &get_row(int i) const;
void set_row(int i, const svec &v);
svec &operator[] (int i);
const svec &operator[] (int i) const;
svec_ref operator() (int i, int j);
const mpq &operator() (int i, int j) const;
smat &operator=(const smat &b);
smat operator*(const smat &b) const;
dvec operator*(const dvec &v) const;
svec operator*(const svec &v) const;
smat inv() const;
void set_transpose(const smat &a);
smat t() const;
mpq det() const;
void fdump(FILE *f, const std::string &name = "") const;
void dump(const std::string &name = "") const;
private:
int m, n;
std::vector<svec> d;
};
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const smat &a);
// **************************************************************************
//
// **************************************************************************
}
#endif

31
qxx/include/qxx/strprintf.hpp Normal file
View File

@@ -0,0 +1,31 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef QXX_STRPRINTF_HPP
#define QXX_STRPRINTF_HPP
#include <string>
namespace q {
void format(std::string &s, const char *fmt, ...);
std::string strprintf(const char *fmt, ...);
} // namespace
#endif

178
qxx/include/qxx/svec.hpp Normal file
View File

@@ -0,0 +1,178 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef QXX_SPARSE_VECTOR_HPP
#define QXX_SPARSE_VECTOR_HPP
#include <cstdio>
#include "rational.hpp"
#include "dvec.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
class sparse_el;
class svec_ref;
class svec;
// **************************************************************************
//
// **************************************************************************
class svec_ref {
public:
svec_ref(svec &vec, int i);
svec_ref(const svec_ref &src);
~svec_ref();
const mpq &get() const;
operator mpq() const;
const mpq &operator=(const svec_ref &src) const;
const mpq &operator=(const mpq &b) const;
const mpq &operator+=(const mpq &b) const;
const mpq &operator-=(const mpq &b) const;
const mpq &operator*=(const mpq &b) const;
const mpq &operator/=(const mpq &b) const;
const mpq &operator+() const;
mpq operator-() const;
mpq operator+(const mpq &b) const;
mpq operator-(const mpq &b) const;
mpq operator*(const mpq &b) const;
mpq operator/(const mpq &b) const;
bool operator<(const mpq &b) const;
bool operator>(const mpq &b) const;
bool operator<=(const mpq &b) const;
bool operator>=(const mpq &b) const;
bool operator==(const mpq &b) const;
bool operator!=(const mpq &b) const;
public:
svec &vector;
int index;
};
// **************************************************************************
//
// **************************************************************************
class svec_el {
public:
svec_el();
svec_el(int idx, const mpq &v);
svec_el(const svec_el &e);
~svec_el();
public:
int index;
mpq value;
};
// **************************************************************************
//
// **************************************************************************
class svec {
public:
svec();
svec(const svec &src);
svec(int n);
svec(const dvec &src);
~svec();
void clear();
int size() const;
void resize(int n);
int nz() const;
bool offs_active() const;
void offs_build();
void offs_clear();
void set_zero();
void gather(const dvec &src);
void spread(dvec &r) const;
dvec dense() const;
int locate(int idx) const;
void remove(int offset);
int &index(int offset);
mpq &value(int offset);
int index(int offset) const;
const mpq &value(int offset) const;
void push_nz(int idx, const mpq &v);
void push(int idx, const mpq &v);
const mpq &get(int idx) const;
const mpq &set(int idx, const mpq &v);
const mpq &operator[](int idx) const;
svec_ref operator[](int idx);
const svec &operator+() const;
svec operator-() const;
svec operator+(const svec &b) const;
svec operator-(const svec &b) const;
mpq operator*(const svec &b) const;
mpq operator*(const dvec &b) const;
svec operator*(const mpq &v);
svec operator/(const mpq &v);
svec &operator=(const svec &b);
svec &operator+=(const svec &b);
svec &operator-=(const svec &b);
svec &operator*=(const mpq &v);
svec &operator/=(const mpq &v);
void sort();
void fdump(FILE *f, const std::string &name = "") const;
void dump(const std::string &name = "") const;
private:
int n;
std::vector<svec_el> d;
std::vector<int> offs;
static const mpq zero;
};
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const svec &v);
// **************************************************************************
//
// **************************************************************************
}
#endif

240
qxx/src/dlu.cpp Normal file
View File

@@ -0,0 +1,240 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cassert>
#include "qxx/dlu.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
dlu::dlu(const dmat &a)
: n(0)
{
factorize(a);
}
dlu::~dlu()
{
}
// **************************************************************************
//
// **************************************************************************
void dlu::pivot(int p)
{
int best_v = -1;
int best_k = -1;
//printf("pivot(%d):\n", p);
for (int k = p; k < n; k++) {
int i = row_bwd[k];
//gmp_printf("\tchoice [%d] row %d: %Qd\n",
// k, i, g(i, p).v);
if (g(i, p).sign()) {
int v = g(i, p).bits();
if ((best_v < 0) || (v < best_v)) {
best_k = k;
best_v = v;
}
}
}
assert((best_k >= 0) && "Singular matrix");
int i = row_bwd[best_k];
row_bwd[best_k] = row_bwd[p];
row_bwd[p] = i;
}
// **************************************************************************
//
// **************************************************************************
void dlu::factorize(const dmat &a)
{
// check
assert(a.rows() == a.cols());
// init
n = a.rows();
g = a;
row_bwd.resize(n);
for (int p = 0; p < n; p++)
row_bwd[p] = p;
// main loop
mpq pv, f;
for (int p = 0; p < n; p++) {
pivot(p);
int pi = row_bwd[p];
pv = g(pi, p);
for (int k = p + 1; k < n; k++) {
int i = row_bwd[k];
f = g(i, p) / pv;
for (int j = p + 1; j < n; j++)
g(i, j) -= f * g(pi, j);
}
//printf("g%d = ", p);
//g.dump(stdout);
}
//printf("row_bwd = [");
//for (int p = 0; p < n; p++)
// printf(" %d", row_bwd[p]);
//printf(" ];\n");
}
// **************************************************************************
//
// **************************************************************************
dmat dlu::L() const
{
dmat l(n, n);
for (int k = 0; k < n; k++) {
for (int j = 0; j < k; j++)
l(k, j) = g(row_bwd[k], j) / g(row_bwd[j], j);
l(k, k) = 1;
}
return(l);
}
dmat dlu::U() const
{
dmat u(n, n);
for (int k = 0; k < n; k++) {
for (int j = k; j < n; j++)
u(k, j) = g(row_bwd[k], j);
}
return(u);
}
// **************************************************************************
//
// **************************************************************************
dvec dlu::solve_Ax(const dvec &b) const
{
assert(b.size() == n);
dvec x(n);
mpq v;
int i0 = row_bwd[0];
x[0] = b[i0] / g(i0, 0);
for (int k = 1; k < n; k++) {
int i = row_bwd[k];
v = 0;
for (int j = 0; j < k; j++)
v += g(i, j) * x[j];
x[k] = (b[i] - v) / g(i, k);
}
for (int k = n - 2; k != -1; k--) {
int i = row_bwd[k];
v = 0;
for (int j = k + 1; j < n; j++)
v += g(i, j) * x[j];
x[k] -= v / g(i, k);
}
return(x);
}
// **************************************************************************
//
// **************************************************************************
dvec dlu::solve_xA(const dvec &b) const
{
assert(b.size() == n);
dvec x(n);
mpq v;
int i0 = row_bwd[0];
x[i0] = b[0] / g(i0, 0);
for (int k = 1; k < n; k++) {
int i = row_bwd[k];
v = 0;
for (int l = 0; l < k; l++) {
int il = row_bwd[l];
v += g(il, k) * x[il];
}
x[i] = (b[k] - v) / g(i, k);
}
for (int k = n - 2; k != -1; k--) {
int i = row_bwd[k];
v = 0;
for (int l = k + 1; l < n; l++) {
int il = row_bwd[l];
v += g(il, k) * x[il];
}
x[i] -= v / g(i, k);
}
return(x);
}
// **************************************************************************
//
// **************************************************************************
mpq dlu::det() const
{
mpq v;
v = 1;
for (int k = 0; k < n; k++) {
int i = row_bwd[k];
v *= g(i, k);
}
return(v);
}
// **************************************************************************
//
// **************************************************************************
};

486
qxx/src/dmat.cpp Normal file
View File

@@ -0,0 +1,486 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cstdio>
#include "qxx/debug.hpp"
#include "qxx/rational.hpp"
#include "qxx/dvec.hpp"
#include "qxx/dlu.hpp"
#include "qxx/dmat.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
dmat::dmat()
: m(0), n(0)
{
}
dmat::dmat(const dmat &b)
: m(b.m), n(b.n), d(b.d)
{
}
dmat::dmat(int m0, int n0)
: m(m0), n(n0), d(m0)
{
for (int i = 0; i < m; i++)
d[i].resize(n);
}
int dmat::rows() const
{
return(m);
}
int dmat::cols() const
{
return(n);
}
void dmat::resize(int m0, int n0)
{
m = m0;
n = n0;
d.resize(m);
for (int i = 0; i < m; i++)
d[i].resize(n);
}
void dmat::clear()
{
m = n = 0;
d.clear();
}
// **************************************************************************
//
// **************************************************************************
const mpq &dmat::get(int i, int j) const
{
Q_RANGE_CHECK(i, 0, m - 1);
Q_RANGE_CHECK(j, 0, n - 1);
return(d[i][j]);
}
void dmat::set(int i, int j, const mpq &v)
{
Q_RANGE_CHECK(i, 0, m - 1);
Q_RANGE_CHECK(j, 0, n - 1);
d[i][j] = v;
}
void dmat::set(const mpq &v)
{
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
d[i][j] = v;
}
void dmat::set_identity()
{
mpq zero(0), one(1);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
d[i][j] = (i == j) ? one : zero;
}
// **************************************************************************
//
// **************************************************************************
dvec dmat::get_col(int j) const
{
Q_RANGE_CHECK(j, 0, n - 1);
dvec r(m);
for (int i = 0; i < m; i++)
r[i] = d[i][j];
return(r);
}
void dmat::set_col(int j, const dvec &v)
{
Q_RANGE_CHECK(j, 0, n - 1);
Q_MATCH(v.size(), m);
for (int i = 0; i < m; i++)
d[i][j] = v[i];
}
const dvec &dmat::get_row(int i) const
{
Q_RANGE_CHECK(i, 0, m - 1);
return(d[i]);
}
void dmat::set_row(int i, const dvec &v)
{
Q_RANGE_CHECK(i, 0, m - 1);
Q_MATCH(v.size(), n);
d[i] = v;
}
// **************************************************************************
//
// **************************************************************************
dvec &dmat::operator[] (int i)
{
Q_RANGE_CHECK(i, 0, m - 1);
return(d[i]);
}
const dvec &dmat::operator[] (int i) const
{
Q_RANGE_CHECK(i, 0, m - 1);
return(d[i]);
}
mpq &dmat::operator() (int i, int j)
{
Q_RANGE_CHECK(i, 0, m - 1);
Q_RANGE_CHECK(j, 0, n - 1);
return(d[i][j]);
}
const mpq &dmat::operator() (int i, int j) const
{
Q_RANGE_CHECK(i, 0, m - 1);
Q_RANGE_CHECK(j, 0, n - 1);
return(d[i][j]);
}
// **************************************************************************
//
// **************************************************************************
const dmat &dmat::operator+() const
{
return(*this);
}
dmat dmat::operator-() const
{
dmat r(m, n);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
r.d[i][j] = -d[i][j];
return(r);
}
dmat dmat::operator+(const dmat &b) const
{
Q_MATCH(m, b.m);
Q_MATCH(n, b.n);
dmat r(m, n);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
r.d[i][j] = d[i][j] + b.d[i][j];
return(r);
}
dmat dmat::operator-(const dmat &b) const
{
Q_MATCH(m, b.m);
Q_MATCH(n, b.n);
dmat r(m, n);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
r.d[i][j] = d[i][j] - b.d[i][j];
return(r);
}
dmat dmat::operator*(const dmat &b) const
{
Q_MATCH(n, b.m);
dmat r(m, b.n);
mpq v;
for (int i = 0; i < m; i++) {
for (int j = 0; j < b.n; j++) {
v = 0;
for (int k = 0; k < n; k++)
v += d[i][k] * b.d[k][j];
r.d[i][j] = v;
}
}
return(r);
}
dvec dmat::operator*(const dvec &b) const
{
Q_MATCH(n, b.size());
dvec r(m);
mpq v;
for (int i = 0; i < m; i++) {
v = 0;
for (int j = 0; j < n; j++)
v += d[i][j] * b[j];
r[i] = v;
}
return(r);
}
dmat dmat::operator*(const mpq &v) const
{
dmat r(m, n);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
r.d[i][j] = d[i][j] * v;
return(r);
}
dmat dmat::operator/(const mpq &v) const
{
dmat r(m, n);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
r.d[i][j] = d[i][j] / v;
return(r);
}
// **************************************************************************
//
// **************************************************************************
dmat dmat::inv() const
{
Q_MATCH(m, n);
dvec e(n);
dmat r(n, n);
dlu lu(*this);
for (int i = 0; i < n; i++) {
e[i] = 1;
r.set_col(i, lu.solve_Ax(e));
//e.dump(stdout, "e");
//r.get_col(i).dump(stdout, "x");
e[i] = 0;
}
return(r);
}
dmat dmat::t() const
{
dmat r(m, n);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
r.d[i][j] = d[j][i];
return(r);
}
mpq dmat::det() const
{
Q_MATCH(m, n);
if (n == 1)
return(d[0][0]);
if (n == 2)
return(d[0][0] * d[1][1] - d[1][0] * d[0][1]);
if (n == 3)
return( d[0][0] * d[1][1] * d[2][2]
+ d[0][1] * d[1][2] * d[2][0]
+ d[0][2] * d[1][0] * d[2][1]
- d[2][0] * d[1][1] * d[0][2]
- d[2][1] * d[1][2] * d[0][0]
- d[2][2] * d[1][0] * d[0][1] );
dlu lu(*this);
return(lu.det());
}
// **************************************************************************
//
// **************************************************************************
bool dmat::operator==(const dmat &b) const
{
Q_MATCH(m, b.m);
Q_MATCH(n, b.n);
for (int i = 0; i < m; i++) {
if (d[i] != b.d[i])
return(false);
}
return(true);
}
bool dmat::operator!=(const dmat &b) const
{
return(!((*this) == b));
}
// **************************************************************************
//
// **************************************************************************
dmat &dmat::operator=(const dmat &b)
{
m = b.m;
n = b.n;
d = b.d;
return(*this);
}
dmat &dmat::operator=(const dvec &v)
{
resize(v.size(), 1);
for (int i = 0; i < m; i++)
d[i][0] = v[i];
return(*this);
}
dmat &dmat::operator+=(const dmat &b)
{
Q_MATCH(m, b.m);
Q_MATCH(n, b.n);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
d[i][j] += b.d[i][j];
return(*this);
}
dmat &dmat::operator-=(const dmat &b)
{
Q_MATCH(m, b.m);
Q_MATCH(n, b.n);
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
d[i][j] -= b.d[i][j];
return(*this);
}
dmat &dmat::operator*=(const dmat &b)
{
dmat r;
r = (*this) * b;
(*this) = r;
return(*this);
}
dmat &dmat::operator*=(const mpq &v)
{
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
d[i][j] *= v;
return(*this);
}
dmat &dmat::operator/=(const mpq &v)
{
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
d[i][j] /= v;
return(*this);
}
// **************************************************************************
//
// **************************************************************************
void dmat::fdump(FILE *f, const std::string &name) const
{
if (name.length())
gmp_fprintf(f, "%s = [\n", name.c_str());
else
gmp_fprintf(f, "[\n");
for (int i = 0; i < m; i++)
d[i].fdump(f);
fprintf(f, "];\n");
}
void dmat::dump(const std::string &name) const
{
fdump(stdout, name);
}
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const dmat &a)
{
for (int i = 0; i < a.rows(); i++)
os << a[i] << std::endl;
return(os);
}
// **************************************************************************
//
// **************************************************************************
}

366
qxx/src/dvec.cpp Normal file
View File

@@ -0,0 +1,366 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#include "qxx/debug.hpp"
#include "qxx/rational.hpp"
#include "qxx/dvec.hpp"
#include "qxx/svec.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
dvec::dvec()
{
}
dvec::dvec(const dvec &b)
: d(b.d)
{
}
dvec::dvec(int n)
: d(n)
{
}
dvec::dvec(const svec &b)
{
b.spread(*this);
}
// **************************************************************************
//
// **************************************************************************
int dvec::size() const
{
return(d.size());
}
void dvec::resize(int n)
{
d.resize(n);
}
void dvec::clear()
{
d.clear();
}
mpq &dvec::operator[] (int i)
{
Q_RANGE_CHECK(i, 0, size() - 1);
return(d[i]);
}
const mpq &dvec::operator[] (int i) const
{
Q_RANGE_CHECK(i, 0, size() - 1);
return(d[i]);
}
// **************************************************************************
//
// **************************************************************************
void dvec::set(const mpq &v)
{
set(0, size() - 1, v);
}
void dvec::set(int i0, int i1, const mpq &v)
{
for (int i = i0; i <= i1; i++)
d[i] = v;
}
void dvec::assign(int n, const mpq &v)
{
resize(n);
set(v);
}
// **************************************************************************
//
// **************************************************************************
dvec dvec::operator+(const dvec &b) const
{
Q_MATCH(size(), b.size());
int n = size();
dvec r(n);
for (int i = 0; i < n; i++)
r.d[i] = d[i] + b.d[i];
return(r);
}
dvec dvec::operator-(const dvec &b) const
{
Q_MATCH(size(), b.size());
int n = size();
dvec r(n);
for (int i = 0; i < n; i++)
r.d[i] = d[i] - b.d[i];
return(r);
}
mpq dvec::operator*(const dvec &b) const
{
Q_MATCH(size(), b.size());
int n = size();
mpq r;
for (int i = 0; i < n; i++)
r += d[i] * b.d[i];
return(r);
}
mpq dvec::operator*(const svec &b) const
{
return(b.operator*(*this));
}
// **************************************************************************
//
// **************************************************************************
const dvec &dvec::operator+() const
{
return(*this);
}
dvec dvec::operator-() const
{
int n = size();
dvec r(n);
for (int i = 0; i < n; i++)
r.d[i] = -d[i];
return(r);
}
dvec dvec::operator+(const mpq &v) const
{
int n = size();
dvec r(n);
for (int i = 0; i < n; i++)
r.d[i] = d[i] + v;
return(r);
}
dvec dvec::operator-(const mpq &v) const
{
int n = size();
dvec r(n);
for (int i = 0; i < n; i++)
r.d[i] = d[i] - v;
return(r);
}
dvec dvec::operator*(const mpq &v) const
{
int n = size();
dvec r(n);
for (int i = 0; i < n; i++)
r.d[i] = d[i] * v;
return(r);
}
dvec dvec::operator/(const mpq &v) const
{
int n = size();
dvec r(n);
for (int i = 0; i < n; i++)
r.d[i] = d[i] / v;
return(r);
}
// **************************************************************************
//
// **************************************************************************
bool dvec::operator==(const dvec &b) const
{
Q_MATCH(size(), b.size());
int n = size();
for (int i = 0; i < n; i++) {
if (d[i] != b.d[i])
return(false);
}
return(true);
}
bool dvec::operator!=(const dvec &b) const
{
return(!((*this) == b));
}
// **************************************************************************
//
// **************************************************************************
dvec &dvec::operator=(const dvec &b)
{
d = b.d;
return(*this);
}
dvec &dvec::operator=(const svec &b)
{
b.spread(*this);
return(*this);
}
// **************************************************************************
//
// **************************************************************************
dvec &dvec::operator+=(const dvec &b)
{
Q_MATCH(size(), b.size());
int n = size();
for (int i = 0; i < n; i++)
d[i] += b.d[i];
return(*this);
}
dvec &dvec::operator-=(const dvec &b)
{
Q_MATCH(size(), b.size());
int n = size();
for (int i = 0; i < n; i++)
d[i] -= b.d[i];
return(*this);
}
dvec &dvec::operator+=(const mpq &v)
{
int n = size();
for (int i = 0; i < n; i++)
d[i] += v;
return(*this);
}
dvec &dvec::operator-=(const mpq &v)
{
int n = size();
for (int i = 0; i < n; i++)
d[i] -= v;
return(*this);
}
dvec &dvec::operator*=(const mpq &v)
{
int n = size();
for (int i = 0; i < n; i++)
d[i] *= v;
return(*this);
}
dvec &dvec::operator/=(const mpq &v)
{
int n = size();
for (int i = 0; i < n; i++)
d[i] /= v;
return(*this);
}
// **************************************************************************
//
// **************************************************************************
void dvec::fdump(FILE *f, const std::string &name) const
{
int n = size();
if (name.length())
gmp_fprintf(f, "%s = [", name.c_str());
for (int i = 0; i < n; i++) {
d[i].fdump(f);
}
if (name.length())
gmp_fprintf(f, " ]';\n");
else
gmp_fprintf(f, ";\n");
}
void dvec::dump(const std::string &name) const
{
fdump(stdout, name);
}
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const dvec &v)
{
int n = v.size();
for (int i = 0; i < n; i++)
os << " " << v[i];
return(os);
}
// **************************************************************************
//
// **************************************************************************
} // namespace

437
qxx/src/rational.cpp Normal file
View File

@@ -0,0 +1,437 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cassert>
#include "qxx/rational.hpp"
#include "qxx/strprintf.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
mpq::mpq()
{
mpq_init(v);
}
mpq::mpq(const mpq &b)
{
mpq_init(v);
mpq_set(v, b.v);
}
mpq::mpq(long num, long den)
{
mpq_init(v);
mpq_set_si(v, num, den);
mpq_canonicalize(v);
}
mpq::mpq(int num, int den)
{
mpq_init(v);
mpq_set_si(v, num, den);
mpq_canonicalize(v);
}
mpq::mpq(long i)
{
mpq_init(v);
mpq_set_si(v, i, 1);
}
mpq::mpq(int i)
{
mpq_init(v);
mpq_set_si(v, i, 1);
}
mpq::mpq(double d)
{
mpq_init(v);
mpq_set_d(v, d);
}
mpq::mpq(const char *s)
{
mpq_init(v);
mpq_set_str(v, s, 0);
mpq_canonicalize(v);
}
mpq::~mpq()
{
mpq_clear(v);
}
// **************************************************************************
//
// **************************************************************************
const mpq &mpq::operator+() const
{
return(*this);
}
mpq mpq::operator-() const
{
mpq r;
mpq_neg(r.v, v);
return(r);
}
mpq mpq::operator+(const mpq &b) const
{
mpq r;
mpq_add(r.v, v, b.v);
return(r);
}
mpq mpq::operator-(const mpq &b) const
{
mpq r;
mpq_sub(r.v, v, b.v);
return(r);
}
mpq mpq::operator*(const mpq &b) const
{
mpq r;
mpq_mul(r.v, v, b.v);
return(r);
}
mpq mpq::operator/(const mpq &b) const
{
mpq r;
mpq_div(r.v, v, b.v);
return(r);
}
mpq &mpq::operator=(const mpq &b)
{
mpq_set(v, b.v);
return(*this);
}
mpq &mpq::operator=(int i)
{
mpq_set_si(v, i, 1);
return(*this);
}
mpq &mpq::operator+=(const mpq &b)
{
mpq_add(v, v, b.v);
return(*this);
}
mpq &mpq::operator-=(const mpq &b)
{
mpq_sub(v, v, b.v);
return(*this);
}
mpq &mpq::operator*=(const mpq &b)
{
mpq_mul(v, v, b.v);
return(*this);
}
mpq &mpq::operator/=(const mpq &b)
{
mpq_div(v, v, b.v);
return(*this);
}
bool mpq::operator<(const mpq &b) const
{
return(mpq_cmp(v, b.v) < 0);
}
bool mpq::operator>(const mpq &b) const
{
return(mpq_cmp(v, b.v) > 0);
}
bool mpq::operator<=(const mpq &b) const
{
return(mpq_cmp(v, b.v) <= 0);
}
bool mpq::operator>=(const mpq &b) const
{
return(mpq_cmp(v, b.v) >= 0);
}
bool mpq::operator==(const mpq &b) const
{
return(mpq_equal(v, b.v));
}
bool mpq::operator!=(const mpq &b) const
{
return(!mpq_equal(v, b.v));
}
// **************************************************************************
//
// **************************************************************************
void mpq::set_neg()
{
mpq_neg(v, v);
}
void mpq::set_neg(const mpq &b)
{
mpq_neg(v, b.v);
}
void mpq::set_inv()
{
mpq_inv(v, v);
}
void mpq::set_inv(const mpq &b)
{
mpq_inv(v, b.v);
}
void mpq::set_abs()
{
mpq_abs(v, v);
}
void mpq::set_abs(const mpq &b)
{
mpq_abs(v, b.v);
}
mpq mpq::neg() const
{
mpq r;
mpq_neg(r.v, v);
return(r);
}
mpq mpq::inv() const
{
mpq r;
mpq_inv(r.v, v);
return(r);
}
mpq mpq::abs() const
{
mpq r;
mpq_abs(r.v, v);
return(r);
}
int mpq::sign() const
{
return(mpq_sgn(v));
}
// **************************************************************************
//
// **************************************************************************
mpq mpq::num() const
{
mpq r;
mpz_set(mpq_numref(r.v), mpq_numref(v));
mpz_set_ui(mpq_denref(r.v), 1);
return(r);
}
mpq mpq::den() const
{
mpq r;
mpz_set(mpq_numref(r.v), mpq_denref(v));
mpz_set_ui(mpq_denref(r.v), 1);
return(r);
}
int mpq::bits() const
{
return(
mpz_sizeinbase(mpq_numref(v), 2)
+ mpz_sizeinbase(mpq_denref(v), 2)
);
}
mpq mpq::frac() const
{
mpq r;
mpz_fdiv_r(mpq_numref(r.v), mpq_numref(v), mpq_denref(v));
mpz_set(mpq_denref(r.v), mpq_denref(v));
return(r);
}
mpq mpq::floor() const
{
mpq r;
mpz_fdiv_q(mpq_numref(r.v), mpq_numref(v), mpq_denref(v));
mpz_set_ui(mpq_denref(r.v), 1);
return(r);
}
mpq mpq::trunc() const
{
mpq r;
mpz_tdiv_q(mpq_numref(r.v), mpq_numref(v), mpq_denref(v));
mpz_set_ui(mpq_denref(r.v), 1);
return(r);
}
mpq mpq::ceil() const
{
mpq r;
mpz_cdiv_q(mpq_numref(r.v), mpq_numref(v), mpq_denref(v));
mpz_set_ui(mpq_denref(r.v), 1);
return(r);
}
// **************************************************************************
//
// **************************************************************************
long mpq::get_long_num() const
{
assert(mpz_fits_slong_p(mpq_numref(v)));
return(mpz_get_si(mpq_numref(v)));
}
long mpq::get_long_den() const
{
assert(mpz_fits_slong_p(mpq_denref(v)));
return(mpz_get_si(mpq_denref(v)));
}
double mpq::get_double() const
{
return(mpq_get_d(v));
}
// **************************************************************************
//
// **************************************************************************
mpq mpq::reduce(mpq max_coeff) const
{
mpq p[3], q[3];
mpq f, ff, r;
p[0] = 0;
q[0] = 1;
p[1] = 1;
q[1] = 0;
f = *this;
while (1) {
ff = f.floor();
r = f - ff;
p[2] = ff * p[1] + p[0];
q[2] = ff * q[1] + q[0];
if ((p[2].abs() > max_coeff) || (q[2] > max_coeff))
break;
if (r.sign() == 0) {
p[1] = p[2];
q[1] = q[2];
break;
}
f = r.inv();
p[0] = p[1];
q[0] = q[1];
p[1] = p[2];
q[1] = q[2];
}
if (q[1].sign() == 0)
q[1] = 1;
return(p[1] / q[1]);
}
// **************************************************************************
//
// **************************************************************************
void mpq::fdump(FILE *f, const std::string &name) const
{
mpq_t t;
mpq_init(t);
mpq_set(t, v);
mpq_canonicalize(t);
if (mpq_cmp(t, v) != 0) {
gmp_fprintf(f, "\nBUG: not canon.: %Qd\n",
mpq_numref(v),
mpq_denref(v));
}
mpq_clear(t);
if (name.length())
gmp_fprintf(f, "%s = %Qd;\n", name.c_str(), v);
else
gmp_fprintf(f, " %Qd", v);
}
void mpq::dump(const std::string &name) const
{
fdump(stdout, name);
}
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const mpq &v)
{
os << strprintf("%Qd", v.v);
return(os);
}
// **************************************************************************
//
// **************************************************************************
}

307
qxx/src/slu.cpp Normal file
View File

@@ -0,0 +1,307 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cassert>
#include "qxx/strprintf.hpp"
#include "qxx/rational.hpp"
#include "qxx/slu.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
perm::perm()
{
}
perm::~perm()
{
}
int perm::size() const
{
return(fwd.size());
}
void perm::resize(int n)
{
fwd.resize(n);
bwd.resize(n);
}
void perm::clear()
{
fwd.clear();
bwd.clear();
}
void perm::id()
{
int n = fwd.size();
for (int i = 0; i < n; i++)
fwd[i] = i;
for (int k = 0; k < n; k++)
bwd[k] = k;
}
void perm::id(int n)
{
resize(n);
id();
}
void perm::pivot(int k, int i)
{
int t = fwd[i];
int u = bwd[k];
fwd[i] = k;
fwd[u] = t;
bwd[k] = i;
bwd[t] = u;
}
// **************************************************************************
//
// **************************************************************************
slu::slu(const smat &a)
{
factorize(a);
}
slu::~slu()
{
}
// **************************************************************************
//
// **************************************************************************
void slu::pivot(int p)
{
double best_v = n * n;
int best_i = -1;
int best_j = -1;
for (int k = p; k < n; k++) {
int i = row.bwd[k];
int inz = w[i].nz();
for (int l = 0; l < inz; l++) {
int j = w[i].index(l);
int jnz = cnz[j];
int kj = col.fwd[j];
double v = (double)(inz - 1) * (jnz - 1);
if ((kj >= p) && (v < best_v)) {
best_v = v;
best_i = i;
best_j = j;
}
}
}
assert((best_i >= 0) && "Singular matrix");
//best_i = best_j = p;
row.pivot(p, best_i);
col.pivot(p, best_j);
}
// **************************************************************************
//
// **************************************************************************
void slu::factorize(const smat &a)
{
// check
assert(a.rows() == a.cols());
// init
n = a.rows();
w = a;
lt.resize(n, n);
u.resize(n, n);
row.id(n);
col.id(n);
cnz.assign(n, 0);
for (int i = 0; i < n; i++) {
for (int l = 0; l < w[i].nz(); l++)
cnz[w[i].index(l)]++;
}
// main loop
svec pr;
mpq pv, f;
std::vector<bool> hit;
for (int p = 0; p < n; p++) {
//w.dense().dump(strprintf("w%d", p).c_str());
pivot(p);
int pi = row.bwd[p];
int pj = col.bwd[p];
// get pivot row
pr = w[pi];
pr.offs_build();
pv = pr[pj];
for (int pl = 0; pl < pr.nz(); pl++)
cnz[pr.index(pl)]--;
// update u, lt, w
u[p] = w[pi];
lt[p] = w.get_col(pj) / pv;
w[pi].set_zero();
// elimination
for (int k = p + 1; k < n; k++) {
int i = row.bwd[k];
svec &ir = w[i];
f = -ir[pj] / pv;
if (f.sign() == 0)
continue;
hit.assign(pr.nz(), false);
for (int l = 0; l < ir.nz(); l++) {
int j = ir.index(l);
int pl = pr.locate(j);
if (pl == -1)
continue;
hit[pl] = true;
ir.value(l) += f * pr.value(pl);
if (ir.value(l).sign() == 0) {
cnz[j]--;
ir.remove(l);
l--;
}
}
for (int pl = 0; pl < pr.nz(); pl++) {
if (hit[pl])
continue;
int j = pr.index(pl);
ir.push(j, f * pr.value(pl));
cnz[j]++;
}
}
}
// postprocessing
for (int p = 0; p < n; p++) {
svec &pr = u[p];
for (int l = 0; l < pr.nz(); l++)
pr.index(l) = col.fwd[pr.index(l)];
}
for (int p = 0; p < n; p++) {
svec &pc = lt[p];
for (int l = 0; l < pc.nz(); l++)
pc.index(l) = row.fwd[pc.index(l)];
}
l.set_transpose(lt);
ut.set_transpose(u);
}
// **************************************************************************
//
// **************************************************************************
dvec slu::solve_gen(const smat &l, const smat &u,
const perm &row, const perm &col, const dvec &b) const
{
dvec x(n), y(n);
mpq v, diag;
// Ly = b
for (int k = 0; k < n; k++) {
v = b[row.bwd[k]];
const svec &lr = l[k];
for (int l = 0; l < lr.nz(); l++) {
int j = lr.index(l);
if (j == k)
diag = lr.value(l);
else
v -= lr.value(l) * y[j];
}
y[k] = v / diag;
}
// Ux = y
for (int k = n - 1; k != -1; k--) {
v = y[k];
const svec &ur = u[k];
for (int l = 0; l < ur.nz(); l++) {
int j = ur.index(l);
if (j == k)
diag = ur.value(l);
else
v -= ur.value(l) * y[j];
}
y[k] = v / diag;
}
// postprocessing
for (int j = 0; j < n; j++)
x[j] = y[col.fwd[j]];
return(x);
}
dvec slu::solve_Ax(const dvec &b) const
{
return(solve_gen(l, u, row, col, b));
}
dvec slu::solve_xA(const dvec &b) const
{
return(solve_gen(ut, lt, col, row, b));
}
// **************************************************************************
//
// **************************************************************************
};

339
qxx/src/smat.cpp Normal file
View File

@@ -0,0 +1,339 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cstdio>
#include "qxx/debug.hpp"
#include "qxx/rational.hpp"
#include "qxx/dvec.hpp"
#include "qxx/slu.hpp"
#include "qxx/smat.hpp"
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
//
// **************************************************************************
smat::smat()
: m(0), n(0)
{
}
smat::smat(const smat &b)
: m(b.m), n(b.n), d(b.d)
{
}
smat::smat(int m0, int n0)
: m(m0), n(n0), d(m0)
{
for (int i = 0; i < m; i++)
d[i].resize(n);
}
int smat::rows() const
{
return(m);
}
int smat::cols() const
{
return(n);
}
void smat::resize(int m0, int n0)
{
m = m0;
n = n0;
d.resize(m);
for (int i = 0; i < m; i++)
d[i].resize(n);
}
void smat::clear()
{
m = n = 0;
d.clear();
}
// **************************************************************************
//
// **************************************************************************
const mpq &smat::get(int i, int j) const
{
Q_RANGE_CHECK(i, 0, m - 1);
Q_RANGE_CHECK(j, 0, n - 1);
return(d[i][j]);
}
void smat::set(int i, int j, const mpq &v)
{
Q_RANGE_CHECK(i, 0, m - 1);
Q_RANGE_CHECK(j, 0, n - 1);
d[i][j] = v;
}
void smat::set_zero()
{
for (int i = 0; i < m; i++)
d[i].set_zero();
}
void smat::set_identity()
{
mpq one(1);
for (int i = 0; i < m; i++) {
d[i].set_zero();
d[i][i] = one;
}
}
// **************************************************************************
//
// **************************************************************************
void smat::gather(const dmat &src)
{
resize(src.rows(), src.cols());
for (int i = 0; i < m; i++)
d[i].gather(src[i]);
}
void smat::spread(dmat &r) const
{
r.resize(m, n);
for (int i = 0; i < m; i++)
d[i].spread(r[i]);
}
dmat smat::dense() const
{
dmat r;
spread(r);
return(r);
}
// **************************************************************************
//
// **************************************************************************
svec smat::get_col(int j) const
{
svec r(m);
for (int i = 0; i < m; i++)
r.push(i, d[i][j]);
return(r);
}
void smat::set_col(int j, const svec &v)
{
Q_MATCH(v.size(), m);
for (int i = 0; i < m; i++)
d[i][j] = v[i];
}
const svec &smat::get_row(int i) const
{
return(d[i]);
}
void smat::set_row(int i, const svec &v)
{
Q_MATCH(v.size(), n);
d[i] = v;
}
// **************************************************************************
//
// **************************************************************************
svec &smat::operator[] (int i)
{
return(d[i]);
}
const svec &smat::operator[] (int i) const
{
return(d[i]);
}
// **************************************************************************
//
// **************************************************************************
svec_ref smat::operator() (int i, int j)
{
Q_RANGE_CHECK(i, 0, m - 1);
Q_RANGE_CHECK(j, 0, n - 1);
return(d[i][j]);
}
const mpq &smat::operator() (int i, int j) const
{
Q_RANGE_CHECK(i, 0, m - 1);
Q_RANGE_CHECK(j, 0, n - 1);
return(d[i][j]);
}
// **************************************************************************
//
// **************************************************************************
smat &smat::operator=(const smat &a)
{
m = a.m;
n = a.n;
d = a.d;
return(*this);
}
// **************************************************************************
//
// **************************************************************************
dvec smat::operator*(const dvec &b) const
{
Q_MATCH(n, b.size());
dvec x(m);
mpq v;
for (int i = 0; i < m; i++) {
v = 0;
const svec &r = d[i];
for (int l = 0; l < r.nz(); l++)
v += r.value(l) * b[r.index(l)];
x[i] = v;
}
return(x);
}
// **************************************************************************
//
// **************************************************************************
void smat::set_transpose(const smat &a)
{
resize(a.n, a.m);
for (int i = 0; i < a.n; i++) {
const svec &ar = a[i];
for (int l = 0; l < ar.nz(); l++)
d[ar.index(l)].push(i, ar.value(l));
}
}
// **************************************************************************
//
// **************************************************************************
smat smat::t() const
{
smat a;
a.set_transpose(*this);
return(a);
}
// **************************************************************************
//
// **************************************************************************
smat smat::inv() const
{
Q_MATCH(m, n);
slu lu(*this);
smat g;
dvec e;
svec x;
g.resize(n, n);
e.assign(n, 0);
for (int i = 0; i < n; i++) {
e[i] = 1;
x = lu.solve_xA(e);
g.set_row(i, x);
e[i] = 0;
}
return(g);
}
// **************************************************************************
//
// **************************************************************************
void smat::fdump(FILE *f, const std::string &name) const
{
if (name.length())
gmp_fprintf(f, "%s = [\n", name.c_str());
else
gmp_fprintf(f, "[\n");
for (int i = 0; i < m; i++)
d[i].fdump(f);
gmp_fprintf(f, "];\n");
}
void smat::dump(const std::string &name) const
{
fdump(stdout, name);
}
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const smat &a)
{
for (int i = 0; i < a.rows(); i++)
os << a[i] << std::endl;
return(os);
}
// **************************************************************************
//
// **************************************************************************
}

85
qxx/src/strprintf.cpp Normal file
View File

@@ -0,0 +1,85 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cstdio>
#include <cstdarg>
#include <gmp.h>
#include "qxx/strprintf.hpp"
namespace q {
// ***************************************************************************
//
// ***************************************************************************
void format(std::string &s, const char *fmt, ...)
{
va_list ap;
int l;
va_start(ap, fmt);
l = gmp_vsnprintf(NULL, 0, fmt, ap);
va_end(ap);
if (l <= 0) {
s = "";
return;
}
char *store = new char[l + 1];
va_start(ap, fmt);
gmp_vsnprintf(store, l + 1, fmt, ap);
va_end(ap);
s = store;
delete[] store;
}
// ***************************************************************************
//
// ***************************************************************************
std::string strprintf(const char *fmt, ...)
{
va_list ap;
int l;
va_start(ap, fmt);
l = gmp_vsnprintf(NULL, 0, fmt, ap);
va_end(ap);
if (l <= 0)
return(std::string());
char *store = new char[l + 1];
va_start(ap, fmt);
gmp_vsnprintf(store, l + 1, fmt, ap);
va_end(ap);
std::string s = store;
delete[] store;
return(s);
}
} // namespace

695
qxx/src/svec.cpp Normal file
View File

@@ -0,0 +1,695 @@
/*
This file is part of qxx -- matrix algebra in exact arithmetic
Copyright (C) 2013-2014 Laurent Poirrier
libp is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with pxx. If not, see <http://www.gnu.org/licenses/>.
*/
#include "qxx/debug.hpp"
#include "qxx/svec.hpp"
#include <algorithm>
// **************************************************************************
//
// **************************************************************************
namespace q {
// **************************************************************************
// svec_ref: construction and destruction
// **************************************************************************
svec_ref::svec_ref(svec &vec, int i)
: vector(vec), index(i)
{
}
svec_ref::svec_ref(const svec_ref &src)
: vector(src.vector), index(src.index)
{
}
svec_ref::~svec_ref()
{
}
// **************************************************************************
// svec_ref: methods
// **************************************************************************
const mpq &svec_ref::get() const
{
return(vector.get(index));
}
// **************************************************************************
// svec_ref: assignment operators
// **************************************************************************
svec_ref::operator mpq() const
{
return(get());
}
const mpq &svec_ref::operator=(const svec_ref &b) const
{
return(vector.set(index, b.get()));
}
const mpq &svec_ref::operator=(const mpq &b) const
{
return(vector.set(index, b));
}
const mpq &svec_ref::operator+=(const mpq &b) const
{
return(vector.set(index, vector.get(index) + b));
}
const mpq &svec_ref::operator-=(const mpq &b) const
{
return(vector.set(index, vector.get(index) - b));
}
const mpq &svec_ref::operator*=(const mpq &b) const
{
return(vector.set(index, vector.get(index) * b));
}
const mpq &svec_ref::operator/=(const mpq &b) const
{
return(vector.set(index, vector.get(index) / b));
}
// **************************************************************************
// svec_ref: read-only operators
// **************************************************************************
const mpq &svec_ref::operator+() const
{
return(get());
}
mpq svec_ref::operator-() const
{
return(-get());
}
mpq svec_ref::operator+(const mpq &b) const
{
return(get() + b);
}
mpq svec_ref::operator-(const mpq &b) const
{
return(get() - b);
}
mpq svec_ref::operator*(const mpq &b) const
{
return(get() * b);
}
mpq svec_ref::operator/(const mpq &b) const
{
return(get() / b);
}
bool svec_ref::operator<(const mpq &b) const
{
return(get() < b);
}
bool svec_ref::operator>(const mpq &b) const
{
return(get() > b);
}
bool svec_ref::operator<=(const mpq &b) const
{
return(get() <= b);
}
bool svec_ref::operator>=(const mpq &b) const
{
return(get() >= b);
}
bool svec_ref::operator==(const mpq &b) const
{
return(get() == b);
}
bool svec_ref::operator!=(const mpq &b) const
{
return(get() != b);
}
// **************************************************************************
// svec_el
// **************************************************************************
svec_el::svec_el()
{
}
svec_el::svec_el(int idx, const mpq &v)
: index(idx), value(v)
{
}
svec_el::svec_el(const svec_el &e)
: index(e.index), value(e.value)
{
}
svec_el::~svec_el()
{
}
// **************************************************************************
// svec
// **************************************************************************
const mpq svec::zero;
svec::svec()
: n(0)
{
}
svec::svec(const svec &src)
: n(src.n), d(src.d)
{
}
svec::svec(int n0)
: n(n0)
{
}
svec::svec(const dvec &src)
: n(src.size())
{
gather(src);
}
svec::~svec()
{
}
// **************************************************************************
//
// **************************************************************************
void svec::clear()
{
n = 0;
d.clear();
}
int svec::size() const
{
return(n);
}
void svec::resize(int n0)
{
n = n0;
}
int svec::nz() const
{
return(d.size());
}
// **************************************************************************
//
// **************************************************************************
bool svec::offs_active() const
{
return(offs.size() != 0);
}
void svec::offs_build()
{
if (offs_active())
return;
offs.resize(n);
for (int i = 0; i < n; i++)
offs[i] = 0;
for (int l = 0; l < nz(); l++)
offs[d[l].index] = l + 1;
}
void svec::offs_clear()
{
offs.resize(0);
}
// **************************************************************************
//
// **************************************************************************
void svec::set_zero()
{
d.clear();
if (offs_active()) {
for (int i = 0; i < n; i++)
offs[i] = 0;
}
}
// **************************************************************************
//
// **************************************************************************
void svec::gather(const dvec &src)
{
n = src.size();
d.clear();
for (int i = 0; i < n; i++)
push(i, src[i]);
}
void svec::spread(dvec &r) const
{
r.resize(n);
r.set(zero);
for (int l = 0; l < nz(); l++)
r[d[l].index] = d[l].value;
}
dvec svec::dense() const
{
dvec r;
spread(r);
return(r);
}
// **************************************************************************
//
// **************************************************************************
int svec::locate(int idx) const
{
Q_RANGE_CHECK(idx, 0, n - 1);
if (offs_active())
return(offs[idx] - 1);
for (int l = 0; l < nz(); l++) {
if (d[l].index == idx)
return(l);
}
return(-1);
}
// **************************************************************************
//
// **************************************************************************
void svec::remove(int l)
{
Q_RANGE_CHECK(l, 0, nz() - 1);
int k = nz() - 1;
if (offs_active())
offs[d[l].index] = 0;
if (l < k)
d[l] = d[k];
d.resize(k);
}
// **************************************************************************
//
// **************************************************************************
int &svec::index(int offset)
{
Q_RANGE_CHECK(offset, 0, nz() - 1);
return(d[offset].index);
}
mpq &svec::value(int offset)
{
Q_RANGE_CHECK(offset, 0, nz() - 1);
return(d[offset].value);
}
int svec::index(int offset) const
{
Q_RANGE_CHECK(offset, 0, nz() - 1);
return(d[offset].index);
}
const mpq &svec::value(int offset) const
{
Q_RANGE_CHECK(offset, 0, nz() - 1);
return(d[offset].value);
}
// **************************************************************************
//
// **************************************************************************
void svec::push_nz(int idx, const mpq &v)
{
Q_RANGE_CHECK(idx, 0, n - 1);
d.push_back(svec_el(idx, v));
if (offs_active())
offs[idx] = nz();
}
void svec::push(int idx, const mpq &v)
{
Q_RANGE_CHECK(idx, 0, n - 1);
if (v.sign())
push_nz(idx, v);
}
const mpq &svec::get(int idx) const
{
Q_RANGE_CHECK(idx, 0, n - 1);
int l = locate(idx);
if (l == -1)
return(zero);
return(d[l].value);
}
const mpq &svec::set(int idx, const mpq &v)
{
Q_RANGE_CHECK(idx, 0, n - 1);
int l = locate(idx);
if (v.sign()) {
if (l == -1) {
push_nz(idx, v);
return(d[nz() - 1].value);
}
d[l].value = v;
return(d[l].value);
} else {
if (l != -1)
remove(l);
return(zero);
}
}
// **************************************************************************
//
// **************************************************************************
const mpq &svec::operator[](int idx) const
{
return(get(idx));
}
svec_ref svec::operator[](int idx)
{
return(svec_ref(*this, idx));
}
// **************************************************************************
//
// **************************************************************************
const svec &svec::operator+() const
{
return(*this);
}
svec svec::operator-() const
{
svec r(*this);
for (int l = 0; l < nz(); l++)
r.d[l].value.set_neg();
return(r);
}
svec svec::operator+(const svec &b) const
{
Q_MATCH(n, b.n);
svec r(*this);
r += b;
return(r);
}
svec svec::operator-(const svec &b) const
{
Q_MATCH(n, b.n);
svec r(*this);
r -= b;
return(r);
}
// **************************************************************************
//
// **************************************************************************
mpq svec::operator*(const svec &b) const
{
Q_MATCH(n, b.n);
mpq v;
if (offs_active()) {
for (int l = 0; l < b.nz(); l++) {
int i = b.d[l].index;
if ((i < n) && (offs[i]))
v += d[offs[i] - 1].value * b.d[l].value;
}
} else if (b.offs_active()) {
for (int l = 0; l < nz(); l++) {
int i = d[l].index;
if ((i < b.n) && (b.offs[i]))
v += d[l].value * b.d[b.offs[i] - 1].value;
}
} else {
dvec db;
b.spread(db);
v = (*this) * db;
}
return(v);
}
mpq svec::operator*(const dvec &b) const
{
Q_MATCH(n, b.size());
mpq v;
for (int l = 0; l < nz(); l++) {
int j = d[l].index;
if (j < b.size())
v += d[l].value * b[j];
}
return(v);
}
// **************************************************************************
//
// **************************************************************************
svec svec::operator*(const mpq &v)
{
svec r(n);
if (!v.sign())
return(r);
for (int l = 0; l < nz(); l++)
r.push(d[l].index, d[l].value * v);
return(r);
}
svec svec::operator/(const mpq &v)
{
svec r(n);
for (int l = 0; l < nz(); l++)
r.push(d[l].index, d[l].value / v);
return(r);
}
// **************************************************************************
//
// **************************************************************************
svec &svec::operator=(const svec &b)
{
n = b.n;
d = b.d;
if ((offs_active()) && (b.offs_active()))
offs = b.offs;
else
offs_clear();
return(*this);
}
svec &svec::operator+=(const svec &b)
{
Q_MATCH(n, b.n);
if (offs_active()) {
for (int l = 0; l < b.nz(); l++) {
int i = b.d[l].index;
if ((i < n) && (offs[i]))
d[offs[i] - 1].value += b.d[l].value;
else
push(i, b.d[l].value);
}
} else {
dvec da, db;
spread(da);
b.spread(db);
da += db;
gather(da);
}
return(*this);
}
svec &svec::operator-=(const svec &b)
{
Q_MATCH(n, b.n);
if (offs_active()) {
for (int l = 0; l < b.nz(); l++) {
int i = b.d[l].index;
if ((i < n) && (offs[i]))
d[offs[i] - 1].value -= b.d[l].value;
else
push(i, -b.d[l].value);
}
} else {
dvec da, db;
spread(da);
b.spread(db);
da += db;
gather(da);
}
return(*this);
}
svec &svec::operator*=(const mpq &v)
{
if (!v.sign()) {
set_zero();
return(*this);
}
for (int l = 0; l < nz(); l++)
d[l].value *= v;
return(*this);
}
svec &svec::operator/=(const mpq &v)
{
for (int l = 0; l < nz(); l++)
d[l].value /= v;
return(*this);
}
// **************************************************************************
//
// **************************************************************************
static bool comp(const svec_el &a, const svec_el &b)
{
return(a.index < b.index);
}
void svec::sort()
{
std::sort(d.begin(), d.end(), comp);
}
// **************************************************************************
//
// **************************************************************************
void svec::fdump(FILE *f, const std::string &name) const
{
if (name.length())
gmp_fprintf(f, "%s = ", name.c_str());
for (int l = 0; l < nz(); l++) {
gmp_fprintf(f, " %d:%Zd/%Zd",
d[l].index,
mpq_numref(d[l].value.v),
mpq_denref(d[l].value.v)
);
}
gmp_fprintf(f, "\n");
}
void svec::dump(const std::string &name) const
{
fdump(stdout, name);
}
// **************************************************************************
//
// **************************************************************************
std::ostream &operator<<(std::ostream &os, const svec &v)
{
for (int l = 0; l < v.nz(); l++)
os << " " << v.index(l) << ":" << v.value(l);
return(os);
}
// **************************************************************************
//
// **************************************************************************
}