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
|
////////////////////////////////////////////////////////////////////////////////
//
// StairCaseMatrix.cc
//
// produced: 13/03/98 jr
// last change: 13/03/98 jr
//
////////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <ctype.h>
#include <string.h>
#include "StairCaseMatrix.hh"
void StairCaseMatrix::_eliminate(const size_type& ridx, const size_type& cidx, const size_type& cno) {
const size_type n = rowdim();
if (is_zero((*this)(ridx,ridx))) {
#ifdef SUPER_VERBOSE
std::cerr << "eraser = " << (*this)(ridx,ridx) << " == ZERO -> try colperm; row : " << ridx << ", col: " << ridx << std::endl;
#endif
for (size_type k = cidx; k < cidx + cno; ++k) {
if (!is_zero((*this)(ridx,k))) {
(*this).swap_cols(ridx,k);
_coefficient *= MINUSONE;
#ifdef SUPER_VERBOSE
std::cerr << "eraser = " << (*this)(ridx,ridx) << " == ZERO -> colperm successful; row : " << ridx << ", col: " << ridx << std::endl;
#endif
return;
}
}
if (is_zero((*this)(ridx,ridx))) {
#ifdef SUPER_VERBOSE
std::cerr << "eraser = " << (*this)(ridx,ridx) << " == ZERO -> colperm unsuccessful; row : " << ridx << ", col: " << ridx << std::endl;
#endif
return;
}
#ifdef SUPER_VERBOSE
else {
std::cerr << "eraser = " << (*this)(ridx,ridx) << " != ZERO; row : " << ridx << ", col: " << ridx << std::endl;
}
#endif
}
#ifdef SUPER_VERBOSE
else {
std::cerr << "eraser = " << (*this)(ridx,ridx) << " != ZERO; row : " << ridx << ", col: " << ridx << std::endl;
}
#endif
const Field& eraser = (*this)(ridx,ridx);
#ifdef SUPER_VERBOSE
std::cerr << "eraser = " << eraser << "; row : " << ridx << ", col: " << ridx << std::endl;
#endif
for (size_type j = cidx; j < cidx + cno; ++j) {
const Field& delinquent = (*this)(ridx,j);
#ifdef SUPER_VERBOSE
std::cerr << "delinquent = " << delinquent << "; row : " << ridx << ", col: " << j << std::endl;
#endif
if (is_zero(delinquent)) {
continue;
}
for (size_type k = ridx + 1; k < n; ++k) {
(*this)(k,j) -= (*this)(k,ridx) * delinquent / eraser;
}
(*this)(ridx,j) = ZERO;
}
#ifdef SUPER_VERBOSE
std::cerr << "after step " << ridx << " of stair case transformation: " << std::endl;
(*this).pretty_print(std::cerr);
#endif
}
StairCaseMatrix& StairCaseMatrix::augment(const Matrix& matrix) {
if (matrix.coldim() == 0) {
return *this;
}
#ifdef INDEX_CHECK
assert(rowdim() == matrix.rowdim());
assert(coldim() + matrix.coldim() <= rowdim());
#endif
const size_type m = coldim();
Matrix::augment(matrix);
#ifdef SUPER_VERBOSE
std::cerr << "before stair case transformation:" << std::endl;
pretty_print(std::cerr);
#endif
for (size_type i = 0; i < m; ++i) {
_eliminate(i, m, matrix.coldim());
}
#ifdef SUPER_VERBOSE
std::cerr << "after stair case transformation: " << std::endl;
(*this).pretty_print(std::cerr);
#endif
return *this;
}
const Field det(const StairCaseMatrix& matrix) {
#ifdef INDEX_CHECK
assert(matrix.coldim() == matrix.rowdim());
#endif
Field result(ONE);
for (size_type i = 0; i < matrix.coldim(); ++i) {
result *= matrix(i,i);
if (is_zero(result)) {
return result;
}
}
return result * matrix._coefficient;
}
std::ostream& StairCaseMatrix::pretty_print(std::ostream& ost) {
ost << "matrix: " << '\n';
Matrix::pretty_print(ost);
ost << "coefficient: " << _coefficient << '\n';
return ost;
}
std::ostream& operator<<(std::ostream& ost, const StairCaseMatrix& matrix) {
ost << (Matrix)matrix << '/' << matrix._coefficient;
return ost;
}
// eof StairCaseMatrix.cc
|
