1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
|
////////////////////////////////////////////////////////////////////////////////
//
// StrictStairCaseMatrixTrans.cc
//
// produced: 30/01/2020 jr
//
////////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <ctype.h>
#include <string.h>
#include "Vector.hh"
#include "StrictStairCaseMatrixTrans.hh"
namespace topcom {
// we eliminate the columns with indices cidx to cidx + cno (exclusive) as follows
// with transformation information:
// ------------------------------------------------------------------------------
// | already eliminated part | uneliminated part
// | 0 1 .. ridx cidx-1 | cidx ... cidx+cno
// ----------------------------------------------------------------------------------
// | * 0 ... 0 ......... 0 | 0 0 ... 0
// | * * 0 . . ......... . | . . .
// | * * * 0 . ......... . | . . .
// | * * * * 0 ......... . | 0 0 0
// rdix | * * * *(*)0 ....... . |(*)(*) (*)
// | * * * * * * 0 ..... . | * * *
// | * * * * * * * 0 ... 0 | * * ... *
// ------------------------------------------------------------------------------
// 0 | the columns in this part contain the coefficients
// | of the representation of the column above as
// | a linear combination of original columns;
// transf | used mainly for the computation of circuit signatures:
// | below the left-most zero column above
// | we find the dependency in this transformation column,
// | and its signs are a circuit signature;
// cidx-1 | initially allocated with rowdim+1 rows, enlarged if needed
void StrictStairCaseMatrixTrans::_eliminate(const parameter_type ridx,
const parameter_type cidx,
const parameter_type cno) {
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "called elimination to eliminate the elements (" << ridx << ", " << cidx << ")"
<< " through (" << ridx << ", " << cidx + cno - 1 << ")"
<< " ..." << std::endl;
MessageStreams::debug() << "starting elimination with start pivot (" << ridx << ", " << ridx << ") out of "
<< rowdim() << ":" << std::endl;
MessageStreams::debug() << message::unlock;
#endif
const parameter_type rdim = rowdim();
const parameter_type cdim = coldim();
parameter_type rowidx = ridx;
while (rowidx < rdim) {
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "current pivot (" << rowidx << ", " << ridx << ")" << std::endl;
MessageStreams::debug() << message::unlock;
#endif
if (is_zero((*this)(rowidx,ridx))) {
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "eraser = " << (*this)(rowidx,ridx) << " == FieldConstants::ZERO -> try colperm; row : " << rowidx << ", col: " << ridx << std::endl;
MessageStreams::debug() << message::unlock;
#endif
for (parameter_type k = cidx; k < cidx + cno; ++k) {
if (!is_zero((*this)(rowidx,k))) {
Matrix::swap_cols(ridx,k);
_coefficient *= FieldConstants::MINUSONE;
_transformation.swap_cols(ridx,k);
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "eraser = " << (*this)(rowidx,ridx) << " != FieldConstants::ZERO -> colperm successful; row : " << rowidx << ", col: " << ridx << std::endl;
pretty_print(MessageStreams::debug());
MessageStreams::debug() << message::unlock;
#endif
return;
}
}
if (is_zero((*this)(rowidx,ridx))) {
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "eraser = " << (*this)(rowidx,ridx) << " == FieldConstants::ZERO -> colperm unsuccessful; row : " << rowidx << ", col: " << ridx << std::endl;
MessageStreams::debug() << message::unlock;
#endif
}
else {
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "eraser = " << (*this)(rowidx,ridx) << " != FieldConstants::ZERO; row : " << rowidx << ", col: " << ridx << std::endl;
MessageStreams::debug() << message::unlock;
#endif
break;
}
}
else {
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "eraser = " << (*this)(rowidx,ridx) << " != FieldConstants::ZERO; row : " << rowidx << ", col: " << ridx << std::endl;
MessageStreams::debug() << message::unlock;
#endif
break;
}
++rowidx;
}
if (rowidx == rowdim()) {
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "all potential erasers == FieldConstants::ZERO for all possible swaps of columns: zero columns; row : all, col: " << ridx << std::endl;
MessageStreams::debug() << message::unlock;
#endif
return;
}
const Field& eraser = (*this)(rowidx,ridx); // do not copy; careful: eraser changes with that matrix element!
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "eraser = " << eraser << "; row : " << rowidx << ", col: " << ridx << std::endl;
MessageStreams::debug() << message::unlock;
#endif
for (parameter_type j = cidx; j < cidx + cno; ++j) {
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "elimination of column " << j << " (" << j - cidx + 1 << " of " << cno << ") many columns to eliminate:" << std::endl;
MessageStreams::debug() << message::unlock;
#endif
const Field& delinquent = (*this)(rowidx,j); // do not copy; careful: delinquent changes with that matrix element!
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "delinquent = " << delinquent << "; row : " << rowidx << ", col: " << j << std::endl;
MessageStreams::debug() << message::unlock;
#endif
if (is_zero(delinquent)) {
continue;
}
const Field multiplier = delinquent / eraser;
for (parameter_type k = 0; k < cdim; ++k) {
_transformation(k,j) -= multiplier * _transformation(k,ridx);
}
for (parameter_type k = rowidx + 1; k < rdim; ++k) {
(*this)(k,j) -= multiplier * (*this)(k,ridx);
}
(*this)(rowidx,j) = FieldConstants::ZERO;
}
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "after step " << ridx << " of stair case transformation: " << std::endl;
(*this).pretty_print(MessageStreams::debug());
MessageStreams::debug() << message::unlock;
#endif
}
StrictStairCaseMatrixTrans& StrictStairCaseMatrixTrans::augment(const Vector& vector) {
#ifdef INDEX_CHECK
assert(rowdim() == vector.dim());
#endif
if (empty()) {
// if the matrix has no columns yet,
// we can assign the StrictStairCaseMatrix
// constructed from a single vector:
return operator=(StrictStairCaseMatrixTrans(vector));
}
const parameter_type cdim(coldim());
const parameter_type rdim(rowdim());
const parameter_type m = cdim < rdim ? cdim : rdim;
Matrix::augment(vector);
_transformation.augment(Vector(vector.dim() + 1, FieldConstants::ZERO));
if (_transformation.rowdim() < cdim + 1) {
for (parameter_type col = 0; col < _transformation.coldim(); ++col) {
_transformation.col(col).resize(cdim + 1);
}
}
_transformation(cdim, cdim) = FieldConstants::ONE;
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "before strict stair case transformation:" << std::endl;
pretty_print(MessageStreams::debug());
MessageStreams::debug() << message::unlock;
#endif
for (parameter_type i = 0; i < m; ++i) {
this->_eliminate(i, cdim, 1);
}
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "after strict stair case transformation: " << std::endl;
(*this).pretty_print(MessageStreams::debug());
MessageStreams::debug() << message::unlock;
#endif
return *this;
}
StrictStairCaseMatrixTrans& StrictStairCaseMatrixTrans::augment(const Matrix& matrix) {
#ifdef INDEX_CHECK
assert(rowdim() == matrix.rowdim());
#endif
for (parameter_type col = 0; col < matrix.coldim(); ++col) {
augment(matrix.col(col));
}
return *this;
}
std::ostream& StrictStairCaseMatrixTrans::pretty_print(std::ostream& ost) const {
ost << "matrix: " << '\n';
Matrix::pretty_print(ost);
ost << "coefficient: " << _coefficient << '\n';
ost << "transformation:" << '\n';
_transformation.pretty_print(ost);
ost << '\n';
return ost;
}
std::istream& StrictStairCaseMatrixTrans::read(std::istream& ist) {
char dash, arrow;
if (!StairCaseMatrix::read(ist)) {
MessageStreams::forced() << message::lock;
MessageStreams::forced() << "StrictStairCaseMatrixTrans::read(std::istream&): "
<< "error while reading StairCaseMatrix - exiting" << std::endl;
MessageStreams::forced() << message::unlock;
exit(1);
}
if (!(ist >> std::ws >> dash >> arrow >> std::ws)) {
MessageStreams::forced() << message::lock;
MessageStreams::forced() << "StrictStairCaseMatrixTrans::read(std::istream&): "
<< "expected `" << ContainerChars::map_chars
<< "' - exiting" << std::endl;
MessageStreams::forced() << message::unlock;
exit(1);
}
if (!(ist >> _transformation)) {
MessageStreams::forced() << message::lock;
MessageStreams::forced() << "StrictStairCaseMatrixTrans::read(std::istream&): "
<< "error while reading _transformation - exiting" << std::endl;
MessageStreams::forced() << message::unlock;
exit(1);
}
return ist;
}
std::ostream& StrictStairCaseMatrixTrans::write(std::ostream& ost) const {
StairCaseMatrix::write(ost);
ost << ContainerChars::map_chars << _transformation;
return ost;
}
}; // namespace topcom
// eof StrictStairCaseMatrixTrans.cc
|
