6663b6c9
 adorian
projet complet av...
|
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
|
/* -*- mode:C++ ; compile-command: "g++-3.4 -I.. -I../include -g -c -Wall TmpLESystemSolver.C" -*- */
// Copyright (c) 2006 Stefan Kaspar
// This file is part of the source of CoCoALib, the CoCoA Library.
// CoCoALib is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License (version 3)
// as published by the Free Software Foundation. A copy of the full
// licence may be found in the file COPYING in this directory.
// CoCoALib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with CoCoA; if not, see <http://www.gnu.org/licenses/>.
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "first.h"
#ifdef USE_GMP_REPLACEMENTS
#undef HAVE_LIBCOCOA
#endif
#ifdef HAVE_LIBCOCOA
#include "TmpLESystemSolver.H"
#include "CoCoA/DenseMatrix.H"
#include "CoCoA/matrix.H"
#include "CoCoA/ring.H"
#include "CoCoA/error.H"
// #include <vector> // Included by DenseMatrix.H
using std::vector;
#include <utility>
using std::pair;
using std::make_pair;
// #include <cstddef> // Included by DenseMatrix.H
using std::size_t;
namespace CoCoADortmund
{
using namespace CoCoA;
// Create a copy of a matrix
void CopyMatrix(matrix& MTarget, const matrix& MSource)
{
const size_t NumRowsMSource = NumRows(MSource);
const size_t NumColsMSource = NumCols(MSource);
for (size_t row = 0; row < NumRowsMSource; ++row)
for (size_t col = 0; col < NumColsMSource; ++col)
SetEntry(MTarget, row, col, MSource(row, col));
}
// Solve the linear system M*x = b by using Gauss' algorithm
bool LESystemSolver(matrix& x0, const matrix& M, const matrix& b)
{
const size_t NumRowsM = NumRows(M);
const size_t NumColsM = NumCols(M);
const size_t NumRowsb = NumRows(b);
const size_t NumColsb = NumCols(b);
// Dimension check
if (NumRowsM != NumRowsb)
CoCoA_ERROR(ERR::BadMatrixSize, "mySolve: M and b must have same number of rows.");
if (NumColsM != NumRows(x0))
CoCoA_ERROR(ERR::BadMatrixSize, "mySolve: M and x0 must have same number of columns.");
if (NumCols(x0) != 1)
CoCoA_ERROR(ERR::BadMatrixSize, "mySolve: NumCols(x0) > 1.");
if (NumColsb != 1)
CoCoA_ERROR(ERR::BadMatrixSize, "mySolve: NumCols(b) > 1.");
// Field check; should we also check if BaseRing(M) = BaseRing(b) = BaseRing(x0)?
ring K(BaseRing(M));
if (!IsField(K))
CoCoA_ERROR(ERR::NotField, "mySolve: Gauss' algorithm over non-fields not yet implemented.");
// Create working copies of M and b
matrix MCopy(NewDenseMat(K, NumRowsM, NumColsM));
CopyMatrix(MCopy, M);
matrix bCopy(NewDenseMat(K, NumRowsb, NumColsb));
CopyMatrix(bCopy, b);
// For solution computation
vector< pair<size_t, size_t> > positions;
// Apply Gauss' algorithm
RingElem c(K);
size_t row = 0;
for (size_t col = 0; col < NumColsM && row < NumRowsM; ++col)
{
// Check if current column contains an element != 0
if (IsZero(MCopy(row, col)))
{
size_t i = row+1;
for ( ; i < NumRowsM; ++i)
{
if (!IsZero(MCopy(i, col)))
{
// Switch MCopy and bCopy rows
MCopy->mySwapRows(i, row);
bCopy->mySwapRows(i, row);
break;
}
}
if (i == NumRowsM)
continue;
}
// For solution computation
positions.push_back(make_pair(row, col));
// Found an element != 0 in current column; apply elemination
c = MCopy(row, col);
for (size_t i = row+1; i < NumRowsM; ++i)
{
// Transform MCopy and bCopy
bCopy->myAddRowMul(i, row, -MCopy(i, col)/c);
MCopy->myAddRowMul(i, row, -MCopy(i, col)/c);
}
++row;
}
// row = rank(MCopy); check if a solution for the equation system exists
for (size_t i = row; i < NumRowsb; ++i)
{
if (!IsZero(bCopy(i, 0)))
return false;
}
// Compute components (x_1, ..., x_n) of vector x0 backwards from x_n to x_1,
// possibly skipping some components
matrix x0Tmp = NewDenseMat(K, NumRows(x0), 1);
while (!positions.empty())
{
const size_t i = positions.back().first, j = positions.back().second;
RingElem x(bCopy(i, 0));
for (size_t k = j + 1; k < NumColsM; ++k)
{
x -= MCopy(i, k) * x0Tmp(k, 0);
}
SetEntry(x0Tmp, j, 0, x/MCopy(i, j));
positions.pop_back();
}
CopyMatrix(x0, x0Tmp);
return true;
}
} // end of namespace CoCoA
#endif
|