-
Notifications
You must be signed in to change notification settings - Fork 63
/
r1updt.c
236 lines (176 loc) · 6.45 KB
/
r1updt.c
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
/* r1updt.f -- translated by f2c (version 20020621).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "cminpack.h"
#include <math.h>
#include "cminpackP.h"
__cminpack_attr__
void __cminpack_func__(r1updt)(int m, int n, real *s, int
ls, const real *u, real *v, real *w, int *sing)
{
/* Initialized data */
#define p5 ((real).5)
#define p25 ((real).25)
/* Local variables */
int i, j, l, jj, nm1;
real tan;
int nmj;
real cos, sin, tau, temp, giant, cotan;
/* ********** */
/* subroutine r1updt */
/* given an m by n lower trapezoidal matrix s, an m-vector u, */
/* and an n-vector v, the problem is to determine an */
/* orthogonal matrix q such that */
/* t */
/* (s + u*v )*q */
/* is again lower trapezoidal. */
/* this subroutine determines q as the product of 2*(n - 1) */
/* transformations */
/* gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1) */
/* where gv(i), gw(i) are givens rotations in the (i,n) plane */
/* which eliminate elements in the i-th and n-th planes, */
/* respectively. q itself is not accumulated, rather the */
/* information to recover the gv, gw rotations is returned. */
/* the subroutine statement is */
/* subroutine r1updt(m,n,s,ls,u,v,w,sing) */
/* where */
/* m is a positive integer input variable set to the number */
/* of rows of s. */
/* n is a positive integer input variable set to the number */
/* of columns of s. n must not exceed m. */
/* s is an array of length ls. on input s must contain the lower */
/* trapezoidal matrix s stored by columns. on output s contains */
/* the lower trapezoidal matrix produced as described above. */
/* ls is a positive integer input variable not less than */
/* (n*(2*m-n+1))/2. */
/* u is an input array of length m which must contain the */
/* vector u. */
/* v is an array of length n. on input v must contain the vector */
/* v. on output v(i) contains the information necessary to */
/* recover the givens rotation gv(i) described above. */
/* w is an output array of length m. w(i) contains information */
/* necessary to recover the givens rotation gw(i) described */
/* above. */
/* sing is a logical output variable. sing is set true if any */
/* of the diagonal elements of the output s are zero. otherwise */
/* sing is set false. */
/* subprograms called */
/* minpack-supplied ... dpmpar */
/* fortran-supplied ... dabs,dsqrt */
/* argonne national laboratory. minpack project. march 1980. */
/* burton s. garbow, kenneth e. hillstrom, jorge j. more, */
/* john l. nazareth */
/* ********** */
/* Parameter adjustments */
--w;
--u;
--v;
--s;
(void)ls;
/* Function Body */
/* giant is the largest magnitude. */
giant = __cminpack_func__(dpmpar)(3);
/* initialize the diagonal element pointer. */
jj = n * ((m << 1) - n + 1) / 2 - (m - n);
/* move the nontrivial part of the last column of s into w. */
l = jj;
for (i = n; i <= m; ++i) {
w[i] = s[l];
++l;
}
/* rotate the vector v into a multiple of the n-th unit vector */
/* in such a way that a spike is introduced into w. */
nm1 = n - 1;
if (nm1 >= 1) {
for (nmj = 1; nmj <= nm1; ++nmj) {
j = n - nmj;
jj -= m - j + 1;
w[j] = 0.;
if (v[j] != 0.) {
/* determine a givens rotation which eliminates the */
/* j-th element of v. */
if (fabs(v[n]) < fabs(v[j])) {
cotan = v[n] / v[j];
sin = p5 / sqrt(p25 + p25 * (cotan * cotan));
cos = sin * cotan;
tau = 1.;
if (fabs(cos) * giant > 1.) {
tau = 1 / cos;
}
} else {
tan = v[j] / v[n];
cos = p5 / sqrt(p25 + p25 * (tan * tan));
sin = cos * tan;
tau = sin;
}
/* apply the transformation to v and store the information */
/* necessary to recover the givens rotation. */
v[n] = sin * v[j] + cos * v[n];
v[j] = tau;
/* apply the transformation to s and extend the spike in w. */
l = jj;
for (i = j; i <= m; ++i) {
temp = cos * s[l] - sin * w[i];
w[i] = sin * s[l] + cos * w[i];
s[l] = temp;
++l;
}
}
}
}
/* add the spike from the rank 1 update to w. */
for (i = 1; i <= m; ++i) {
w[i] += v[n] * u[i];
}
/* eliminate the spike. */
*sing = FALSE_;
if (nm1 >= 1) {
for (j = 1; j <= nm1; ++j) {
if (w[j] != 0.) {
/* determine a givens rotation which eliminates the */
/* j-th element of the spike. */
if (fabs(s[jj]) < fabs(w[j])) {
cotan = s[jj] / w[j];
sin = p5 / sqrt(p25 + p25 * (cotan * cotan));
cos = sin * cotan;
tau = 1.;
if (fabs(cos) * giant > 1.) {
tau = 1 / cos;
}
} else {
tan = w[j] / s[jj];
cos = p5 / sqrt(p25 + p25 * (tan * tan));
sin = cos * tan;
tau = sin;
}
/* apply the transformation to s and reduce the spike in w. */
l = jj;
for (i = j; i <= m; ++i) {
temp = cos * s[l] + sin * w[i];
w[i] = -sin * s[l] + cos * w[i];
s[l] = temp;
++l;
}
/* store the information necessary to recover the */
/* givens rotation. */
w[j] = tau;
}
/* test for zero diagonal elements in the output s. */
if (s[jj] == 0.) {
*sing = TRUE_;
}
jj += m - j + 1;
}
}
/* move w back into the last column of the output s. */
l = jj;
for (i = n; i <= m; ++i) {
s[l] = w[i];
++l;
}
if (s[jj] == 0.) {
*sing = TRUE_;
}
/* last card of subroutine r1updt. */
} /* __minpack_func__(r1updt) */