1/* linbox/algorithms/localsmith.h
2 * Copyright (C) LinBox
3 *
4 * Written by David Saunders
5 *
6 * ------------------------------------
7 *
8 *
9 * ========LICENCE========
10 * This file is part of the library LinBox.
11 *
12 * LinBox is free software: you can redistribute it and/or modify
13 * it under the terms of the  GNU Lesser General Public
14 * License as published by the Free Software Foundation; either
15 * version 2.1 of the License, or (at your option) any later version.
16 *
17 * This library is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
20 * Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with this library; if not, write to the Free Software
24 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
25 * ========LICENCE========
26 *.
27 */
28
29#ifndef __LINBOX_smith_form_local2_H
30#define __LINBOX_smith_form_local2_H
31
32
33#include <vector>
34#include <list>
35//#include <algorithm>
36
37#include "linbox/ring/local2_32.h"
38
39namespace LinBox
40{
41
42	/**
43	  \brief Smith normal form (invariant factors) of a matrix over a local ring.
44	  */
45	template<class LocalRing>
46	class SmithFormLocal;
47
48	template <>
49	class SmithFormLocal<Local2_32>
50	{
51	public:
52		typedef Local2_32 LocalPIR;
53		typedef LocalPIR::Element Elt;
54
55		template<class Matrix>
56		std::list<Elt>& operator()(std::list<Elt>& L, Matrix& A, const LocalPIR& R)
57		{
58			Elt d; R.assign(d, R.one);
59			Elt *p = &(A[0][0]);
60			return smithStep(L, d, p, A.rowdim(), A.coldim(), A.getStride(), R);
61		}
62
63		std::list<Elt>&
64		smithStep(std::list<Elt>& L, Elt& d, Elt* Ap, size_t m, size_t n, size_t stride, const LocalPIR& R)
65		{
66			if ( m == 0 || n == 0 )
67				return L;
68
69			LocalPIR::Exponent g = LocalPIR::Exponent(32); //R.init(g); // must change to 2^31 maybe.
70			size_t i, j;
71			/* Arguably this search order should be reversed to increase the likelyhood of no col swap,
72			   assuming row swaps cheaper.  Not so, however on my example. -bds 11Nov */
73			for ( i = 0; i != m; ++i)
74			{
75				for (j = 0; j != n; ++j)
76				{
77					R.gcdin(g, Ap[i*stride + j]);
78					if ( R.isUnit(g) ) break;
79				}
80				if ( R.isUnit(g) ) break;
81			}
82			if ( R.isZero(g) )
83			{
84				L.insert(L.end(), (m < n) ? m : n, 0);
85				return L;
86			}
87			if ( i != m ) // g is a unit and, because this is a local ring,
88			// value at which this first happened also is a unit.
89			{ // put pivot in 0,0 position
90				if ( i != 0 ) // swap rows
91					std::swap_ranges(Ap, Ap+n, Ap + i*stride);
92				if ( j != 0 ) // swap cols
93					for(size_t k = 0; k != m; ++k)
94						std::swap(Ap[k*stride + 0], Ap[k*stride + j]);
95
96				// elimination step - crude and for dense only - fix later
97				// Want to use a block method or "left looking" elimination.
98				Elt f; R.inv(f, Ap[0*stride + 0] );
99				R.negin(f);
100
101				// normalize first row to -1, ...
102				for ( j = 0; j != n; ++j)
103					R.mulin(Ap[0*stride + j], f);
104
105				// eliminate in subsequent rows
106				for ( i = 1; i != m; ++i)
107				{
108					f = Ap[i*stride + 0];
109					for ( j = 0; j != n; ++j)
110						R.axpyin( Ap[i*stride +j], f, Ap[0*stride +j] );
111				}
112				L.push_back(d);
113				return smithStep(L, d, Ap + stride+1,m-1, n-1, stride, R);
114			}
115			else
116			{
117				for ( i = 0; i != m; ++i)
118					for ( j = 0; j != n; ++j)
119					{
120						R.divin(Ap[i*stride + j], g);
121					}
122				return smithStep(L, R.mulin(d, g), Ap, m, n, stride, R);
123			}
124		}
125
126	}; // end SmithFormLocal
127
128} // end LinBox
129
130#endif // __LINBOX_smith_form_local2_H
131
132// Local Variables:
133// mode: C++
134// tab-width: 4
135// indent-tabs-mode: nil
136// c-basic-offset: 4
137// End:
138// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
139