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
|
//static const char* sccsid = "@(#)math_Householder.cxx 3.2 95/01/10"; // Do not delete this line. Used by sccs.
//#ifndef DEB
#define No_Standard_RangeError
#define No_Standard_OutOfRange
#define No_Standard_DimensionError
//#endif
#include <math_Householder.ixx>
#include <Standard_DimensionError.hxx>
#include <Standard_OutOfRange.hxx>
// Cette classe decrit la methode de Householder qui transforme A en un
// produit de matrice orthogonale par une triangulaire superieure. Les seconds
// membres sont modifies dans le meme temps.
// Les references sur le cote sont celles de l'algorithme explique en page
// 90 du livre "Introduction a l'analyse numerique matricielle et a
// l'optimisation." par P.G. CIARLET, edition MASSON. Les secondes
// references sont celles du sous-programme HOUSEO d'Euclid.
// A la difference du sous-programme Houseo, la premiere colonne n'est pas
// traitee separement. Les tests effectues ont montre que le code effectue
// specialement pour celle-ci etait plus long qu'une simple recopie. C'est
// donc cette solution de recopie initiale qui a ete retenue.
math_Householder::math_Householder(const math_Matrix& A, const math_Vector& B,
const Standard_Real EPS):
Sol(1, A.ColNumber(), 1, 1),
Q(1, A.RowNumber(),
1, A.ColNumber()) {
mylowerArow = A.LowerRow();
mylowerAcol = A.LowerCol();
myupperArow = A.UpperRow();
myupperAcol = A.UpperCol();
math_Matrix B1(1, B.Length(), 1, 1);
B1.SetCol(1, B);
Perform(A, B1, EPS);
}
math_Householder::math_Householder(const math_Matrix& A, const math_Matrix& B,
const Standard_Real EPS):
Sol(1, A.ColNumber(),
1, B.ColNumber()),
Q(1, A.RowNumber(),
A.LowerCol(), A.UpperCol()) {
mylowerArow = A.LowerRow();
mylowerAcol = A.LowerCol();
myupperArow = A.UpperRow();
myupperAcol = A.UpperCol();
Perform(A, B, EPS);
}
math_Householder::math_Householder(const math_Matrix& A, const math_Matrix& B,
const Standard_Integer lowerArow,
const Standard_Integer upperArow,
const Standard_Integer lowerAcol,
const Standard_Integer upperAcol,
const Standard_Real EPS):
Sol(1, upperAcol-lowerAcol+1,
1, B.ColNumber()),
Q(1, upperArow-lowerArow+1,
1, upperAcol-lowerAcol+1) {
mylowerArow = lowerArow;
myupperArow = upperArow;
mylowerAcol = lowerAcol;
myupperAcol = upperAcol;
Perform(A, B, EPS);
}
void math_Householder::Perform(const math_Matrix& A, const math_Matrix& B,
const Standard_Real EPS) {
Standard_Integer i, j, k, n, l, m;
Standard_Real scale, f, g, h = 0., alfaii;
Standard_Real qki;
Standard_Real cj;
n = Q.ColNumber();
l = Q.RowNumber();
m = B.ColNumber();
math_Matrix B2(1, l, 1, m);
Standard_Integer lbrow = B.LowerRow();
for (i = 1; i <= l; i++) {
for (j = 1; j <= n; j++) {
Q(i, j) = A(i+mylowerArow-1, j+mylowerAcol-1);
}
for (j=1; j <= m; j++) {
B2(i, j) = B(i+lbrow-1, j);
}
}
Standard_DimensionError_Raise_if(l != B.RowNumber() || n > l, " ");
// Traitement de chaque colonne de A:
for (i = 1; i <= n; i++) {
h = scale = 0.0;
for (k = i; k <= l; k++) {
qki = Q(k, i);
h += qki*qki; // = ||a||*||a|| = EUAI
}
f = Q(i,i); // = a1 = AII
g = f<1.e-15 ? Sqrt(h) : -Sqrt(h);
if (fabs(g) <= EPS) {
Done = Standard_False;
return;
}
h -= f*g; // = (v*v)/2 = C1
alfaii = g-f; // = v = ALFAII
for (j =i+1; j <= n; j++) {
scale = 0.0;
for (k = i; k <= l; k++) {
scale += Q(k,i)* Q(k,j); // = SCAL
}
cj = (g*Q(i,j) - scale)/h;
Q(i,j) = Q(i,j) - alfaii*cj;
for(k= i+1; k <= l; k++) {
Q(k,j) = Q(k, j) + cj * Q(k,i);
}
}
// Modification de B:
for (j = 1; j <= m; j++) {
scale= Q(i,i)*B2(i,j);
for (k = i+1; k <= l; k++) {
scale += Q(k,i)*B2(k,j);
}
cj = (g*B2(i,j) - scale)/h;
B2(i,j) = B2(i,j) - cj*alfaii;
for (k = i+1; k <= l; k++) {
B2(k,j) = B2(k,j) + cj*Q(k,i);
}
}
Q(i,i) = g;
}
// Remontee:
for (j = 1; j <= m; j++) {
Sol(n,j) = B2(n,j)/Q(n,n);
for (i = n -1; i >=1; i--) {
scale= 0.0;
for(k = i+1; k <= n; k++) {
scale += Q(i,k) * Sol(k,j);
}
Sol(i,j) = (B2(i,j) - scale)/Q(i,i);
}
}
Done = Standard_True;
}
void math_Householder::Dump(Standard_OStream& o) const {
o <<"math_Householder ";
if (Done) {
o << " Status = Done \n";
}
else {
o << "Status = not Done \n";
}
}
|
