Various Smith form algorithms may be used for matrices over the integers or over Z_m. Moduli greater than 2^32 are not supported. Several types of example matrices may be constructed or matrix read from file. Run the program with no arguments for a synopsis of the command line parameters.
For the "adaptive" method, the matrix must be over the integers. This is expected to work best for large matrices.
For the "2local" method, the computaattion is done mod 2^32.
For the "local" method, the modulus must be a prime power.
For the "ilio" method, the modulus may be arbitrary composite. If the modulus is a multiple of the integer determinant, the intege Smith form is obtained. Determinant plus ilio may be best for smaller matrices.
This example was used during the design process of the adaptive algorithm.
#include <iostream>
#include <string>
#include <vector>
#include <list>
#include <linbox/algorithms/smith-form-sparseelim-local.h>
#include <linbox/ring/local2_32.h>
#include <linbox/ring/pir-modular-int32.h>
#include <linbox/algorithms/smith-form-local.h>
#include <linbox/algorithms/smith-form-iliopoulos.h>
template<class PIR>
string src, string file, string format);
template<
class I1,
class Lp>
void distinct (I1 a, I1 b, Lp& c);
template <class I> void display(I b, I e);
int main(int argc, char* argv[])
{
typedef PIRModular<int32_t> PIR;
if (argc < 5) {
cout << "usage: " << argv[0] << " alg m n source format \n" << endl;
cout << "alg = `adaptive', `ilio', `local', or `2local', \n"
<< "m is modulus (ignored by 2local, adaptive), "
<< "n is matrix order, \n"
<< "source is `random', `random-rough', `fib', `tref', or a filename \n"
<< "format is `dense' or `sparse' (if matrix from a file)\n"
<< "compile with -DBIG if you want big integers used.\n";
return 0;
}
string algo = argv[1];
unsigned long m = atoi(argv[2]);
int n = atoi(argv[3]);
string src = argv[4];
string file = src;
string format = (argc >= 6 ? argv[5] : "");
UserTimer T;
if (algo == "adaptive")
{
typedef Givaro::ZRing<Integer> Ints;
Ints Z;
std::ifstream input (file);
M.read(input);
DenseVector<Givaro::ZRing<Integer> > v(Z,(size_t)n);
T.start();
T.stop();
list<pair<integer, size_t> > p;
cout << "#";
display(p.begin(), p.end());
cout << "# adaptive, Ints, n = " << n << endl;
cout << "T" << n << "adaptive" << m << " := ";
}
else if (algo == "ilio") {
PIR R( (int32_t)m);
Mat(M, R, n, src, file, format);
T.start();
SmithFormIliopoulos::smithFormIn (M);
T.stop();
typedef list< PIR::Element > List;
List L;
for (size_t i = 0; i < M.rowdim(); ++i)
L.push_back(M[(size_t)i][(size_t)i]);
list<pair<PIR::Element, size_t> > p;
cout << "#";
display(p.begin(), p.end());
cout << "# ilio, PIR-Modular-int32_t(" << m << "), n = " << n << endl;
cout << "T" << n << "ilio" << m << " := ";
}
else if (algo == "local") {
if (format == "sparse" ) {
typedef Givaro::Modular<int32_t> Field;
Field F(m);
std::ifstream input (argv[4]);
if (!input) { std::cerr << "Error opening matrix file: " << argv[1] << std::endl; return -1; }
SparseMatrix<Field, SparseMatrixFormat::SparseSeq > B (ms);
std::cout << "B is " << B.rowdim() << " by " << B.coldim() << std::endl;
if (B.rowdim() <= 20 && B.coldim() <= 20) B.write(std::cout) << std::endl;
Integer p(m), im(m);
Givaro::IntPrimeDom IPD;
for(unsigned int k = 2; ( ( ! IPD.isprime(p) ) && (p > 1) ); ++k)
Givaro::root( p, im, k );
std::vector<std::pair<size_t,Field::Element> > local;
PGD(local, B, (int32_t)m, (int32_t)p);
typedef list< Field::Element > List;
List L;
for ( auto p_it = local.begin(); p_it != local.end(); ++p_it) {
for(size_t i = 0; i < (size_t) p_it->first; ++i)
L.push_back((Field::Element)p_it->second);
}
size_t M = (B.rowdim() > B.coldim() ? B.coldim() : B.rowdim());
for (size_t i = L.size(); i < M; ++i)
L.push_back(0);
list<pair<Field::Element, size_t> > pl;
std::cout << "#";
display(pl.begin(), pl.end());
cout << "# local, PowerGaussDomain<int32_t>(" << M << "), n = " << n << endl;
}
else {
PIR R( (int32_t)m);
Mat(M, R, n, src, file, format);
typedef list< PIR::Element > List;
List L;
T.start();
SmithForm( L, M, R );
T.stop();
list<pair<PIR::Element, size_t> > p;
cout << "#";
display(p.begin(), p.end());
cout << "# local, PIR-Modular-int32_t(" << m << "), n = " << n << endl;
}
cout << "T" << n << "local" << m << " := ";
}
else if (algo == "2local") {
Mat(M, R, n, src, file, format);
typedef list< Local2_32::Element > List;
List L;
T.start();
SmithForm( L, M, R );
T.stop();
list<pair<Local2_32::Element, size_t> > p;
cout << "#";
display(p.begin(), p.end());
cout << "# 2local, Local2_32, n = " << n << endl;
cout << "T" << n << "local2_32 := ";
}
else {
printf ("Unknown algorithms\n");
exit (-1);
}
T.print(cout); cout << ";" << endl;
return 0 ;
}
template < class Ring >
{
N = n;
for (int k = 0; k < N; ++k) {
int i = rand()%(int)M.
rowdim();
int j = rand()%(int)M.
coldim();
if (i == j) continue;
int a = 0;
for (
size_t l = 0; l < M.
rowdim(); ++l) {
if (a)
R.subin(M[(size_t)l][(size_t)i], M[(size_t)l][(size_t)j]);
else
R.addin(M[(size_t)l][(size_t)i], M[(size_t)l][(size_t)j]);
}
for (
size_t l = 0; l < M.
coldim(); ++l) {
if (a)
R.subin(M[(size_t)i][l], M[(size_t)j][l]);
else
R.addin(M[(size_t)i][l], M[(size_t)j][l]);
}
}
std::ofstream out("matrix", std::ios::out);
out << n << " " << n << "\n";
for (int i = 0; i < n; ++ i) {
for ( int j = 0; j < n; ++ j) {
R. write(out, M[(size_t)i][(size_t)j]);
out << " ";
}
out << "\n";
}
}
template <class PIR>
M.
resize((
size_t)n, (
size_t)n, R.zero);
if (n > 10000) {cerr << "n too big" << endl; exit(-1);}
int jth_factor[130] =
{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149,
151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313,
317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691,
701, 709, 719, 727, 733};
for (int j= 0, i = 0 ; i < n; ++j)
{
typename PIR::Element v; R.init(v, jth_factor[25+j]);
for (int k = j ; k > 0 && i < n ; --k)
{ M[(size_t)i][(size_t)i] = v; ++i;
if (i < n) {M[(size_t)i][(size_t)i] = v; ++i;}
}
}
scramble(M);
}
template <class PIR>
M.
resize((
size_t)n,(
size_t) n, R.zero);
for (int i= 0 ; i < n; ++i)
R.init(M[(size_t)i][(size_t)i], i % 1000 + 1);
scramble(M);
}
template <class PIR>
M.
resize((
size_t)n,(
size_t) n, R.zero);
typename PIR::Element pmone; R.assign(pmone, R.one);
for (int i= 0 ; i < n; ++i) M[(size_t)i][(size_t)i] = R.one;
int j = 1, k = 0;
for (int i= 0 ; i < n-1; ++i) {
if ( i == k) {
M[(size_t)i][(
size_t)i+1] = R.
zero;
k += ++j;
}
else {
M[(size_t)i][(size_t)i+1] = pmone;
R.negin(pmone);
}
R.neg(M[(size_t)i+1][(size_t)i], M[(size_t)i][(size_t)i+1]);
}
scramble(M);
}
template <class PIR>
M.
resize((
size_t)n, (
size_t)n, R.zero);
std::vector<int> power2;
int i = 1;
do {
power2. push_back(i);
i *= 2;
} while (i < n);
std::ifstream in ("prime", std::ios::in);
for ( i = 0; i < n; ++ i)
in >> M[(size_t)i][(size_t)i];
std::vector<int>::iterator p;
for ( i = 0; i < n; ++ i) {
for ( p = power2. begin(); (p != power2. end()) && (*p <= i); ++ p)
M[(size_t)i][(size_t)(i - *p)] = 1;
for ( p = power2. begin(); (p != power2. end()) && (*p < n - i); ++ p)
M[(size_t)i][(size_t)(i + *p)] = 1;
}
}
struct pwrlist
{
vector<integer> m;
{ m.push_back(1); m.push_back(q);
}
{
for (int i = (int)m.size(); i <= e; ++i) m.push_back(m[1]*m[(size_t)i-1]);
return m[(size_t)e];
}
};
template <class num>
num& qread(num& Val, pwrlist& M, istream& in)
{
char c;
in >> c;
if (c == '0') return Val = 0;
if (c == '1') return Val = 1;
if (c != 'p' && c != 'q') { cout << "exiting due to unknown char " << c << endl; exit(-1);}
in.get(c);
if (c !='^') {in.putback(c); return Val = M[1];}
else
{ int expt; in >> expt;
return Val = M[expt];
};
}
template <class PIR>
{
pwrlist pwrs(q);
for (
unsigned int i = 0; i < M.
rowdim(); ++ i)
for (
unsigned int j = 0; j < M.
coldim(); ++ j) {
int Val;
qread(Val, pwrs, in);
R. init (M[(size_t)i][(size_t)j], Val);
}
}
template<class I1, class Lp>
{
typename iterator_traits<I1>::value_type e;
size_t count = 0;
if (a != b) {e = *a; ++a; count = 1;}
else return;
while (a != b)
{ if (*a == e) ++count;
else
{ c.push_back(typename Lp::value_type(e, count));
e = *a; count = 1;
}
++a;
}
c.push_back(typename Lp::value_type(e, count));
return;
}
template <class I>
void display(I b, I e)
{ cout << "(";
for (I p = b; p != e; ++p) cout << p->first << " " << p->second << ", ";
cout << ")" << endl;
}
template <class PIR>
string src, string file, string format) {
if (src == "random-rough") RandomRoughMat(M, R, n);
else if (src == "random") RandomFromDiagMat(M, R, n);
else if (src == "fib") RandomFibMat(M, R, n);
else if (src == "tref") TrefMat(M, R, n);
else
{
int rdim, cdim;
std::ifstream in (file.c_str(), std::ios::in);
if (! in) { cerr << "error: unable to open file" << endl; exit(-1); }
in >> rdim >> cdim;
M. resize ((size_t)rdim, (size_t)cdim);
if (format == "dense" ) {
for (int i = 0; i < rdim; ++ i)
for ( int j = 0; j < cdim; ++ j) {
in >> Val;
R. init (M[(size_t)i][(size_t)j], Val);
}
}
else if (format == "sparse") {
int i, j;
char mark;
in >> mark;
do {
in >> i >> j;
in. ignore (1);
in >> val;
if ( i == 0) break;
R. init (M[(size_t)i-1][(size_t)j-1], val);
} while (true);
}
else if (format == "kdense") KratMat(M, R, n, in);
else {
cout << "Format: " << format << " Unknown\n";
exit (-1);
}
}
}