1 /* Copyright (C) LinBox
2 *
3 * Author: bds
4 *
5 * ========LICENCE========
6 * This file is part of the library LinBox.
7 *
8 * LinBox is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU Lesser General Public
10 * License as published by the Free Software Foundation; either
11 * version 2.1 of the License, or (at your option) any later version.
12 *
13 * This library is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with this library; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 * ========LICENCE========
22 */
23
24 /*! @file tests/test-smith-form.h
25 * @ingroup tests
26 * @brief tools for making matrix with known SNF.
27 */
28
29 #ifndef __TEST_SMITH_FORM_H
30 #define __TEST_SMITH_FORM_H
31 #include <linbox/linbox-config.h>
32
33 #include "linbox/util/commentator.h"
34 #include "linbox/matrix/dense-matrix.h"
35 #include "linbox/vector/blas-vector.h"
36 #include "linbox/solutions/smith-form.h"
37
38 using std::endl;
39 using namespace LinBox;
40
41 template <class PIR> // This is for PIR = Z or Z_n
makeBumps(BlasVector<PIR> & b,int choice)42 BlasVector<PIR> & makeBumps(BlasVector<PIR> & b, int choice) {
43 const PIR & R = b.field();
44 typename PIR::Element two, three, nine, x;
45 R.init(two,2);
46 R.init(three,3);
47 R.init(nine,9);
48 R.init(x,202);
49 // b is a single row
50 size_t n = b.size();
51 switch (choice) {
52 case 0: // all zero
53 for(size_t i = 0; i < n; ++i) b.setEntry(i,R.zero);
54 break;
55 case 1: // identity
56 for(size_t i = 0; i < n; ++i) b.setEntry(i,R.one);
57 break;
58 case 2: // powers of 2
59 for(size_t i = 0; i < n; ++i) b.setEntry(i,two);
60 break;
61 case 3: // random followed by 202,0. Random part is largely 1's.
62 for(size_t i = 0; i < n-2; ++i) {
63 size_t r = rand()%20;
64 if (r < 17) b.setEntry(i,R.one);
65 if (r == 17) b.setEntry(i,two);
66 if (r == 18) b.setEntry(i,three);
67 if (r == 19) b.setEntry(i,nine);
68 }
69 b.setEntry(n-2,x);
70 b.setEntry(n-1,R.zero);
71 }
72 // negate a few
73 for (size_t k = rand()%4; k > 0; --k){
74 size_t i = rand()%n;
75 b.getEntry(x,i);
76 b.setEntry(i,R.negin(x));
77 }
78
79 return b;
80 }
81
82
83 // For any PIR, build an increasing sequence of smith invariants d from "bumps" b.
84 template <class PIR>
prefixProduct(BlasVector<PIR> & d,const BlasVector<PIR> & b)85 BlasVector<PIR> & prefixProduct (BlasVector<PIR> & d, const BlasVector<PIR> & b) {
86 const PIR& R = d.field();
87 typename PIR::Element x,y; R.init(x); R.init(y);
88 d.setEntry(0,b.getEntry(x,0));
89 for (size_t i = 1; i < d.size(); ++i){
90 d.getEntry(x,i-1);
91 b.getEntry(y,i);
92 d.setEntry(i,R.mulin(x, y));
93 }
94 return d;
95 }
96
97 // Generate A with snf = diag(d) (up to sign), based on the bumps.
98 // Think of bumps[i] as s_i/s_{i-1}, quotient of smith invariants.
99 // The lumps are used for off diagonal entries in L,U (triangular scramblers).
100 template <class PIR>
makeSNFExample(DenseMatrix<PIR> & A,BlasVector<PIR> & d,const BlasVector<PIR> & bumps,const BlasVector<PIR> & lumps)101 void makeSNFExample(DenseMatrix<PIR>& A,
102 BlasVector<PIR> & d,
103 const BlasVector<PIR> & bumps,
104 const BlasVector<PIR> & lumps) {
105 //LinBox::VectorWrapper::ensureDim (d, bumps.size());
106 //LinBox::VectorWrapper::ensureDim (d, std::min(A.rowdim(), A.coldim()));
107 prefixProduct(d, bumps);
108
109 // make A = UDL for random unimodular L,U
110 const PIR & R = A.field();
111 DenseMatrix<PIR> L(R, A.coldim(), A.coldim()),
112 U(R, A.rowdim(), A.rowdim());
113 typename PIR::Element x; R.init(x);
114 size_t i, j, k;
115 k = lumps.size();
116 A.zero();
117 for(i = 0; i < d.size(); ++i) A.setEntry(i,i,d.getEntry(x,i));
118
119
120 L.zero();
121 for(i = 0; i < L.rowdim(); ++i) L.setEntry(i,i,R.one);
122 for (i = 0; i < L.rowdim(); ++ i)
123 for (j = 0; j < i; ++ j) L.setEntry(i,j,lumps[rand()%k]);
124
125 U.zero();
126 for(i = 0; i < U.rowdim(); ++i) U.setEntry(i,i,R.one);
127 for (i = 0; i < U.rowdim(); ++ i)
128 for (j = i+1; j < U.coldim(); ++ j) U.setEntry(i,j,lumps[rand()%k]);
129
130
131
132 // A <- UAL
133 BlasMatrixDomain<PIR> MD(R);
134 MD.mulin_left(A,L);
135 MD.mulin_right(U,A);
136
137 for (i = 0; i < d.size(); ++ i)
138 d.setEntry(i,R.abs(x, d.getEntry(x,i)));
139 // Now A is matrix equivalent to diag prefix product of bumps.
140 // Now d is SNF diagonal (vector of invariants) for A.
141 }
142
143 template <class PIR>
checkSNFExample(const BlasVector<PIR> & d,const BlasVector<PIR> & x)144 bool checkSNFExample( const BlasVector<PIR>& d, const BlasVector<PIR>& x ){
145
146 std::ostream & report =
147 #ifndef DISABLE_COMMENTATOR
148 commentator().report()
149 #else
150 std::cerr
151 #endif
152 ;
153
154 VectorDomain<PIR> VD(d.field());
155
156 report << "Expected smith form: ";
157 VD.write (report, d) << endl;
158
159 report << "Computed Smith form: ";
160 VD. write (report, x) << endl;
161
162 if (not VD.areEqual (d, x)) {
163 report << "ERROR: Computed not as Expected." << endl;
164 return false;
165 }
166
167 report << "PASSED." << endl;
168
169 return true;
170 }
171
172 template <class PIR>
checkSNFExample(const LinBox::SmithList<PIR> & d,const LinBox::SmithList<PIR> & x,const PIR & R)173 bool checkSNFExample( const LinBox::SmithList<PIR>& d, const LinBox::SmithList<PIR>& x, const PIR& R){
174
175 std::ostream & report =
176 #ifndef DISABLE_COMMENTATOR
177 commentator().report()
178 #else
179 std::clog
180 #endif
181 ;
182 report << "Expected smith form SL: " << '{';
183 for(auto const & sit: d) report << '{' << sit.first << ',' << sit.second << '}';
184 report << '}' << std::endl;
185
186 report << "Computed Smith form SL: " << '{';
187 for(auto const & sit: x) report << '{' << sit.first << ',' << sit.second << '}';
188 report << '}' << std::endl;
189
190
191 auto dit=d.begin();
192 auto xit=x.begin();
193 bool pass=true;
194 for( ; dit != d.end(); ++dit, ++xit) {
195 if (!R.areEqual(dit->first,xit->first)
196 || dit->second!=xit->second) {
197 R.write( R.write(
198 report << "ERROR: Computed not as Expected. "
199 << '(', dit->first) << ',' << dit->second << ')'
200 << " != "
201 << '(', xit->first) << ',' << xit->second << ')'
202 << endl;
203 pass = false;
204 }
205 }
206
207 report << "PASSED." << endl;
208
209 return pass;
210 }
211 #endif // __TEST_SMITH_FORM_H
212
213 // Local Variables:
214 // mode: C++
215 // tab-width: 4
216 // indent-tabs-mode: nil
217 // c-basic-offset: 4
218 // End:
219 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
220