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
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
|
/* r1updt.f -- translated by f2c (version 20020621).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "minpack.h"
#include <math.h>
#include "minpackP.h"
__minpack_attr__
void __minpack_func__(r1updt)(const int *m, const int *n, real *s, const int *
ls, const real *u, real *v, real *w, int *sing)
{
/* Initialized data */
#define p5 .5
#define p25 .25
const int c__3 = 3;
/* System generated locals */
int i__1, i__2;
real d__1, d__2;
/* 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 = __minpack_func__(dpmpar)(&c__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;
i__1 = *m;
for (i__ = *n; i__ <= i__1; ++i__) {
w[i__] = s[l];
++l;
/* L10: */
}
/* 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) {
goto L70;
}
i__1 = nm1;
for (nmj = 1; nmj <= i__1; ++nmj) {
j = *n - nmj;
jj -= *m - j + 1;
w[j] = 0.;
if (v[j] == 0.) {
goto L50;
}
/* determine a givens rotation which eliminates the */
/* j-th element of v. */
if ((d__1 = v[*n], abs(d__1)) >= (d__2 = v[j], abs(d__2))) {
goto L20;
}
cotan = v[*n] / v[j];
/* Computing 2nd power */
d__1 = cotan;
sin__ = p5 / sqrt(p25 + p25 * (d__1 * d__1));
cos__ = sin__ * cotan;
tau = 1.;
if (abs(cos__) * giant > 1.) {
tau = 1. / cos__;
}
goto L30;
L20:
tan__ = v[j] / v[*n];
/* Computing 2nd power */
d__1 = tan__;
cos__ = p5 / sqrt(p25 + p25 * (d__1 * d__1));
sin__ = cos__ * tan__;
tau = sin__;
L30:
/* 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;
i__2 = *m;
for (i__ = j; i__ <= i__2; ++i__) {
temp = cos__ * s[l] - sin__ * w[i__];
w[i__] = sin__ * s[l] + cos__ * w[i__];
s[l] = temp;
++l;
/* L40: */
}
L50:
/* L60: */
;
}
L70:
/* add the spike from the rank 1 update to w. */
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
w[i__] += v[*n] * u[i__];
/* L80: */
}
/* eliminate the spike. */
*sing = FALSE_;
if (nm1 < 1) {
goto L140;
}
i__1 = nm1;
for (j = 1; j <= i__1; ++j) {
if (w[j] == 0.) {
goto L120;
}
/* determine a givens rotation which eliminates the */
/* j-th element of the spike. */
if ((d__1 = s[jj], abs(d__1)) >= (d__2 = w[j], abs(d__2))) {
goto L90;
}
cotan = s[jj] / w[j];
/* Computing 2nd power */
d__1 = cotan;
sin__ = p5 / sqrt(p25 + p25 * (d__1 * d__1));
cos__ = sin__ * cotan;
tau = 1.;
if (abs(cos__) * giant > 1.) {
tau = 1. / cos__;
}
goto L100;
L90:
tan__ = w[j] / s[jj];
/* Computing 2nd power */
d__1 = tan__;
cos__ = p5 / sqrt(p25 + p25 * (d__1 * d__1));
sin__ = cos__ * tan__;
tau = sin__;
L100:
/* apply the transformation to s and reduce the spike in w. */
l = jj;
i__2 = *m;
for (i__ = j; i__ <= i__2; ++i__) {
temp = cos__ * s[l] + sin__ * w[i__];
w[i__] = -sin__ * s[l] + cos__ * w[i__];
s[l] = temp;
++l;
/* L110: */
}
/* store the information necessary to recover the */
/* givens rotation. */
w[j] = tau;
L120:
/* test for zero diagonal elements in the output s. */
if (s[jj] == 0.) {
*sing = TRUE_;
}
jj += *m - j + 1;
/* L130: */
}
L140:
/* move w back into the last column of the output s. */
l = jj;
i__1 = *m;
for (i__ = *n; i__ <= i__1; ++i__) {
s[l] = w[i__];
++l;
/* L150: */
}
if (s[jj] == 0.) {
*sing = TRUE_;
}
return;
/* last card of subroutine r1updt. */
} /* r1updt_ */
|
