example matrices that were chosen for Smith form testing.
example matrices that were chosen for Smith form testing.
Various Smith form algorithms may be used for matrices over the integers or over Z_m. Moduli greater than 2^32 are not supported here. Several types of example matrices may be constructed or the matrix be read from a 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 computation 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 integer 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 <linbox/ring/pir-ntl-zz_p.h>
template <class PIR>
int main(
int argc,
char* argv[])
{
if (argc < 3 or argc > 4) {
cout << "usage: " << argv[0] << " type n [filename]" << endl;
cout << " type = `random', `random-rough', `tref', or `fib',"
<< " and n is the dimension" << endl;
cout << " If filename is present, matrix is written there, else to cout." << endl;
return 0;
}
string type = argv[1];
int n = atoi(argv[2]);
typedef Givaro::ZRing<Integer> PIR;
PIR R;
Mat(M,R,n,type);
if (argc == 4) {
ofstream out(argv[3]);
M.write(out) << endl;
} else {
M.write(cout) << endl;
}
}
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]);
}
}
}
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, cin);
R. init (M[(size_t)i][(size_t)j], Val);
}
}
template <class PIR>
string src) {
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 if (src == "krat") KratMat(M, R, n);
else {
}
}
const _Field & field() const
Retrieve a reference to a row.
Definition blas-matrix.inl:1535
size_t rowdim() const
Get the number of rows in the matrix.
Definition blas-matrix.inl:502
void resize(const size_t &m, const size_t &n, const Element &val=Element())
Resize the matrix to the given dimensions.
Definition blas-matrix.inl:526
size_t coldim() const
Get the number of columns in the matrix.
Definition blas-matrix.inl:508
std::istream & read(std::istream &file)
Read the matrix from an input stream.
Definition blas-matrix.inl:423
Ring (in fact, a unique factorization domain) of polynomial with coefficients in class NTL_zz_p (inte...
Definition ntl-lzz_px.h:78
int main()
no command line args
Definition ex-fields-archetype.C:70
Givaro::Integer integer
Integers in LinBox.
Definition integer.h:55
linbox base configuration file
Namespace in which all linbox code resides.
Definition alt-blackbox-block-container.h:4
LinBox timer is Givaro's.