1 /* rwupdt.f -- translated by f2c (version 20020621).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "cminpack.h"
7 #include <math.h>
8 #include "cminpackP.h"
9
10 __cminpack_attr__
__cminpack_func__(rwupdt)11 void __cminpack_func__(rwupdt)(int n, real *r, int ldr,
12 const real *w, real *b, real *alpha, real *cos,
13 real *sin)
14 {
15 /* Initialized data */
16
17 #define p5 .5
18 #define p25 .25
19
20 /* System generated locals */
21 int r_dim1, r_offset;
22
23 /* Local variables */
24 int i, j, jm1;
25 real tan, temp, rowj, cotan;
26
27 /* ********** */
28
29 /* subroutine rwupdt */
30
31 /* given an n by n upper triangular matrix r, this subroutine */
32 /* computes the qr decomposition of the matrix formed when a row */
33 /* is added to r. if the row is specified by the vector w, then */
34 /* rwupdt determines an orthogonal matrix q such that when the */
35 /* n+1 by n matrix composed of r augmented by w is premultiplied */
36 /* by (q transpose), the resulting matrix is upper trapezoidal. */
37 /* the matrix (q transpose) is the product of n transformations */
38
39 /* g(n)*g(n-1)* ... *g(1) */
40
41 /* where g(i) is a givens rotation in the (i,n+1) plane which */
42 /* eliminates elements in the (n+1)-st plane. rwupdt also */
43 /* computes the product (q transpose)*c where c is the */
44 /* (n+1)-vector (b,alpha). q itself is not accumulated, rather */
45 /* the information to recover the g rotations is supplied. */
46
47 /* the subroutine statement is */
48
49 /* subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin) */
50
51 /* where */
52
53 /* n is a positive integer input variable set to the order of r. */
54
55 /* r is an n by n array. on input the upper triangular part of */
56 /* r must contain the matrix to be updated. on output r */
57 /* contains the updated triangular matrix. */
58
59 /* ldr is a positive integer input variable not less than n */
60 /* which specifies the leading dimension of the array r. */
61
62 /* w is an input array of length n which must contain the row */
63 /* vector to be added to r. */
64
65 /* b is an array of length n. on input b must contain the */
66 /* first n elements of the vector c. on output b contains */
67 /* the first n elements of the vector (q transpose)*c. */
68
69 /* alpha is a variable. on input alpha must contain the */
70 /* (n+1)-st element of the vector c. on output alpha contains */
71 /* the (n+1)-st element of the vector (q transpose)*c. */
72
73 /* cos is an output array of length n which contains the */
74 /* cosines of the transforming givens rotations. */
75
76 /* sin is an output array of length n which contains the */
77 /* sines of the transforming givens rotations. */
78
79 /* subprograms called */
80
81 /* fortran-supplied ... dabs,dsqrt */
82
83 /* argonne national laboratory. minpack project. march 1980. */
84 /* burton s. garbow, dudley v. goetschel, kenneth e. hillstrom, */
85 /* jorge j. more */
86
87 /* ********** */
88 /* Parameter adjustments */
89 --sin;
90 --cos;
91 --b;
92 --w;
93 r_dim1 = ldr;
94 r_offset = 1 + r_dim1 * 1;
95 r -= r_offset;
96
97 /* Function Body */
98
99 for (j = 1; j <= n; ++j) {
100 rowj = w[j];
101 jm1 = j - 1;
102
103 /* apply the previous transformations to */
104 /* r(i,j), i=1,2,...,j-1, and to w(j). */
105
106 if (jm1 >= 1) {
107 for (i = 1; i <= jm1; ++i) {
108 temp = cos[i] * r[i + j * r_dim1] + sin[i] * rowj;
109 rowj = -sin[i] * r[i + j * r_dim1] + cos[i] * rowj;
110 r[i + j * r_dim1] = temp;
111 }
112 }
113
114 /* determine a givens rotation which eliminates w(j). */
115
116 cos[j] = 1.;
117 sin[j] = 0.;
118 if (rowj != 0.) {
119 if (fabs(r[j + j * r_dim1]) < fabs(rowj)) {
120 cotan = r[j + j * r_dim1] / rowj;
121 sin[j] = p5 / sqrt(p25 + p25 * (cotan * cotan));
122 cos[j] = sin[j] * cotan;
123 } else {
124 tan = rowj / r[j + j * r_dim1];
125 cos[j] = p5 / sqrt(p25 + p25 * (tan * tan));
126 sin[j] = cos[j] * tan;
127 }
128
129 /* apply the current transformation to r(j,j), b(j), and alpha. */
130
131 r[j + j * r_dim1] = cos[j] * r[j + j * r_dim1] + sin[j] * rowj;
132 temp = cos[j] * b[j] + sin[j] * *alpha;
133 *alpha = -sin[j] * b[j] + cos[j] * *alpha;
134 b[j] = temp;
135 }
136 }
137
138 /* last card of subroutine rwupdt. */
139
140 } /* rwupdt_ */
141
142