-
Notifications
You must be signed in to change notification settings - Fork 61
/
rwupdt_.c
161 lines (124 loc) · 4.51 KB
/
rwupdt_.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
/* rwupdt.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__(rwupdt)(const int *n, real *r__, const int *ldr,
const real *w, real *b, real *alpha, real *cos__,
real *sin__)
{
/* Initialized data */
#define p5 ((real).5)
#define p25 ((real).25)
/* System generated locals */
int r_dim1, r_offset, i__1, i__2;
real d__1;
/* Local variables */
int i__, j, jm1;
real tan__, temp, rowj, cotan;
/* ********** */
/* subroutine rwupdt */
/* given an n by n upper triangular matrix r, this subroutine */
/* computes the qr decomposition of the matrix formed when a row */
/* is added to r. if the row is specified by the vector w, then */
/* rwupdt determines an orthogonal matrix q such that when the */
/* n+1 by n matrix composed of r augmented by w is premultiplied */
/* by (q transpose), the resulting matrix is upper trapezoidal. */
/* the matrix (q transpose) is the product of n transformations */
/* g(n)*g(n-1)* ... *g(1) */
/* where g(i) is a givens rotation in the (i,n+1) plane which */
/* eliminates elements in the (n+1)-st plane. rwupdt also */
/* computes the product (q transpose)*c where c is the */
/* (n+1)-vector (b,alpha). q itself is not accumulated, rather */
/* the information to recover the g rotations is supplied. */
/* the subroutine statement is */
/* subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin) */
/* where */
/* n is a positive integer input variable set to the order of r. */
/* r is an n by n array. on input the upper triangular part of */
/* r must contain the matrix to be updated. on output r */
/* contains the updated triangular matrix. */
/* ldr is a positive integer input variable not less than n */
/* which specifies the leading dimension of the array r. */
/* w is an input array of length n which must contain the row */
/* vector to be added to r. */
/* b is an array of length n. on input b must contain the */
/* first n elements of the vector c. on output b contains */
/* the first n elements of the vector (q transpose)*c. */
/* alpha is a variable. on input alpha must contain the */
/* (n+1)-st element of the vector c. on output alpha contains */
/* the (n+1)-st element of the vector (q transpose)*c. */
/* cos is an output array of length n which contains the */
/* cosines of the transforming givens rotations. */
/* sin is an output array of length n which contains the */
/* sines of the transforming givens rotations. */
/* subprograms called */
/* fortran-supplied ... dabs,dsqrt */
/* argonne national laboratory. minpack project. march 1980. */
/* burton s. garbow, dudley v. goetschel, kenneth e. hillstrom, */
/* jorge j. more */
/* ********** */
/* Parameter adjustments */
--sin__;
--cos__;
--b;
--w;
r_dim1 = *ldr;
r_offset = 1 + r_dim1 * 1;
r__ -= r_offset;
/* Function Body */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
rowj = w[j];
jm1 = j - 1;
/* apply the previous transformations to */
/* r(i,j), i=1,2,...,j-1, and to w(j). */
if (jm1 < 1) {
goto L20;
}
i__2 = jm1;
for (i__ = 1; i__ <= i__2; ++i__) {
temp = cos__[i__] * r__[i__ + j * r_dim1] + sin__[i__] * rowj;
rowj = -sin__[i__] * r__[i__ + j * r_dim1] + cos__[i__] * rowj;
r__[i__ + j * r_dim1] = temp;
/* L10: */
}
L20:
/* determine a givens rotation which eliminates w(j). */
cos__[j] = 1.;
sin__[j] = 0.;
if (rowj == 0.) {
goto L50;
}
if ((d__1 = r__[j + j * r_dim1], abs(d__1)) >= abs(rowj)) {
goto L30;
}
cotan = r__[j + j * r_dim1] / rowj;
/* Computing 2nd power */
d__1 = cotan;
sin__[j] = p5 / sqrt(p25 + p25 * (d__1 * d__1));
cos__[j] = sin__[j] * cotan;
goto L40;
L30:
tan__ = rowj / r__[j + j * r_dim1];
/* Computing 2nd power */
d__1 = tan__;
cos__[j] = p5 / sqrt(p25 + p25 * (d__1 * d__1));
sin__[j] = cos__[j] * tan__;
L40:
/* apply the current transformation to r(j,j), b(j), and alpha. */
r__[j + j * r_dim1] = cos__[j] * r__[j + j * r_dim1] + sin__[j] *
rowj;
temp = cos__[j] * b[j] + sin__[j] * *alpha;
*alpha = -sin__[j] * b[j] + cos__[j] * *alpha;
b[j] = temp;
L50:
/* L60: */
;
}
return;
/* last card of subroutine rwupdt. */
} /* rwupdt_ */