1 /* ./src_f77/dlarrf.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
dlarrf_(integer * n,doublereal * d__,doublereal * l,doublereal * ld,doublereal * lld,integer * ifirst,integer * ilast,doublereal * w,doublereal * dplus,doublereal * lplus,doublereal * work,integer * iwork,integer * info)8 /* Subroutine */ int dlarrf_(integer *n, doublereal *d__, doublereal *l,
9 doublereal *ld, doublereal *lld, integer *ifirst, integer *ilast,
10 doublereal *w, doublereal *dplus, doublereal *lplus, doublereal *work,
11 integer *iwork, integer *info)
12 {
13 /* System generated locals */
14 integer i__1;
15
16 /* Local variables */
17 static integer i__;
18 static doublereal s, eps, delta, sigma;
19 extern doublereal dlamch_(char *, ftnlen);
20
21
22 /* -- LAPACK auxiliary routine (version 3.0) -- */
23 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
24 /* Courant Institute, Argonne National Lab, and Rice University */
25 /* June 30, 1999 */
26
27 /* .. Scalar Arguments .. */
28 /* .. */
29 /* .. Array Arguments .. */
30 /* .. */
31
32 /* Purpose */
33 /* ======= */
34
35 /* Given the initial representation L D L^T and its cluster of close */
36 /* eigenvalues (in a relative measure), W( IFIRST ), W( IFIRST+1 ), ... */
37 /* W( ILAST ), DLARRF finds a new relatively robust representation */
38 /* L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the */
39 /* eigenvalues of L(+) D(+) L(+)^T is relatively isolated. */
40
41 /* Arguments */
42 /* ========= */
43
44 /* N (input) INTEGER */
45 /* The order of the matrix. */
46
47 /* D (input) DOUBLE PRECISION array, dimension (N) */
48 /* The n diagonal elements of the diagonal matrix D. */
49
50 /* L (input) DOUBLE PRECISION array, dimension (N-1) */
51 /* The (n-1) subdiagonal elements of the unit bidiagonal */
52 /* matrix L. */
53
54 /* LD (input) DOUBLE PRECISION array, dimension (N-1) */
55 /* The n-1 elements L(i)*D(i). */
56
57 /* LLD (input) DOUBLE PRECISION array, dimension (N-1) */
58 /* The n-1 elements L(i)*L(i)*D(i). */
59
60 /* IFIRST (input) INTEGER */
61 /* The index of the first eigenvalue in the cluster. */
62
63 /* ILAST (input) INTEGER */
64 /* The index of the last eigenvalue in the cluster. */
65
66 /* W (input/output) DOUBLE PRECISION array, dimension (N) */
67 /* On input, the eigenvalues of L D L^T in ascending order. */
68 /* W( IFIRST ) through W( ILAST ) form the cluster of relatively */
69 /* close eigenalues. */
70 /* On output, W( IFIRST ) thru' W( ILAST ) are estimates of the */
71 /* corresponding eigenvalues of L(+) D(+) L(+)^T. */
72
73 /* SIGMA (input) DOUBLE PRECISION */
74 /* The shift used to form L(+) D(+) L(+)^T. */
75
76 /* DPLUS (output) DOUBLE PRECISION array, dimension (N) */
77 /* The n diagonal elements of the diagonal matrix D(+). */
78
79 /* LPLUS (output) DOUBLE PRECISION array, dimension (N) */
80 /* The first (n-1) elements of LPLUS contain the subdiagonal */
81 /* elements of the unit bidiagonal matrix L(+). LPLUS( N ) is */
82 /* set to SIGMA. */
83
84 /* WORK (input) DOUBLE PRECISION array, dimension (???) */
85 /* Workspace. */
86
87 /* Further Details */
88 /* =============== */
89
90 /* Based on contributions by */
91 /* Inderjit Dhillon, IBM Almaden, USA */
92 /* Osni Marques, LBNL/NERSC, USA */
93
94 /* ===================================================================== */
95
96 /* .. Parameters .. */
97 /* .. */
98 /* .. Local Scalars .. */
99 /* .. */
100 /* .. External Functions .. */
101 /* .. */
102 /* .. Intrinsic Functions .. */
103 /* .. */
104 /* .. Executable Statements .. */
105
106 /* Parameter adjustments */
107 --iwork;
108 --work;
109 --lplus;
110 --dplus;
111 --w;
112 --lld;
113 --ld;
114 --l;
115 --d__;
116
117 /* Function Body */
118 *info = 0;
119 eps = dlamch_("Precision", (ftnlen)9);
120 if (*ifirst == 1) {
121 sigma = w[*ifirst];
122 } else if (*ilast == *n) {
123 sigma = w[*ilast];
124 } else {
125 *info = 1;
126 return 0;
127 }
128
129 /* Compute the new relatively robust representation (RRR) */
130
131 delta = eps * 2.;
132 L10:
133 if (*ifirst == 1) {
134 sigma -= abs(sigma) * delta;
135 } else {
136 sigma += abs(sigma) * delta;
137 }
138 s = -sigma;
139 i__1 = *n - 1;
140 for (i__ = 1; i__ <= i__1; ++i__) {
141 dplus[i__] = d__[i__] + s;
142 lplus[i__] = ld[i__] / dplus[i__];
143 s = s * lplus[i__] * l[i__] - sigma;
144 /* L20: */
145 }
146 dplus[*n] = d__[*n] + s;
147 if (*ifirst == 1) {
148 i__1 = *n;
149 for (i__ = 1; i__ <= i__1; ++i__) {
150 if (dplus[i__] < 0.) {
151 delta *= 2.;
152 goto L10;
153 }
154 /* L30: */
155 }
156 } else {
157 i__1 = *n;
158 for (i__ = 1; i__ <= i__1; ++i__) {
159 if (dplus[i__] > 0.) {
160 delta *= 2.;
161 goto L10;
162 }
163 /* L40: */
164 }
165 }
166 i__1 = *ilast;
167 for (i__ = *ifirst; i__ <= i__1; ++i__) {
168 w[i__] -= sigma;
169 /* L50: */
170 }
171 lplus[*n] = sigma;
172
173 return 0;
174
175 /* End of DLARRF */
176
177 } /* dlarrf_ */
178
179