<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<title>linbox: examples/smith.C</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css"/>
</head>
<body>
<!-- Generated by Doxygen 1.7.4 -->
<script type="text/javascript">
function hasClass(ele,cls) {
return ele.className.match(new RegExp('(\\s|^)'+cls+'(\\s|$)'));
}
function addClass(ele,cls) {
if (!this.hasClass(ele,cls)) ele.className += " "+cls;
}
function removeClass(ele,cls) {
if (hasClass(ele,cls)) {
var reg = new RegExp('(\\s|^)'+cls+'(\\s|$)');
ele.className=ele.className.replace(reg,' ');
}
}
function toggleVisibility(linkObj) {
var base = linkObj.getAttribute('id');
var summary = document.getElementById(base + '-summary');
var content = document.getElementById(base + '-content');
var trigger = document.getElementById(base + '-trigger');
if ( hasClass(linkObj,'closed') ) {
summary.style.display = 'none';
content.style.display = 'block';
trigger.src = 'open.png';
removeClass(linkObj,'closed');
addClass(linkObj,'opened');
} else if ( hasClass(linkObj,'opened') ) {
summary.style.display = 'block';
content.style.display = 'none';
trigger.src = 'closed.png';
removeClass(linkObj,'opened');
addClass(linkObj,'closed');
}
return false;
}
</script>
<div id="top">
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
<tbody>
<tr style="height: 56px;">
<td style="padding-left: 0.5em;">
<div id="projectname">linbox</div>
</td>
</tr>
</tbody>
</table>
</div>
<div id="navrow1" class="tabs">
<ul class="tablist">
<li><a href="index.html"><span>Main Page</span></a></li>
<li><a href="pages.html"><span>Related Pages</span></a></li>
<li><a href="modules.html"><span>Modules</span></a></li>
<li><a href="namespaces.html"><span>Namespaces</span></a></li>
<li><a href="annotated.html"><span>Data Structures</span></a></li>
<li><a href="files.html"><span>Files</span></a></li>
<li><a href="dirs.html"><span>Directories</span></a></li>
<li><a href="examples.html"><span>Examples</span></a></li>
</ul>
</div>
</div>
<div class="header">
<div class="headertitle">
<div class="title">examples/smith.C</div> </div>
</div>
<div class="contents">
<p>mod m Smith form by elmination</p>
<dl class="author"><dt><b>Author:</b></dt><dd>bds & zw</dd></dl>
<p>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.</p>
<p>For the "adaptive" method, the matrix must be over the integers. This is expected to work best for large matrices.</p>
<p>For the "2local" method, the computaattion is done mod 2^32.</p>
<p>For the "local" method, the modulus must be a prime power.</p>
<p>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.</p>
<p>This example was used during the design process of the adaptive algorithm.</p>
<div class="fragment"><pre class="fragment"><span class="comment">/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */</span>
<span class="comment">// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s</span>
<span class="comment">/*</span>
<span class="comment"> * examples/smith.C</span>
<span class="comment"> *</span>
<span class="comment"> * Copyright (C) 2005, 2010 D. Saunders, Z. Wang, J-G Dumas</span>
<span class="comment"> *</span>
<span class="comment"> * This file is part of LinBox.</span>
<span class="comment"> *</span>
<span class="comment"> * LinBox is free software: you can redistribute it and/or modify</span>
<span class="comment"> * it under the terms of the GNU Lesser General Public License as</span>
<span class="comment"> * published by the Free Software Foundation, either version 2 of</span>
<span class="comment"> * the License, or (at your option) any later version.</span>
<span class="comment"> *</span>
<span class="comment"> * LinBox is distributed in the hope that it will be useful,</span>
<span class="comment"> * but WITHOUT ANY WARRANTY; without even the implied warranty of</span>
<span class="comment"> * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the</span>
<span class="comment"> * GNU Lesser General Public License for more details.</span>
<span class="comment"> *</span>
<span class="comment"> * You should have received a copy of the GNU Lesser General Public</span>
<span class="comment"> * License along with LinBox. If not, see</span>
<span class="comment"> * <http://www.gnu.org/licenses/>.</span>
<span class="comment"> */</span>
<span class="preprocessor">#include <iostream></span>
<span class="preprocessor">#include <string></span>
<span class="preprocessor">#include <vector></span>
<span class="preprocessor">#include <list></span>
<span class="keyword">using namespace </span>std;
<span class="preprocessor">#include "linbox/field/modular-int32.h"</span>
<span class="preprocessor">#include "linbox/blackbox/sparse.h"</span>
<span class="preprocessor">#include "linbox/algorithms/smith-form-sparseelim-local.h"</span>
<span class="preprocessor">#include "<a class="code" href="timer_8h.html" title="LinBox timer is Givaro's.">linbox/util/timer.h</a>"</span>
<span class="preprocessor">#include "linbox/field/unparametric.h"</span>
<span class="preprocessor">#include "linbox/field/local2_32.h"</span>
<span class="preprocessor">#include "linbox/field/PIR-modular-int32.h"</span>
<span class="preprocessor">#include "linbox/algorithms/smith-form-local.h"</span>
<span class="preprocessor">#include "linbox/algorithms/smith-form-local2.h"</span>
<span class="preprocessor">#include <linbox/algorithms/smith-form-iliopoulos.h></span>
<span class="preprocessor">#include "<a class="code" href="smith-form-adaptive_8h.html" title="Implement the adaptive algorithm for Smith form computation.">linbox/algorithms/smith-form-adaptive.h</a>"</span>
<span class="preprocessor">#include "<a class="code" href="blackbox_2dense_8h.html" title="NO DOC.">linbox/blackbox/dense.h</a>"</span>
<span class="keyword">using namespace </span>LinBox;
<span class="comment">// #ifndef BIG</span>
<span class="keyword">template</span><<span class="keyword">class</span> PIR>
<span class="keywordtype">void</span> <a name="a0"></a><a class="code" href="smith_8_c.html#ac924a785b0083e0f26562df322262428" title="Output matrix is determined by src which may be: "random-rough" This mat will have s...">Mat</a>(<a name="_a1"></a><a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<PIR></a>& M, PIR& R, <span class="keywordtype">int</span> n,
<span class="keywordtype">string</span> src, <span class="keywordtype">string</span> file, <span class="keywordtype">string</span> format);
<span class="keyword">template</span><<span class="keyword">class</span> I1, <span class="keyword">class</span> Lp> <span class="keywordtype">void</span> distinct (I1 a, I1 b, Lp& c);
<span class="keyword">template</span> <<span class="keyword">class</span> I> <span class="keywordtype">void</span> display(I b, I e);
<span class="keywordtype">int</span> main(<span class="keywordtype">int</span> argc, <span class="keywordtype">char</span>* argv[])
{
<span class="keyword">typedef</span> PIRModular<int32_t> PIR;
<span class="keywordflow">if</span> (argc < 5) {
cout << <span class="stringliteral">"usage: "</span> << argv[0] << <span class="stringliteral">" alg m n source format \n"</span> << endl;
cout << <span class="stringliteral">"alg = `adaptive', `ilio', `local', or `2local', \n"</span>
<< <span class="stringliteral">"m is modulus (ignored by 2local, adaptive), "</span>
<< <span class="stringliteral">"n is matrix order, \n"</span>
<< <span class="stringliteral">"source is `random', `random-rough', `fib', `tref', or a filename \n"</span>
<< <span class="stringliteral">"format is `dense' or `sparse' (if matrix from a file)\n"</span>
<< <span class="stringliteral">"compile with -DBIG if you want big integers used.\n"</span>;
<span class="keywordflow">return</span> 0;
}
<span class="keywordtype">string</span> algo = argv[1];
<span class="keywordtype">int</span> m = atoi(argv[2]);
<span class="keywordtype">int</span> n = atoi(argv[3]);
<span class="keywordtype">string</span> src = argv[4];
<span class="keywordtype">string</span> file = src;
<span class="keywordtype">string</span> format = (argc >= 6 ? argv[5] : <span class="stringliteral">""</span>);
UserTimer T;
<span class="keywordflow">if</span> (algo == <span class="stringliteral">"adaptive"</span>)
{
<span class="keyword">typedef</span> <a name="_a2"></a><a class="code" href="class_lin_box_1_1_p_i_d__integer.html" title="Domain for integer operations.">PID_integer</a> Ints;
Ints Z;
<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<Ints></a> M(Z);
<a class="code" href="smith_8_c.html#ac924a785b0083e0f26562df322262428" title="Output matrix is determined by src which may be: "random-rough" This mat will have s...">Mat</a>(M, Z, n, src, file, format);
vector<integer> v(n);
T.start();
<a name="a3"></a><a class="code" href="group__solutions.html#gadebecfaca29a9dbb684a65b9e1e4bc72" title="Compute the Smith form of A.">SmithFormAdaptive::smithForm</a>(v, M);
T.stop();
list<pair<integer, size_t> > p;
distinct(v.begin(), v.end(), p);
cout << <span class="stringliteral">"#"</span>;
display(p.begin(), p.end());
cout << <span class="stringliteral">"# adaptive, Ints, n = "</span> << n << endl;
cout << <span class="stringliteral">"T"</span> << n << <span class="stringliteral">"adaptive"</span> << m << <span class="stringliteral">" := "</span>;
}
<span class="keywordflow">else</span> <span class="keywordflow">if</span> (algo == <span class="stringliteral">"ilio"</span>) {
PIR R(m);
<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<PIR></a> M(R);
<a class="code" href="smith_8_c.html#ac924a785b0083e0f26562df322262428" title="Output matrix is determined by src which may be: "random-rough" This mat will have s...">Mat</a>(M, R, n, src, file, format);
T.start();
SmithFormIliopoulos::smithFormIn (M);
T.stop();
<span class="keyword">typedef</span> list< PIR::Element > List;
List L;
<span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> i = 0; i < M.rowdim(); ++i) L.push_back(M[i][i]);
list<pair<PIR::Element, size_t> > p;
distinct(L.begin(), L.end(), p);
cout << <span class="stringliteral">"#"</span>;
display(p.begin(), p.end());
cout << <span class="stringliteral">"# ilio, PIR-Modular-int32_t("</span> << m << <span class="stringliteral">"), n = "</span> << n << endl;
cout << <span class="stringliteral">"T"</span> << n << <span class="stringliteral">"ilio"</span> << m << <span class="stringliteral">" := "</span>;
}
<span class="keywordflow">else</span> <span class="keywordflow">if</span> (algo == <span class="stringliteral">"local"</span>) { <span class="comment">// m must be a prime power</span>
<span class="keywordflow">if</span> (format == <span class="stringliteral">"sparse"</span> ) {
<span class="keyword">typedef</span> <a name="_a4"></a><a class="code" href="class_lin_box_1_1_modular_3_01int32__t_01_4.html" title="Specialization of Modular to int32_t element type with efficient dot product.">Modular<int32_t></a> <a name="_a5"></a><a class="code" href="class_lin_box_1_1_modular_3_01uint32__t_01_4.html" title="Specialization of class Modular for uint32_t element type.">Field</a>;
Field F(m);
std::ifstream input (argv[4]);
<span class="keywordflow">if</span> (!input) { std::cerr << <span class="stringliteral">"Error opening matrix file: "</span> << argv[1] << std::endl; <span class="keywordflow">return</span> -1; }
<a name="_a6"></a><a class="code" href="class_lin_box_1_1_matrix_stream.html" title="MatrixStream.">MatrixStream<Field></a> ms( F, input );
<a name="_a7"></a><a class="code" href="class_lin_box_1_1_sparse_matrix.html" title="vector of sparse rows.">SparseMatrix<Field, Vector<Field>::SparseSeq</a> > B (ms);
std::cout << <span class="stringliteral">"B is "</span> << B.rowdim() << <span class="stringliteral">" by "</span> << B.coldim() << std::endl;
<span class="keywordflow">if</span> (B.rowdim() <= 20 && B.coldim() <= 20) B.write(std::cout) << std::endl;
Integer p(m), im(m);
<span class="comment">// Should better ask user to give the prime !!!</span>
<span class="keywordflow">for</span>(<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> k = 2; ( ( ! ::Givaro::probab_prime(p) ) && (p > 1) ); ++k)
::Givaro::root( p, im, k );
<span class="comment">// using Sparse Elimination</span>
<a name="_a8"></a><a class="code" href="class_lin_box_1_1_power_gauss_domain.html" title="Repository of functions for rank modulo a prime power by elimination on sparse matrices.">LinBox::PowerGaussDomain< Field ></a> PGD( F );
std::vector<std::pair<size_t,size_t> > local;
PGD(local, B, m, (<span class="keywordtype">int</span>)p);
<span class="keyword">typedef</span> list< Field::Element > List;
List L;
<span class="keywordflow">for</span> ( std::vector<std::pair<size_t,size_t> >::iterator
p_it = local.begin(); p_it != local.end(); ++p_it) {
<span class="keywordflow">for</span>(<span class="keywordtype">size_t</span> i = 0; i < (size_t) p_it->first; ++i)
L.push_back((Field::Element)p_it->second);
}
<span class="keywordtype">size_t</span> M = (B.rowdim() > B.coldim() ? B.coldim() : B.rowdim());
<span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> i = L.size(); i < M; ++i)
L.push_back(0);
list<pair<Field::Element, size_t> > pl;
distinct(L.begin(), L.end(), pl);
std::cout << <span class="stringliteral">"#"</span>;
<span class="comment">//display(local.begin(), local.end());</span>
display(pl.begin(), pl.end());
cout << <span class="stringliteral">"# local, PowerGaussDomain<int32_t>("</span> << M << <span class="stringliteral">"), n = "</span> << n << endl;
}
<span class="keywordflow">else</span> {
PIR R(m);
<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<PIR></a> M(R);
<a class="code" href="smith_8_c.html#ac924a785b0083e0f26562df322262428" title="Output matrix is determined by src which may be: "random-rough" This mat will have s...">Mat</a>(M, R, n, src, file, format);
<span class="keyword">typedef</span> list< PIR::Element > List;
List L;
<a name="_a9"></a><a class="code" href="class_lin_box_1_1_smith_form_local.html" title="Smith normal form (invariant factors) of a matrix over a local ring.">SmithFormLocal<PIR></a> SmithForm;
T.start();
SmithForm( L, M, R );
T.stop();
list<pair<PIR::Element, size_t> > p;
distinct(L.begin(), L.end(), p);
cout << <span class="stringliteral">"#"</span>;
display(p.begin(), p.end());
cout << <span class="stringliteral">"# local, PIR-Modular-int32_t("</span> << m << <span class="stringliteral">"), n = "</span> << n << endl;
}
cout << <span class="stringliteral">"T"</span> << n << <span class="stringliteral">"local"</span> << m << <span class="stringliteral">" := "</span>;
}
<span class="keywordflow">else</span> <span class="keywordflow">if</span> (algo == <span class="stringliteral">"2local"</span>) {
<a name="_a10"></a><a class="code" href="struct_lin_box_1_1_local2__32.html" title="Fast arithmetic mod 2^32, including gcd.">Local2_32</a> R;
<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<Local2_32></a> M(R);
<a class="code" href="smith_8_c.html#ac924a785b0083e0f26562df322262428" title="Output matrix is determined by src which may be: "random-rough" This mat will have s...">Mat</a>(M, R, n, src, file, format);
<span class="keyword">typedef</span> list< Local2_32::Element > List;
List L;
<a class="code" href="class_lin_box_1_1_smith_form_local.html" title="Smith normal form (invariant factors) of a matrix over a local ring.">SmithFormLocal<Local2_32></a> SmithForm;
T.start();
SmithForm( L, M, R );
T.stop();
list<pair<Local2_32::Element, size_t> > p;
distinct(L.begin(), L.end(), p);
cout << <span class="stringliteral">"#"</span>;
display(p.begin(), p.end());
cout << <span class="stringliteral">"# 2local, Local2_32, n = "</span> << n << endl;
cout << <span class="stringliteral">"T"</span> << n << <span class="stringliteral">"local2_32 := "</span>;
}
<span class="keywordflow">else</span> {
printf (<span class="stringliteral">"Unknown algorithms\n"</span>);
exit (-1);
}
T.print(cout); cout << <span class="stringliteral">";"</span> << endl;
<span class="keywordflow">return</span> 0 ;
}
<span class="keyword">template</span> <<span class="keyword">class</span> PIR>
<span class="keywordtype">void</span> <a class="code" href="smith_8_c.html#ac924a785b0083e0f26562df322262428" title="Output matrix is determined by src which may be: "random-rough" This mat will have s...">Mat</a>(<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<PIR></a>& M, PIR& R, <span class="keywordtype">int</span> n,
<span class="keywordtype">string</span> src, <span class="keywordtype">string</span> file, <span class="keywordtype">string</span> format) {
<span class="keyword">typename</span> PIR::Element zero; R.init(zero, 0);
<span class="keywordflow">if</span> (src == <span class="stringliteral">"random-rough"</span>) RandomRoughMat(M, R, n);
<span class="keywordflow">else</span> <span class="keywordflow">if</span> (src == <span class="stringliteral">"random"</span>) RandomFromDiagMat(M, R, n);
<span class="keywordflow">else</span> <span class="keywordflow">if</span> (src == <span class="stringliteral">"fib"</span>) RandomFibMat(M, R, n);
<span class="keywordflow">else</span> <span class="keywordflow">if</span> (src == <span class="stringliteral">"tref"</span>) TrefMat(M, R, n);
<span class="keywordflow">else</span> <span class="comment">// from file</span>
{
<span class="keywordtype">int</span> rdim, cdim;
std::ifstream in (file.c_str(), std::ios::in);
<span class="keywordflow">if</span> (! in) { cerr << <span class="stringliteral">"error: unable to open file"</span> << endl; exit(-1); }
in >> rdim >> cdim;
M. resize (rdim, cdim);
<a class="code" href="group__integers.html#gad62eceb96963b157a2357aba991f6d6e" title="Integers in LinBox.">integer</a> Val;
<span class="keywordflow">if</span> (format == <span class="stringliteral">"dense"</span> ) {
<span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i < rdim; ++ i)
<span class="keywordflow">for</span> ( <span class="keywordtype">int</span> j = 0; j < cdim; ++ j) {
in >> Val;
R. init (M[i][j], Val);
}
}
<span class="keywordflow">else</span> <span class="keywordflow">if</span> (format == <span class="stringliteral">"sparse"</span>) {
<span class="keywordtype">int</span> i, j;
<span class="keywordtype">char</span> mark;
in >> mark;
<a class="code" href="group__integers.html#gad62eceb96963b157a2357aba991f6d6e" title="Integers in LinBox.">LinBox::integer</a> val;
<span class="keywordflow">do</span> {
in >> i >> j;
in. ignore (1);
in >> val;
<span class="keywordflow">if</span> ( i == 0) <span class="keywordflow">break</span>;
R. init (M[i-1][j-1], val);
} <span class="keywordflow">while</span> (<span class="keyword">true</span>);
}
<span class="comment">//Krattenthaler's q^e matrices, given by exponent</span>
<span class="keywordflow">else</span> <span class="keywordflow">if</span> (format == <span class="stringliteral">"kdense"</span>) KratMat(M, R, n, in);
<span class="keywordflow">else</span> {
cout << <span class="stringliteral">"Format: "</span> << format << <span class="stringliteral">" Unknown\n"</span>;
exit (-1);
}
}
<span class="comment">/*show some entries</span>
<span class="comment"> for (int k = 0; k < 10; ++k)</span>
<span class="comment"> cout << M.getEntry(0,k) << " " << M.getEntry(M.rowdim()-1, M.coldim()-1 - k) << endl;</span>
<span class="comment"> cout << endl << M.rowdim() << " " << M.coldim() << endl;</span>
<span class="comment"> */</span>
<span class="comment">/* some row ops and some col ops */</span>
} <span class="comment">// Mat</span>
<span class="comment">// This mat will have s, near sqrt(n), distinct invariant factors,</span>
<span class="comment">// each repeated twice), involving the s primes 101, 103, ...</span>
<span class="keyword">template</span> <<span class="keyword">class</span> PIR>
<span class="keywordtype">void</span> RandomRoughMat(<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<PIR></a>& M, PIR& R, <span class="keywordtype">int</span> n) {
<span class="keyword">typename</span> PIR::Element zero; R.init(zero, 0);
M.<a name="a11"></a><a class="code" href="class_lin_box_1_1_dense_matrix_base.html#a0a4488dda5f370f56b9eab4fe8a54683" title="Resize the matrix to the given dimensions.">resize</a>(n, n, zero);
<span class="keywordflow">if</span> (n > 10000) {cerr << <span class="stringliteral">"n too big"</span> << endl; exit(-1);}
<span class="keywordtype">int</span> 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};
<span class="keywordflow">for</span> (<span class="keywordtype">int</span> j= 0, i = 0 ; i < n; ++j)
{
<span class="keyword">typename</span> PIR::Element v; R.init(v, jth_factor[25+j]);
<span class="keywordflow">for</span> (<span class="keywordtype">int</span> k = j ; k > 0 && i < n ; --k)
{ M[i][i] = v; ++i;
<span class="keywordflow">if</span> (i < n) {M[i][i] = v; ++i;}
}
}
scramble(M);
}
<span class="comment">// This mat will have the same nontrivial invariant factors as</span>
<span class="comment">// diag(1,2,3,5,8, ... 999, 0, 1, 2, ...).</span>
<span class="keyword">template</span> <<span class="keyword">class</span> PIR>
<span class="keywordtype">void</span> RandomFromDiagMat(<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<PIR></a>& M, PIR& R, <span class="keywordtype">int</span> n) {
<span class="keyword">typename</span> PIR::Element zero; R.init(zero, 0);
M.<a class="code" href="class_lin_box_1_1_dense_matrix_base.html#a0a4488dda5f370f56b9eab4fe8a54683" title="Resize the matrix to the given dimensions.">resize</a>(n, n, zero);
<span class="keywordflow">for</span> (<span class="keywordtype">int</span> i= 0 ; i < n; ++i)
R.init(M[i][i], i % 1000 + 1);
scramble(M);
}
<span class="comment">// This mat will have the same nontrivial invariant factors as</span>
<span class="comment">// diag(1,2,3,5,8, ... fib(k)), where k is about sqrt(n).</span>
<span class="comment">// The basic matrix is block diagonal with i-th block of order i and</span>
<span class="comment">// being a tridiagonal {-1,0,1} matrix whose snf = diag(i-1 1's, fib(i)),</span>
<span class="comment">// where fib(1) = 1, fib(2) = 2. But note that, depending on n,</span>
<span class="comment">// the last block may be truncated, thus repeating an earlier fibonacci number.</span>
<span class="keyword">template</span> <<span class="keyword">class</span> PIR>
<span class="keywordtype">void</span> RandomFibMat(<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<PIR></a>& M, PIR& R, <span class="keywordtype">int</span> n) {
<span class="keyword">typename</span> PIR::Element zero; R.init(zero, 0);
M.<a class="code" href="class_lin_box_1_1_dense_matrix_base.html#a0a4488dda5f370f56b9eab4fe8a54683" title="Resize the matrix to the given dimensions.">resize</a>(n, n, zero);
<span class="keyword">typename</span> PIR::Element one; R.init(one, 1);
<span class="keywordflow">for</span> (<span class="keywordtype">int</span> i= 0 ; i < n; ++i) M[i][i] = one;
<span class="keywordtype">int</span> j = 1, k = 0;
<span class="keywordflow">for</span> (<span class="keywordtype">int</span> i= 0 ; i < n-1; ++i) {
<span class="keywordflow">if</span> ( i == k) {
M[i][i+1] = zero;
k += ++j;
}
<span class="keywordflow">else</span> {
M[i][i+1] = one;
R.negin(one);
}
R.neg(M[i+1][i], M[i][i+1]);
}
scramble(M);
}
<span class="keyword">template</span> < <span class="keyword">class</span> Ring >
<span class="keywordtype">void</span> scramble(<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<Ring></a>& M)
{
Ring R = M.<a name="a12"></a><a class="code" href="class_lin_box_1_1_dense_matrix.html#a66be394bace6f92c18729a119581d7bd" title="Retrieve the field over which this matrix is defined.">field</a>();
<span class="keywordtype">int</span> N,n = (int)M.<a name="a13"></a><a class="code" href="class_lin_box_1_1_dense_matrix.html#a0fee90e6b8ef7fcd8d76a7f3e1f268e6" title="Get the number of rows in the matrix.">rowdim</a>(); <span class="comment">// number of random basic row and col ops.</span>
N = n;
<span class="keywordflow">for</span> (<span class="keywordtype">int</span> k = 0; k < N; ++k) {
<span class="keywordtype">int</span> i = rand()%(int)M.<a class="code" href="class_lin_box_1_1_dense_matrix.html#a0fee90e6b8ef7fcd8d76a7f3e1f268e6" title="Get the number of rows in the matrix.">rowdim</a>();
<span class="keywordtype">int</span> j = rand()%(int)M.<a name="a14"></a><a class="code" href="class_lin_box_1_1_dense_matrix.html#a32edb490d3597f5553152d14b102e227" title="Get the number of columns in the matrix.">coldim</a>();
<span class="keywordflow">if</span> (i == j) <span class="keywordflow">continue</span>;
<span class="comment">// M*i += alpha M*j and Mi* += beta Mj</span>
<span class="comment">//int a = rand()%2;</span>
<span class="keywordtype">int</span> a = 0;
<span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> l = 0; l < M.<a class="code" href="class_lin_box_1_1_dense_matrix.html#a0fee90e6b8ef7fcd8d76a7f3e1f268e6" title="Get the number of rows in the matrix.">rowdim</a>(); ++l) {
<span class="keywordflow">if</span> (a)
R.subin(M[l][i], M[l][j]);
<span class="keywordflow">else</span>
R.addin(M[l][i], M[l][j]);
<span class="comment">//K.axpy(c, M.getEntry(l, i), x, M.getEntry(l, j));</span>
<span class="comment">//M.setEntry(l, i, c);</span>
}
<span class="comment">//a = rand()%2;</span>
<span class="keywordflow">for</span> (<span class="keywordtype">size_t</span> l = 0; l < M.<a class="code" href="class_lin_box_1_1_dense_matrix.html#a32edb490d3597f5553152d14b102e227" title="Get the number of columns in the matrix.">coldim</a>(); ++l) {
<span class="keywordflow">if</span> (a)
R.subin(M[i][l], M[j][l]);
<span class="keywordflow">else</span>
R.addin(M[i][l], M[j][l]);
}
}
std::ofstream out(<span class="stringliteral">"matrix"</span>, std::ios::out);
out << n << <span class="stringliteral">" "</span> << n << <span class="stringliteral">"\n"</span>;
<span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i < n; ++ i) {
<span class="keywordflow">for</span> ( <span class="keywordtype">int</span> j = 0; j < n; ++ j) {
R. write(out, M[i][j]);
out << <span class="stringliteral">" "</span>;
}
out << <span class="stringliteral">"\n"</span>;
}
<span class="comment">//}</span>
}
<span class="comment">// special mats tref and krat</span>
<span class="comment">// Trefethen's challenge #7 mat (primes on diag, 1's on 2^e bands).</span>
<span class="keyword">template</span> <<span class="keyword">class</span> PIR>
<span class="keywordtype">void</span> TrefMat(<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<PIR></a>& M, PIR& R, <span class="keywordtype">int</span> n) {
<span class="keyword">typename</span> PIR::Element zero; R.init(zero, 0);
M.<a class="code" href="class_lin_box_1_1_dense_matrix_base.html#a0a4488dda5f370f56b9eab4fe8a54683" title="Resize the matrix to the given dimensions.">resize</a>(n, n, zero);
std::vector<int> power2;
<span class="keywordtype">int</span> i = 1;
<span class="keywordflow">do</span> {
power2. push_back(i);
i *= 2;
} <span class="keywordflow">while</span> (i < n);
std::ifstream in (<span class="stringliteral">"prime"</span>, std::ios::in);
<span class="keywordflow">for</span> ( i = 0; i < n; ++ i)
in >> M[i][i];
std::vector<int>::iterator p;
<span class="keywordflow">for</span> ( i = 0; i < n; ++ i) {
<span class="keywordflow">for</span> ( p = power2. begin(); (p != power2. end()) && (*p <= i); ++ p)
M[i][i - *p] = 1;
<span class="keywordflow">for</span> ( p = power2. begin(); (p != power2. end()) && (*p < n - i); ++ p)
M[i][i + *p] = 1;
}
}
<span class="keyword">struct </span>pwrlist
{
vector<integer> m;
pwrlist(<a class="code" href="group__integers.html#gad62eceb96963b157a2357aba991f6d6e" title="Integers in LinBox.">integer</a> q)
{ m.push_back(1); m.push_back(q); <span class="comment">//cout << "pwrlist " << m[0] << " " << m[1] << endl;</span>
}
<a class="code" href="group__integers.html#gad62eceb96963b157a2357aba991f6d6e" title="Integers in LinBox.">integer</a> operator[](<span class="keywordtype">int</span> e)
{
<span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = (<span class="keywordtype">int</span>)m.size(); i <= e; ++i) m.push_back(m[1]*m[i-1]);
<span class="keywordflow">return</span> m[e];
}
};
<span class="comment">// Read "1" or "q" or "q^e", for some (small) exponent e.</span>
<span class="comment">// Return value of the power of q at q = _q.</span>
<span class="keyword">template</span> <<span class="keyword">class</span> num>
num& qread(num& Val, pwrlist& M, istream& in)
{
<span class="keywordtype">char</span> c;
in >> c; <span class="comment">// next nonwhitespace</span>
<span class="keywordflow">if</span> (c == <span class="charliteral">'0'</span>) <span class="keywordflow">return</span> Val = 0;
<span class="keywordflow">if</span> (c == <span class="charliteral">'1'</span>) <span class="keywordflow">return</span> Val = 1;
<span class="keywordflow">if</span> (c != <span class="charliteral">'p'</span> && c != <span class="charliteral">'q'</span>) { cout << <span class="stringliteral">"exiting due to unknown char "</span> << c << endl; exit(-1);}
in.get(c);
<span class="keywordflow">if</span> (c !=<span class="charliteral">'^'</span>) {in.putback(c); <span class="keywordflow">return</span> Val = M[1];}
<span class="keywordflow">else</span>
{ <span class="keywordtype">int</span> expt; in >> expt;
<span class="keywordflow">return</span> Val = M[expt];
};
}
<span class="keyword">template</span> <<span class="keyword">class</span> PIR>
<span class="keywordtype">void</span> KratMat(<a class="code" href="class_lin_box_1_1_dense_matrix.html" title="Blackbox interface to dense matrix representation.">DenseMatrix<PIR></a>& M, PIR& R, <span class="keywordtype">int</span> q, istream& in)
{
pwrlist pwrs(q);
<span class="keywordflow">for</span> (<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < M.<a name="a15"></a><a class="code" href="class_lin_box_1_1_dense_matrix.html#a0fee90e6b8ef7fcd8d76a7f3e1f268e6" title="Get the number of rows in the matrix.">rowdim</a>(); ++ i)
<span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < M.<a name="a16"></a><a class="code" href="class_lin_box_1_1_dense_matrix.html#a32edb490d3597f5553152d14b102e227" title="Get the number of columns in the matrix.">coldim</a>(); ++ j) {
<span class="keywordtype">int</span> Val;
qread(Val, pwrs, in);
R. init (M[i][j], Val);
}
}
<span class="keyword">template</span><<span class="keyword">class</span> I1, <span class="keyword">class</span> Lp>
<span class="keywordtype">void</span> distinct (I1 a, I1 b, Lp& c)
{ <span class="keyword">typename</span> I1::value_type e;
<span class="keywordtype">size_t</span> count = 0;
<span class="keywordflow">if</span> (a != b) {e = *a; ++a; count = 1;}
<span class="keywordflow">else</span> <span class="keywordflow">return</span>;
<span class="keywordflow">while</span> (a != b)
{ <span class="keywordflow">if</span> (*a == e) ++count;
<span class="keywordflow">else</span>
{ c.push_back(<span class="keyword">typename</span> Lp::value_type(e, count));
e = *a; count = 1;
}
++a;
}
c.push_back(<span class="keyword">typename</span> Lp::value_type(e, count));
<span class="keywordflow">return</span>;
}
<span class="keyword">template</span> <<span class="keyword">class</span> I>
<span class="keywordtype">void</span> display(I b, I e)
{ cout << <span class="stringliteral">"("</span>;
<span class="keywordflow">for</span> (I p = b; p != e; ++p) cout << p->first << <span class="stringliteral">" "</span> << p->second << <span class="stringliteral">", "</span>;
cout << <span class="stringliteral">")"</span> << endl;
}
</pre></div> </div>
</div>
<hr class="footer"/><address class="footer"><small>Generated on Tue Aug 30 2011 for linbox by 
<a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.7.4 </small></address>
</body>
</html>
