linbox
examples/samplebb.C

generate an example matrix with specified frobenius form.

generate an example matrix with specified frobenius form.samplebb takes options and any number of argument triples denoting companion matrix blocks. For example, the call "samplebb -r 7 2 3 a3 1 1" generates a sparsely randomized matrix (because of the '-r' option) matrix which is similar (because of the two triples '7 2 3' and 'a2 1 1') to a direct sum of 3 companion matrices for (x-7)^2 plus one companion matrix for x^3 + x + 2, the polynomial denoted by 'a3'.

In general, in the first position of each triple 'aK' denotes the polynomial x^k + x + K-1 and a number n denotes the polynomial x-n. The second number in the triple specifies a power of the polynomial and the third specifies how many companion matrix blocks for that power of that polynomial.

Possible options are -r lightly randomized similarity transform, matrix remains sparse. -R fully randomized similarity transform, matrix becomes dense.

The matrix is written to standard out in SMS format (triples).

For some other examples: "samplebb 1 1 2 2 1 2 4 1 2 12 1 1 0 1 1" is a 8 by 8 diagonal matrix in smith form, diag(1,1,2,2,4,4,12,0)

/*
* examples/samplebb.C
*
* Copyright (C) 2005, 2010 D Saunders
* ========LICENCE========
* This file is part of the library LinBox.
*
* LinBox is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library 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
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
* ========LICENCE========
*/
#include <iostream>
#include <string>
#include <list>
#include <vector>
#include <linbox/blackbox/direct-sum.h>
#include <linbox/blackbox/companion.h>
#include <linbox/ring/ntl.h>
#include <NTL/ZZX.h>
using std::string;
using std::list;
using NTL::ZZX;
void stripOptions(int& acp, char* avp[], string& opts, const int ac, char** av)
{
acp = 0;
for (int i = 1; i < ac; ++i)
{
//std::cout << av[i] << " ";
if (av[i][0] == '-') {
for (const char* j = av[i]+1; *j != 0; ++j)
opts.push_back(*j);
}
else {
avp[acp] = av[i];
++acp;
}
}
}
template <class List, class Ring>
void augmentBB(List& L, char* code, int e, int k, const Ring& R)
{
typedef typename Ring::Element Int;
Int a;
ZZX p;
// build poly p
if ( *code != 'a') // build linear poly
{
R.init(a, -atoi(code));
p += ZZX(0, a);
p += ZZX(1, R.one);
}
else // build long poly
{
int n = atoi(code+1);
R.init(a, n-1);
p += ZZX(n, R.one);
p += ZZX(1, R.one);
p += ZZX(0, a);
}
//std::cout << "(code, e, k) =(" << code << ", " << e << ", " << k << ")" << std::endl;
//std::cout << "Correspoding poly: " << p << std::endl;
// compute q = p^e
ZZX q(0, R.one);
for(int i = 0; i < e; ++i) q *= p;
//std::cout <<"Polynomial: " << q << std::endl;
LinBox::DenseVector<Ring> v(R,deg(q)+1);
for (int i = 0; i < v.size(); ++i) v[i] = coeff(q, i);
// companion matrix of q
Companion<Ring>* C = new Companion<Ring>(R, v);
for(int i = 0; i < k; ++i) L.push_back(C);
}
template < class Ring >
void scramble(DenseMatrix<Ring>& M)
{
Ring R = M.field();
int N,n = M.rowdim(); // number of random basic row and col ops.
N = 2*n;
for (int k = 0; k < N; ++k) {
int i = rand()%M.rowdim();
int j = rand()%M.coldim();
if (i == j) continue;
// M*i += alpha M*j and Mj* -= alpha Mi*
typename Ring::Element alpha, beta, x;
R.init(alpha, rand()%5 - 2);
R.neg(beta, alpha);
for (size_t l = 0; l < M.rowdim(); ++l) if (!R.isZero(alpha)) {
R.mul(x, alpha, M[l][j]);
R.addin(M[l][i], x);
}
for (size_t l = 0; l < M.rowdim(); ++l) if (!R.isZero(alpha)) {
R.mul(x, beta, M[i][l]);
R.addin(M[j][l], x);
}
}
/*
std::ofstream out("matrix", std::ios::out);
//M. write(std::cout);
out << n << " " << n << "\n";
for (int i = 0; i < n; ++ i) {
for ( int j = 0; j < n; ++ j) {
R. write(out, M[i][j]);
out << " ";
}
out << "\n";
}
*/
}
template <class Matrix>
void printMatrix (const Matrix& A)
{
int m = A. rowdim();
int n = A. coldim();
typedef typename Matrix::Field Ring;
typedef typename Ring::Element Element;
const Ring &r = A.field();
LinBox::DenseVector<Ring> x(r,m), y(r,n);
std::cout << m << " " << n << " M" << std::endl;
typename LinBox::DenseVector<Ring>::iterator y_p;
for (int i = 0; i < m; ++ i) {
r. assign (x[i], r.one);
A. applyTranspose(y, x);
for(y_p = y. begin(); y_p != y. end(); ++ y_p)
if (! r.isZero(*y_p))
std::cout << i+1 << " " << y_p - y.begin() + 1 << " " << *y_p << std::endl;
r. assign (x[i], r.zero);
}
std::cout << "0 0 0" << std::endl;
}
int main(int ac, char* av[])
{
if (ac < 2)
{ std::cout << "usage: " << av[0] <<
" options block-groups." << std::endl;
std::cout << av[0] << " -r 1 2 3 a4 1 1" << std::endl;
std::cout <<
"for lightly randomized matrix similar to direct sum of 3 copies of companion " << std::endl
<< "matrix of (x-1)^2 and one copy of companion matrix of (x^4 + x + 3)^1." << std::endl;
}
typedef NTL_ZZ Ring;
Ring Z;
typedef Companion<Ring> BB;
int acp; char* avp[ac];
string opts;
stripOptions(acp, avp, opts, ac, av);
//std::cout << "number of triples: " << acp << std::endl;
//for (int i = 0; i < acp; ++ i)
// std::cout << avp[i];
//std::cout << std::endl;
//std::cout << "Begin to ....\n";
list<BB*> L;
for (int i = 0; i < acp; i += 3)
augmentBB(L, avp[i], atoi(avp[i+1]), atoi(avp[i+2]), Z);
DirectSum<BB> A(L);
//std::cout <<"Option: " << opts.c_str() << std::endl;
if (opts.size() >= 1)
{ if (opts[0] == 'r')
{
// into sparse matrix, then 3n row ops with corresponding col ops
DenseMatrix<Ring> B(Z,A.rowdim(), A.coldim());
//MatrixDomain<Ring> MD(Z);
LinBox::MatrixHom::map (B, A);
scramble(B);
printMatrix(B);
// delete B;
}
if (opts[0] == 'R') { throw LinBox::NotImplementedYet(); }
// into dense matrix, then many row ops
//...
}
else {
printMatrix (A);
}
return 0 ;
}
// Local Variables:
// mode: C++
// tab-width: 4
// indent-tabs-mode: nil
// c-basic-offset: 4
// End:
// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
A BlasVector<_Field > represents a vector as an array of _Field::Elements.
If C = DirectSum(A, B) and y = xA and z = wB, then (y,z) = (x,w)C.
Definition: direct-sum.h:58
the integer ring.
Definition: ring/ntl/ntl-zz.h:58
Ring (in fact, a unique factorization domain) of polynomial with coefficients in class NTL_zz_p (inte...
Definition: ntl-lzz_px.h:78
bool isZero(const Element &x) const
Test if an element equals zero.
Definition: ntl-lzz_px.h:234
int main()
no command line args
Definition: ex-fields-archetype.C:70
Matrix Homomorphism A map function converts a matrix on a field/ring to its natural image in another ...
Companion matrix of a monic polynomial.
Definition: companion.h:41