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
|
////////////////////////////////////////////////////////////////////////////////
//
// StrictStairCaseMatrix.cc
//
// produced: 30/01/2020 jr
//
////////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <ctype.h>
#include <string.h>
#include "StrictStairCaseMatrix.hh"
namespace topcom {
// we eliminate the columns with indices cidx to cidx + cno (exclusive) as follows:
// ------------------------------------------------------------------------------
// already eliminated part | uneliminated part
// 0 1 .. ridx cidx-1 | cidx ... cidx+cno-1
// ------------------------------------------------------------------------------
// * 0 ... 0 ......... 0 | 0 0 ... 0
// * * 0 . . ......... . | . . .
// * * * 0 . ......... . | . . .
// * * * * 0 ......... . | 0 0 0
// * * * *(*)0 ....... . |(*)(*) (*)
// * * * * * * 0 ..... . | * * *
// * * * * * * * 0 ... . | * * ... *
void StrictStairCaseMatrix::_eliminate(const parameter_type ridx,
const parameter_type cidx,
const parameter_type cno) {
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "starting elimination of columns "
<< "[" << cidx << ", " << cidx + cno << ")"
<< " with start pivot (" << ridx << ", " << ridx << ") out of "
<< rowdim() << ":" << std::endl;
MessageStreams::debug() << message::unlock;
#endif
const parameter_type n = rowdim();
parameter_type rowidx = ridx;
while (rowidx < n) {
#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;
#ifdef SUPER_VERBOSE
MessageStreams::debug() << message::lock;
MessageStreams::debug() << "eraser = " << (*this)(rowidx,ridx) << " != FieldConstants::ZERO -> colperm successful; row : " << rowidx << ", col: " << ridx << std::endl;
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) {
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 = rowidx + 1; k < n; ++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
}
StrictStairCaseMatrix& StrictStairCaseMatrix::augment(const Matrix& matrix) {
if (matrix.coldim() == 0) {
return *this;
}
#ifdef INDEX_CHECK
assert(rowdim() == matrix.rowdim());
assert(coldim() + matrix.coldim() <= rowdim() + 1);
#endif
const parameter_type first_uneliminated_col(coldim());
// eliminate only up to lower triangular form:
parameter_type max_elim_rows = first_uneliminated_col;
if (max_elim_rows > rowdim()) {
// more than rowdim() many rows cannot be eliminated:
max_elim_rows = rowdim();
}
Matrix::augment(matrix);
#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 < max_elim_rows; ++i) {
this->_eliminate(i, first_uneliminated_col, matrix.coldim());
}
#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;
}
// determinant form induced by the first rowdim - 1 columns:
const Field StrictStairCaseMatrix::valuation(const Vector& vector) const {
StrictStairCaseMatrix tmp(*this, 0, rowdim() - 1);
tmp.augment(vector);
return tmp.left_upper_det();
}
}; // namespace topcom
// eof StrictStairCaseMatrix.cc
|
