1: #region Translated by Jose Antonio De Santiago-Castillo.
2:
3: //Translated by Jose Antonio De Santiago-Castillo.
4: //E-mail:JAntonioDeSantiago@gmail.com
5: //Web: www.DotNumerics.com
6: //
7: //Fortran to C# Translation.
8: //Translated by:
9: //F2CSharp Version 0.71 (November 10, 2009)
10: //Code Optimizations: None
11: //
12: #endregion
13:
14: using System;
15: using DotNumerics.FortranLibrary;
16:
17: namespace DotNumerics.CSLapack
18: {
19: /// <summary>
20: /// -- LAPACK routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// This subroutine computes the I-th updated eigenvalue of a symmetric
27: /// rank-one modification to a diagonal matrix whose elements are
28: /// given in the array d, and that
29: ///
30: /// D(i) .LT. D(j) for i .LT. j
31: ///
32: /// and that RHO .GT. 0. This is arranged by the calling routine, and is
33: /// no loss in generality. The rank-one modified system is thus
34: ///
35: /// diag( D ) + RHO * Z * Z_transpose.
36: ///
37: /// where we assume the Euclidean norm of Z is 1.
38: ///
39: /// The method consists of approximating the rational functions in the
40: /// secular equation by simpler interpolating rational functions.
41: ///
42: ///</summary>
43: public class DLAED4
44: {
45:
46:
47: #region Dependencies
48:
49: DLAMCH _dlamch; DLAED5 _dlaed5; DLAED6 _dlaed6;
50:
51: #endregion
52:
53:
54: #region Fields
55:
56: const int MAXIT = 30; const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0;
57: const double THREE = 3.0E0;const double FOUR = 4.0E0; const double EIGHT = 8.0E0; const double TEN = 10.0E0;
58: bool ORGATI = false;bool SWTCH = false; bool SWTCH3 = false; int II = 0; int IIM1 = 0; int IIP1 = 0; int IP1 = 0;
59: int ITER = 0;int J = 0; int NITER = 0; double A = 0; double B = 0; double C = 0; double DEL = 0; double DLTLB = 0;
60: double DLTUB = 0;double DPHI = 0; double DPSI = 0; double DW = 0; double EPS = 0; double ERRETM = 0; double ETA = 0;
61: double MIDPT = 0;double PHI = 0; double PREW = 0; double PSI = 0; double RHOINV = 0; double TAU = 0; double TEMP = 0;
62: double TEMP1 = 0;double W = 0; double[] ZZ = new double[3]; int offset_zz = 0; int o_zz = -1;
63:
64: #endregion
65:
66: public DLAED4(DLAMCH dlamch, DLAED5 dlaed5, DLAED6 dlaed6)
67: {
68:
69:
70: #region Set Dependencies
71:
72: this._dlamch = dlamch; this._dlaed5 = dlaed5; this._dlaed6 = dlaed6;
73:
74: #endregion
75:
76: }
77:
78: public DLAED4()
79: {
80:
81:
82: #region Dependencies (Initialization)
83:
84: LSAME lsame = new LSAME();
85: DLAMC3 dlamc3 = new DLAMC3();
86: DLAED5 dlaed5 = new DLAED5();
87: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
88: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
89: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
90: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
91: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
92: DLAED6 dlaed6 = new DLAED6(dlamch);
93:
94: #endregion
95:
96:
97: #region Set Dependencies
98:
99: this._dlamch = dlamch; this._dlaed5 = dlaed5; this._dlaed6 = dlaed6;
100:
101: #endregion
102:
103: }
104: /// <summary>
105: /// Purpose
106: /// =======
107: ///
108: /// This subroutine computes the I-th updated eigenvalue of a symmetric
109: /// rank-one modification to a diagonal matrix whose elements are
110: /// given in the array d, and that
111: ///
112: /// D(i) .LT. D(j) for i .LT. j
113: ///
114: /// and that RHO .GT. 0. This is arranged by the calling routine, and is
115: /// no loss in generality. The rank-one modified system is thus
116: ///
117: /// diag( D ) + RHO * Z * Z_transpose.
118: ///
119: /// where we assume the Euclidean norm of Z is 1.
120: ///
121: /// The method consists of approximating the rational functions in the
122: /// secular equation by simpler interpolating rational functions.
123: ///
124: ///</summary>
125: /// <param name="N">
126: /// (input) INTEGER
127: /// The length of all arrays.
128: ///</param>
129: /// <param name="I">
130: /// (input) INTEGER
131: /// The index of the eigenvalue to be computed. 1 .LE. I .LE. N.
132: ///</param>
133: /// <param name="D">
134: /// (input) DOUBLE PRECISION array, dimension (N)
135: /// The original eigenvalues. It is assumed that they are in
136: /// order, D(I) .LT. D(J) for I .LT. J.
137: ///</param>
138: /// <param name="Z">
139: /// (input) DOUBLE PRECISION array, dimension (N)
140: /// The components of the updating vector.
141: ///</param>
142: /// <param name="DELTA">
143: /// (output) DOUBLE PRECISION array, dimension (N)
144: /// If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th
145: /// component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5
146: /// for detail. The vector DELTA contains the information necessary
147: /// to construct the eigenvectors by DLAED3 and DLAED9.
148: ///</param>
149: /// <param name="RHO">
150: /// (input) DOUBLE PRECISION
151: /// The scalar in the symmetric updating formula.
152: ///</param>
153: /// <param name="DLAM">
154: /// (output) DOUBLE PRECISION
155: /// The computed lambda_I, the I-th updated eigenvalue.
156: ///</param>
157: /// <param name="INFO">
158: /// (output) INTEGER
159: /// = 0: successful exit
160: /// .GT. 0: if INFO = 1, the updating process failed.
161: ///</param>
162: public void Run(int N, int I, double[] D, int offset_d, double[] Z, int offset_z, ref double[] DELTA, int offset_delta, double RHO
163: , ref double DLAM, ref int INFO)
164: {
165:
166: #region Array Index Correction
167:
168: int o_d = -1 + offset_d; int o_z = -1 + offset_z; int o_delta = -1 + offset_delta;
169:
170: #endregion
171:
172:
173: #region Prolog
174:
175: // *
176: // * -- LAPACK routine (version 3.1) --
177: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
178: // * November 2006
179: // *
180: // * .. Scalar Arguments ..
181: // * ..
182: // * .. Array Arguments ..
183: // * ..
184: // *
185: // * Purpose
186: // * =======
187: // *
188: // * This subroutine computes the I-th updated eigenvalue of a symmetric
189: // * rank-one modification to a diagonal matrix whose elements are
190: // * given in the array d, and that
191: // *
192: // * D(i) < D(j) for i < j
193: // *
194: // * and that RHO > 0. This is arranged by the calling routine, and is
195: // * no loss in generality. The rank-one modified system is thus
196: // *
197: // * diag( D ) + RHO * Z * Z_transpose.
198: // *
199: // * where we assume the Euclidean norm of Z is 1.
200: // *
201: // * The method consists of approximating the rational functions in the
202: // * secular equation by simpler interpolating rational functions.
203: // *
204: // * Arguments
205: // * =========
206: // *
207: // * N (input) INTEGER
208: // * The length of all arrays.
209: // *
210: // * I (input) INTEGER
211: // * The index of the eigenvalue to be computed. 1 <= I <= N.
212: // *
213: // * D (input) DOUBLE PRECISION array, dimension (N)
214: // * The original eigenvalues. It is assumed that they are in
215: // * order, D(I) < D(J) for I < J.
216: // *
217: // * Z (input) DOUBLE PRECISION array, dimension (N)
218: // * The components of the updating vector.
219: // *
220: // * DELTA (output) DOUBLE PRECISION array, dimension (N)
221: // * If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th
222: // * component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5
223: // * for detail. The vector DELTA contains the information necessary
224: // * to construct the eigenvectors by DLAED3 and DLAED9.
225: // *
226: // * RHO (input) DOUBLE PRECISION
227: // * The scalar in the symmetric updating formula.
228: // *
229: // * DLAM (output) DOUBLE PRECISION
230: // * The computed lambda_I, the I-th updated eigenvalue.
231: // *
232: // * INFO (output) INTEGER
233: // * = 0: successful exit
234: // * > 0: if INFO = 1, the updating process failed.
235: // *
236: // * Internal Parameters
237: // * ===================
238: // *
239: // * Logical variable ORGATI (origin-at-i?) is used for distinguishing
240: // * whether D(i) or D(i+1) is treated as the origin.
241: // *
242: // * ORGATI = .true. origin at i
243: // * ORGATI = .false. origin at i+1
244: // *
245: // * Logical variable SWTCH3 (switch-for-3-poles?) is for noting
246: // * if we are working with THREE poles!
247: // *
248: // * MAXIT is the maximum number of iterations allowed for each
249: // * eigenvalue.
250: // *
251: // * Further Details
252: // * ===============
253: // *
254: // * Based on contributions by
255: // * Ren-Cang Li, Computer Science Division, University of California
256: // * at Berkeley, USA
257: // *
258: // * =====================================================================
259: // *
260: // * .. Parameters ..
261: // * ..
262: // * .. Local Scalars ..
263: // * ..
264: // * .. Local Arrays ..
265: // * ..
266: // * .. External Functions ..
267: // * ..
268: // * .. External Subroutines ..
269: // * ..
270: // * .. Intrinsic Functions ..
271: // INTRINSIC ABS, MAX, MIN, SQRT;
272: // * ..
273: // * .. Executable Statements ..
274: // *
275: // * Since this routine is called in an inner loop, we do no argument
276: // * checking.
277: // *
278: // * Quick return for N=1 and 2.
279: // *
280:
281: #endregion
282:
283:
284: #region Body
285:
286: INFO = 0;
287: if (N == 1)
288: {
289: // *
290: // * Presumably, I=1 upon entry
291: // *
292: DLAM = D[1 + o_d] + RHO * Z[1 + o_z] * Z[1 + o_z];
293: DELTA[1 + o_delta] = ONE;
294: return;
295: }
296: if (N == 2)
297: {
298: this._dlaed5.Run(I, D, offset_d, Z, offset_z, ref DELTA, offset_delta, RHO, ref DLAM);
299: return;
300: }
301: // *
302: // * Compute machine epsilon
303: // *
304: EPS = this._dlamch.Run("Epsilon");
305: RHOINV = ONE / RHO;
306: // *
307: // * The case I = N
308: // *
309: if (I == N)
310: {
311: // *
312: // * Initialize some basic variables
313: // *
314: II = N - 1;
315: NITER = 1;
316: // *
317: // * Calculate initial guess
318: // *
319: MIDPT = RHO / TWO;
320: // *
321: // * If ||Z||_2 is not one, then TEMP should be set to
322: // * RHO * ||Z||_2^2 / TWO
323: // *
324: for (J = 1; J <= N; J++)
325: {
326: DELTA[J + o_delta] = (D[J + o_d] - D[I + o_d]) - MIDPT;
327: }
328: // *
329: PSI = ZERO;
330: for (J = 1; J <= N - 2; J++)
331: {
332: PSI = PSI + Z[J + o_z] * Z[J + o_z] / DELTA[J + o_delta];
333: }
334: // *
335: C = RHOINV + PSI;
336: W = C + Z[II + o_z] * Z[II + o_z] / DELTA[II + o_delta] + Z[N + o_z] * Z[N + o_z] / DELTA[N + o_delta];
337: // *
338: if (W <= ZERO)
339: {
340: TEMP = Z[N - 1 + o_z] * Z[N - 1 + o_z] / (D[N + o_d] - D[N - 1 + o_d] + RHO) + Z[N + o_z] * Z[N + o_z] / RHO;
341: if (C <= TEMP)
342: {
343: TAU = RHO;
344: }
345: else
346: {
347: DEL = D[N + o_d] - D[N - 1 + o_d];
348: A = - C * DEL + Z[N - 1 + o_z] * Z[N - 1 + o_z] + Z[N + o_z] * Z[N + o_z];
349: B = Z[N + o_z] * Z[N + o_z] * DEL;
350: if (A < ZERO)
351: {
352: TAU = TWO * B / (Math.Sqrt(A * A + FOUR * B * C) - A);
353: }
354: else
355: {
356: TAU = (A + Math.Sqrt(A * A + FOUR * B * C)) / (TWO * C);
357: }
358: }
359: // *
360: // * It can be proved that
361: // * D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO
362: // *
363: DLTLB = MIDPT;
364: DLTUB = RHO;
365: }
366: else
367: {
368: DEL = D[N + o_d] - D[N - 1 + o_d];
369: A = - C * DEL + Z[N - 1 + o_z] * Z[N - 1 + o_z] + Z[N + o_z] * Z[N + o_z];
370: B = Z[N + o_z] * Z[N + o_z] * DEL;
371: if (A < ZERO)
372: {
373: TAU = TWO * B / (Math.Sqrt(A * A + FOUR * B * C) - A);
374: }
375: else
376: {
377: TAU = (A + Math.Sqrt(A * A + FOUR * B * C)) / (TWO * C);
378: }
379: // *
380: // * It can be proved that
381: // * D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2
382: // *
383: DLTLB = ZERO;
384: DLTUB = MIDPT;
385: }
386: // *
387: for (J = 1; J <= N; J++)
388: {
389: DELTA[J + o_delta] = (D[J + o_d] - D[I + o_d]) - TAU;
390: }
391: // *
392: // * Evaluate PSI and the derivative DPSI
393: // *
394: DPSI = ZERO;
395: PSI = ZERO;
396: ERRETM = ZERO;
397: for (J = 1; J <= II; J++)
398: {
399: TEMP = Z[J + o_z] / DELTA[J + o_delta];
400: PSI = PSI + Z[J + o_z] * TEMP;
401: DPSI = DPSI + TEMP * TEMP;
402: ERRETM = ERRETM + PSI;
403: }
404: ERRETM = Math.Abs(ERRETM);
405: // *
406: // * Evaluate PHI and the derivative DPHI
407: // *
408: TEMP = Z[N + o_z] / DELTA[N + o_delta];
409: PHI = Z[N + o_z] * TEMP;
410: DPHI = TEMP * TEMP;
411: ERRETM = EIGHT * ( - PHI - PSI) + ERRETM - PHI + RHOINV + Math.Abs(TAU) * (DPSI + DPHI);
412: // *
413: W = RHOINV + PHI + PSI;
414: // *
415: // * Test for convergence
416: // *
417: if (Math.Abs(W) <= EPS * ERRETM)
418: {
419: DLAM = D[I + o_d] + TAU;
420: goto LABEL250;
421: }
422: // *
423: if (W <= ZERO)
424: {
425: DLTLB = Math.Max(DLTLB, TAU);
426: }
427: else
428: {
429: DLTUB = Math.Min(DLTUB, TAU);
430: }
431: // *
432: // * Calculate the new step
433: // *
434: NITER = NITER + 1;
435: C = W - DELTA[N - 1 + o_delta] * DPSI - DELTA[N + o_delta] * DPHI;
436: A = (DELTA[N - 1 + o_delta] + DELTA[N + o_delta]) * W - DELTA[N - 1 + o_delta] * DELTA[N + o_delta] * (DPSI + DPHI);
437: B = DELTA[N - 1 + o_delta] * DELTA[N + o_delta] * W;
438: if (C < ZERO) C = Math.Abs(C);
439: if (C == ZERO)
440: {
441: // * ETA = B/A
442: // * ETA = RHO - TAU
443: ETA = DLTUB - TAU;
444: }
445: else
446: {
447: if (A >= ZERO)
448: {
449: ETA = (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);
450: }
451: else
452: {
453: ETA = TWO * B / (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));
454: }
455: }
456: // *
457: // * Note, eta should be positive if w is negative, and
458: // * eta should be negative otherwise. However,
459: // * if for some reason caused by roundoff, eta*w > 0,
460: // * we simply use one Newton step instead. This way
461: // * will guarantee eta*w < 0.
462: // *
463: if (W * ETA > ZERO) ETA = - W / (DPSI + DPHI);
464: TEMP = TAU + ETA;
465: if (TEMP > DLTUB || TEMP < DLTLB)
466: {
467: if (W < ZERO)
468: {
469: ETA = (DLTUB - TAU) / TWO;
470: }
471: else
472: {
473: ETA = (DLTLB - TAU) / TWO;
474: }
475: }
476: for (J = 1; J <= N; J++)
477: {
478: DELTA[J + o_delta] = DELTA[J + o_delta] - ETA;
479: }
480: // *
481: TAU = TAU + ETA;
482: // *
483: // * Evaluate PSI and the derivative DPSI
484: // *
485: DPSI = ZERO;
486: PSI = ZERO;
487: ERRETM = ZERO;
488: for (J = 1; J <= II; J++)
489: {
490: TEMP = Z[J + o_z] / DELTA[J + o_delta];
491: PSI = PSI + Z[J + o_z] * TEMP;
492: DPSI = DPSI + TEMP * TEMP;
493: ERRETM = ERRETM + PSI;
494: }
495: ERRETM = Math.Abs(ERRETM);
496: // *
497: // * Evaluate PHI and the derivative DPHI
498: // *
499: TEMP = Z[N + o_z] / DELTA[N + o_delta];
500: PHI = Z[N + o_z] * TEMP;
501: DPHI = TEMP * TEMP;
502: ERRETM = EIGHT * ( - PHI - PSI) + ERRETM - PHI + RHOINV + Math.Abs(TAU) * (DPSI + DPHI);
503: // *
504: W = RHOINV + PHI + PSI;
505: // *
506: // * Main loop to update the values of the array DELTA
507: // *
508: ITER = NITER + 1;
509: // *
510: for (NITER = ITER; NITER <= MAXIT; NITER++)
511: {
512: // *
513: // * Test for convergence
514: // *
515: if (Math.Abs(W) <= EPS * ERRETM)
516: {
517: DLAM = D[I + o_d] + TAU;
518: goto LABEL250;
519: }
520: // *
521: if (W <= ZERO)
522: {
523: DLTLB = Math.Max(DLTLB, TAU);
524: }
525: else
526: {
527: DLTUB = Math.Min(DLTUB, TAU);
528: }
529: // *
530: // * Calculate the new step
531: // *
532: C = W - DELTA[N - 1 + o_delta] * DPSI - DELTA[N + o_delta] * DPHI;
533: A = (DELTA[N - 1 + o_delta] + DELTA[N + o_delta]) * W - DELTA[N - 1 + o_delta] * DELTA[N + o_delta] * (DPSI + DPHI);
534: B = DELTA[N - 1 + o_delta] * DELTA[N + o_delta] * W;
535: if (A >= ZERO)
536: {
537: ETA = (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);
538: }
539: else
540: {
541: ETA = TWO * B / (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));
542: }
543: // *
544: // * Note, eta should be positive if w is negative, and
545: // * eta should be negative otherwise. However,
546: // * if for some reason caused by roundoff, eta*w > 0,
547: // * we simply use one Newton step instead. This way
548: // * will guarantee eta*w < 0.
549: // *
550: if (W * ETA > ZERO) ETA = - W / (DPSI + DPHI);
551: TEMP = TAU + ETA;
552: if (TEMP > DLTUB || TEMP < DLTLB)
553: {
554: if (W < ZERO)
555: {
556: ETA = (DLTUB - TAU) / TWO;
557: }
558: else
559: {
560: ETA = (DLTLB - TAU) / TWO;
561: }
562: }
563: for (J = 1; J <= N; J++)
564: {
565: DELTA[J + o_delta] = DELTA[J + o_delta] - ETA;
566: }
567: // *
568: TAU = TAU + ETA;
569: // *
570: // * Evaluate PSI and the derivative DPSI
571: // *
572: DPSI = ZERO;
573: PSI = ZERO;
574: ERRETM = ZERO;
575: for (J = 1; J <= II; J++)
576: {
577: TEMP = Z[J + o_z] / DELTA[J + o_delta];
578: PSI = PSI + Z[J + o_z] * TEMP;
579: DPSI = DPSI + TEMP * TEMP;
580: ERRETM = ERRETM + PSI;
581: }
582: ERRETM = Math.Abs(ERRETM);
583: // *
584: // * Evaluate PHI and the derivative DPHI
585: // *
586: TEMP = Z[N + o_z] / DELTA[N + o_delta];
587: PHI = Z[N + o_z] * TEMP;
588: DPHI = TEMP * TEMP;
589: ERRETM = EIGHT * ( - PHI - PSI) + ERRETM - PHI + RHOINV + Math.Abs(TAU) * (DPSI + DPHI);
590: // *
591: W = RHOINV + PHI + PSI;
592: }
593: // *
594: // * Return with INFO = 1, NITER = MAXIT and not converged
595: // *
596: INFO = 1;
597: DLAM = D[I + o_d] + TAU;
598: goto LABEL250;
599: // *
600: // * End for the case I = N
601: // *
602: }
603: else
604: {
605: // *
606: // * The case for I < N
607: // *
608: NITER = 1;
609: IP1 = I + 1;
610: // *
611: // * Calculate initial guess
612: // *
613: DEL = D[IP1 + o_d] - D[I + o_d];
614: MIDPT = DEL / TWO;
615: for (J = 1; J <= N; J++)
616: {
617: DELTA[J + o_delta] = (D[J + o_d] - D[I + o_d]) - MIDPT;
618: }
619: // *
620: PSI = ZERO;
621: for (J = 1; J <= I - 1; J++)
622: {
623: PSI = PSI + Z[J + o_z] * Z[J + o_z] / DELTA[J + o_delta];
624: }
625: // *
626: PHI = ZERO;
627: for (J = N; J >= I + 2; J += - 1)
628: {
629: PHI = PHI + Z[J + o_z] * Z[J + o_z] / DELTA[J + o_delta];
630: }
631: C = RHOINV + PSI + PHI;
632: W = C + Z[I + o_z] * Z[I + o_z] / DELTA[I + o_delta] + Z[IP1 + o_z] * Z[IP1 + o_z] / DELTA[IP1 + o_delta];
633: // *
634: if (W > ZERO)
635: {
636: // *
637: // * d(i)< the ith eigenvalue < (d(i)+d(i+1))/2
638: // *
639: // * We choose d(i) as origin.
640: // *
641: ORGATI = true;
642: A = C * DEL + Z[I + o_z] * Z[I + o_z] + Z[IP1 + o_z] * Z[IP1 + o_z];
643: B = Z[I + o_z] * Z[I + o_z] * DEL;
644: if (A > ZERO)
645: {
646: TAU = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));
647: }
648: else
649: {
650: TAU = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);
651: }
652: DLTLB = ZERO;
653: DLTUB = MIDPT;
654: }
655: else
656: {
657: // *
658: // * (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1)
659: // *
660: // * We choose d(i+1) as origin.
661: // *
662: ORGATI = false;
663: A = C * DEL - Z[I + o_z] * Z[I + o_z] - Z[IP1 + o_z] * Z[IP1 + o_z];
664: B = Z[IP1 + o_z] * Z[IP1 + o_z] * DEL;
665: if (A < ZERO)
666: {
667: TAU = TWO * B / (A - Math.Sqrt(Math.Abs(A * A + FOUR * B * C)));
668: }
669: else
670: {
671: TAU = - (A + Math.Sqrt(Math.Abs(A * A + FOUR * B * C))) / (TWO * C);
672: }
673: DLTLB = - MIDPT;
674: DLTUB = ZERO;
675: }
676: // *
677: if (ORGATI)
678: {
679: for (J = 1; J <= N; J++)
680: {
681: DELTA[J + o_delta] = (D[J + o_d] - D[I + o_d]) - TAU;
682: }
683: }
684: else
685: {
686: for (J = 1; J <= N; J++)
687: {
688: DELTA[J + o_delta] = (D[J + o_d] - D[IP1 + o_d]) - TAU;
689: }
690: }
691: if (ORGATI)
692: {
693: II = I;
694: }
695: else
696: {
697: II = I + 1;
698: }
699: IIM1 = II - 1;
700: IIP1 = II + 1;
701: // *
702: // * Evaluate PSI and the derivative DPSI
703: // *
704: DPSI = ZERO;
705: PSI = ZERO;
706: ERRETM = ZERO;
707: for (J = 1; J <= IIM1; J++)
708: {
709: TEMP = Z[J + o_z] / DELTA[J + o_delta];
710: PSI = PSI + Z[J + o_z] * TEMP;
711: DPSI = DPSI + TEMP * TEMP;
712: ERRETM = ERRETM + PSI;
713: }
714: ERRETM = Math.Abs(ERRETM);
715: // *
716: // * Evaluate PHI and the derivative DPHI
717: // *
718: DPHI = ZERO;
719: PHI = ZERO;
720: for (J = N; J >= IIP1; J += - 1)
721: {
722: TEMP = Z[J + o_z] / DELTA[J + o_delta];
723: PHI = PHI + Z[J + o_z] * TEMP;
724: DPHI = DPHI + TEMP * TEMP;
725: ERRETM = ERRETM + PHI;
726: }
727: // *
728: W = RHOINV + PHI + PSI;
729: // *
730: // * W is the value of the secular function with
731: // * its ii-th element removed.
732: // *
733: SWTCH3 = false;
734: if (ORGATI)
735: {
736: if (W < ZERO) SWTCH3 = true;
737: }
738: else
739: {
740: if (W > ZERO) SWTCH3 = true;
741: }
742: if (II == 1 || II == N) SWTCH3 = false;
743: // *
744: TEMP = Z[II + o_z] / DELTA[II + o_delta];
745: DW = DPSI + DPHI + TEMP * TEMP;
746: TEMP = Z[II + o_z] * TEMP;
747: W = W + TEMP;
748: ERRETM = EIGHT * (PHI - PSI) + ERRETM + TWO * RHOINV + THREE * Math.Abs(TEMP) + Math.Abs(TAU) * DW;
749: // *
750: // * Test for convergence
751: // *
752: if (Math.Abs(W) <= EPS * ERRETM)
753: {
754: if (ORGATI)
755: {
756: DLAM = D[I + o_d] + TAU;
757: }
758: else
759: {
760: DLAM = D[IP1 + o_d] + TAU;
761: }
762: goto LABEL250;
763: }
764: // *
765: if (W <= ZERO)
766: {
767: DLTLB = Math.Max(DLTLB, TAU);
768: }
769: else
770: {
771: DLTUB = Math.Min(DLTUB, TAU);
772: }
773: // *
774: // * Calculate the new step
775: // *
776: NITER = NITER + 1;
777: if (!SWTCH3)
778: {
779: if (ORGATI)
780: {
781: C = W - DELTA[IP1 + o_delta] * DW - (D[I + o_d] - D[IP1 + o_d]) * Math.Pow(Z[I + o_z] / DELTA[I + o_delta],2);
782: }
783: else
784: {
785: C = W - DELTA[I + o_delta] * DW - (D[IP1 + o_d] - D[I + o_d]) * Math.Pow(Z[IP1 + o_z] / DELTA[IP1 + o_delta],2);
786: }
787: A = (DELTA[I + o_delta] + DELTA[IP1 + o_delta]) * W - DELTA[I + o_delta] * DELTA[IP1 + o_delta] * DW;
788: B = DELTA[I + o_delta] * DELTA[IP1 + o_delta] * W;
789: if (C == ZERO)
790: {
791: if (A == ZERO)
792: {
793: if (ORGATI)
794: {
795: A = Z[I + o_z] * Z[I + o_z] + DELTA[IP1 + o_delta] * DELTA[IP1 + o_delta] * (DPSI + DPHI);
796: }
797: else
798: {
799: A = Z[IP1 + o_z] * Z[IP1 + o_z] + DELTA[I + o_delta] * DELTA[I + o_delta] * (DPSI + DPHI);
800: }
801: }
802: ETA = B / A;
803: }
804: else
805: {
806: if (A <= ZERO)
807: {
808: ETA = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);
809: }
810: else
811: {
812: ETA = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));
813: }
814: }
815: }
816: else
817: {
818: // *
819: // * Interpolation using THREE most relevant poles
820: // *
821: TEMP = RHOINV + PSI + PHI;
822: if (ORGATI)
823: {
824: TEMP1 = Z[IIM1 + o_z] / DELTA[IIM1 + o_delta];
825: TEMP1 = TEMP1 * TEMP1;
826: C = TEMP - DELTA[IIP1 + o_delta] * (DPSI + DPHI) - (D[IIM1 + o_d] - D[IIP1 + o_d]) * TEMP1;
827: ZZ[1 + o_zz] = Z[IIM1 + o_z] * Z[IIM1 + o_z];
828: ZZ[3 + o_zz] = DELTA[IIP1 + o_delta] * DELTA[IIP1 + o_delta] * ((DPSI - TEMP1) + DPHI);
829: }
830: else
831: {
832: TEMP1 = Z[IIP1 + o_z] / DELTA[IIP1 + o_delta];
833: TEMP1 = TEMP1 * TEMP1;
834: C = TEMP - DELTA[IIM1 + o_delta] * (DPSI + DPHI) - (D[IIP1 + o_d] - D[IIM1 + o_d]) * TEMP1;
835: ZZ[1 + o_zz] = DELTA[IIM1 + o_delta] * DELTA[IIM1 + o_delta] * (DPSI + (DPHI - TEMP1));
836: ZZ[3 + o_zz] = Z[IIP1 + o_z] * Z[IIP1 + o_z];
837: }
838: ZZ[2 + o_zz] = Z[II + o_z] * Z[II + o_z];
839: this._dlaed6.Run(NITER, ORGATI, C, DELTA, IIM1 + o_delta, ZZ, offset_zz, W
840: , ref ETA, ref INFO);
841: if (INFO != 0) goto LABEL250;
842: }
843: // *
844: // * Note, eta should be positive if w is negative, and
845: // * eta should be negative otherwise. However,
846: // * if for some reason caused by roundoff, eta*w > 0,
847: // * we simply use one Newton step instead. This way
848: // * will guarantee eta*w < 0.
849: // *
850: if (W * ETA >= ZERO) ETA = - W / DW;
851: TEMP = TAU + ETA;
852: if (TEMP > DLTUB || TEMP < DLTLB)
853: {
854: if (W < ZERO)
855: {
856: ETA = (DLTUB - TAU) / TWO;
857: }
858: else
859: {
860: ETA = (DLTLB - TAU) / TWO;
861: }
862: }
863: // *
864: PREW = W;
865: // *
866: for (J = 1; J <= N; J++)
867: {
868: DELTA[J + o_delta] = DELTA[J + o_delta] - ETA;
869: }
870: // *
871: // * Evaluate PSI and the derivative DPSI
872: // *
873: DPSI = ZERO;
874: PSI = ZERO;
875: ERRETM = ZERO;
876: for (J = 1; J <= IIM1; J++)
877: {
878: TEMP = Z[J + o_z] / DELTA[J + o_delta];
879: PSI = PSI + Z[J + o_z] * TEMP;
880: DPSI = DPSI + TEMP * TEMP;
881: ERRETM = ERRETM + PSI;
882: }
883: ERRETM = Math.Abs(ERRETM);
884: // *
885: // * Evaluate PHI and the derivative DPHI
886: // *
887: DPHI = ZERO;
888: PHI = ZERO;
889: for (J = N; J >= IIP1; J += - 1)
890: {
891: TEMP = Z[J + o_z] / DELTA[J + o_delta];
892: PHI = PHI + Z[J + o_z] * TEMP;
893: DPHI = DPHI + TEMP * TEMP;
894: ERRETM = ERRETM + PHI;
895: }
896: // *
897: TEMP = Z[II + o_z] / DELTA[II + o_delta];
898: DW = DPSI + DPHI + TEMP * TEMP;
899: TEMP = Z[II + o_z] * TEMP;
900: W = RHOINV + PHI + PSI + TEMP;
901: ERRETM = EIGHT * (PHI - PSI) + ERRETM + TWO * RHOINV + THREE * Math.Abs(TEMP) + Math.Abs(TAU + ETA) * DW;
902: // *
903: SWTCH = false;
904: if (ORGATI)
905: {
906: if ( - W > Math.Abs(PREW) / TEN) SWTCH = true;
907: }
908: else
909: {
910: if (W > Math.Abs(PREW) / TEN) SWTCH = true;
911: }
912: // *
913: TAU = TAU + ETA;
914: // *
915: // * Main loop to update the values of the array DELTA
916: // *
917: ITER = NITER + 1;
918: // *
919: for (NITER = ITER; NITER <= MAXIT; NITER++)
920: {
921: // *
922: // * Test for convergence
923: // *
924: if (Math.Abs(W) <= EPS * ERRETM)
925: {
926: if (ORGATI)
927: {
928: DLAM = D[I + o_d] + TAU;
929: }
930: else
931: {
932: DLAM = D[IP1 + o_d] + TAU;
933: }
934: goto LABEL250;
935: }
936: // *
937: if (W <= ZERO)
938: {
939: DLTLB = Math.Max(DLTLB, TAU);
940: }
941: else
942: {
943: DLTUB = Math.Min(DLTUB, TAU);
944: }
945: // *
946: // * Calculate the new step
947: // *
948: if (!SWTCH3)
949: {
950: if (!SWTCH)
951: {
952: if (ORGATI)
953: {
954: C = W - DELTA[IP1 + o_delta] * DW - (D[I + o_d] - D[IP1 + o_d]) * Math.Pow(Z[I + o_z] / DELTA[I + o_delta],2);
955: }
956: else
957: {
958: C = W - DELTA[I + o_delta] * DW - (D[IP1 + o_d] - D[I + o_d]) * Math.Pow(Z[IP1 + o_z] / DELTA[IP1 + o_delta],2);
959: }
960: }
961: else
962: {
963: TEMP = Z[II + o_z] / DELTA[II + o_delta];
964: if (ORGATI)
965: {
966: DPSI = DPSI + TEMP * TEMP;
967: }
968: else
969: {
970: DPHI = DPHI + TEMP * TEMP;
971: }
972: C = W - DELTA[I + o_delta] * DPSI - DELTA[IP1 + o_delta] * DPHI;
973: }
974: A = (DELTA[I + o_delta] + DELTA[IP1 + o_delta]) * W - DELTA[I + o_delta] * DELTA[IP1 + o_delta] * DW;
975: B = DELTA[I + o_delta] * DELTA[IP1 + o_delta] * W;
976: if (C == ZERO)
977: {
978: if (A == ZERO)
979: {
980: if (!SWTCH)
981: {
982: if (ORGATI)
983: {
984: A = Z[I + o_z] * Z[I + o_z] + DELTA[IP1 + o_delta] * DELTA[IP1 + o_delta] * (DPSI + DPHI);
985: }
986: else
987: {
988: A = Z[IP1 + o_z] * Z[IP1 + o_z] + DELTA[I + o_delta] * DELTA[I + o_delta] * (DPSI + DPHI);
989: }
990: }
991: else
992: {
993: A = DELTA[I + o_delta] * DELTA[I + o_delta] * DPSI + DELTA[IP1 + o_delta] * DELTA[IP1 + o_delta] * DPHI;
994: }
995: }
996: ETA = B / A;
997: }
998: else
999: {
1000: if (A <= ZERO)
1001: {
1002: ETA = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);
1003: }
1004: else
1005: {
1006: ETA = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));
1007: }
1008: }
1009: }
1010: else
1011: {
1012: // *
1013: // * Interpolation using THREE most relevant poles
1014: // *
1015: TEMP = RHOINV + PSI + PHI;
1016: if (SWTCH)
1017: {
1018: C = TEMP - DELTA[IIM1 + o_delta] * DPSI - DELTA[IIP1 + o_delta] * DPHI;
1019: ZZ[1 + o_zz] = DELTA[IIM1 + o_delta] * DELTA[IIM1 + o_delta] * DPSI;
1020: ZZ[3 + o_zz] = DELTA[IIP1 + o_delta] * DELTA[IIP1 + o_delta] * DPHI;
1021: }
1022: else
1023: {
1024: if (ORGATI)
1025: {
1026: TEMP1 = Z[IIM1 + o_z] / DELTA[IIM1 + o_delta];
1027: TEMP1 = TEMP1 * TEMP1;
1028: C = TEMP - DELTA[IIP1 + o_delta] * (DPSI + DPHI) - (D[IIM1 + o_d] - D[IIP1 + o_d]) * TEMP1;
1029: ZZ[1 + o_zz] = Z[IIM1 + o_z] * Z[IIM1 + o_z];
1030: ZZ[3 + o_zz] = DELTA[IIP1 + o_delta] * DELTA[IIP1 + o_delta] * ((DPSI - TEMP1) + DPHI);
1031: }
1032: else
1033: {
1034: TEMP1 = Z[IIP1 + o_z] / DELTA[IIP1 + o_delta];
1035: TEMP1 = TEMP1 * TEMP1;
1036: C = TEMP - DELTA[IIM1 + o_delta] * (DPSI + DPHI) - (D[IIP1 + o_d] - D[IIM1 + o_d]) * TEMP1;
1037: ZZ[1 + o_zz] = DELTA[IIM1 + o_delta] * DELTA[IIM1 + o_delta] * (DPSI + (DPHI - TEMP1));
1038: ZZ[3 + o_zz] = Z[IIP1 + o_z] * Z[IIP1 + o_z];
1039: }
1040: }
1041: this._dlaed6.Run(NITER, ORGATI, C, DELTA, IIM1 + o_delta, ZZ, offset_zz, W
1042: , ref ETA, ref INFO);
1043: if (INFO != 0) goto LABEL250;
1044: }
1045: // *
1046: // * Note, eta should be positive if w is negative, and
1047: // * eta should be negative otherwise. However,
1048: // * if for some reason caused by roundoff, eta*w > 0,
1049: // * we simply use one Newton step instead. This way
1050: // * will guarantee eta*w < 0.
1051: // *
1052: if (W * ETA >= ZERO) ETA = - W / DW;
1053: TEMP = TAU + ETA;
1054: if (TEMP > DLTUB || TEMP < DLTLB)
1055: {
1056: if (W < ZERO)
1057: {
1058: ETA = (DLTUB - TAU) / TWO;
1059: }
1060: else
1061: {
1062: ETA = (DLTLB - TAU) / TWO;
1063: }
1064: }
1065: // *
1066: for (J = 1; J <= N; J++)
1067: {
1068: DELTA[J + o_delta] = DELTA[J + o_delta] - ETA;
1069: }
1070: // *
1071: TAU = TAU + ETA;
1072: PREW = W;
1073: // *
1074: // * Evaluate PSI and the derivative DPSI
1075: // *
1076: DPSI = ZERO;
1077: PSI = ZERO;
1078: ERRETM = ZERO;
1079: for (J = 1; J <= IIM1; J++)
1080: {
1081: TEMP = Z[J + o_z] / DELTA[J + o_delta];
1082: PSI = PSI + Z[J + o_z] * TEMP;
1083: DPSI = DPSI + TEMP * TEMP;
1084: ERRETM = ERRETM + PSI;
1085: }
1086: ERRETM = Math.Abs(ERRETM);
1087: // *
1088: // * Evaluate PHI and the derivative DPHI
1089: // *
1090: DPHI = ZERO;
1091: PHI = ZERO;
1092: for (J = N; J >= IIP1; J += - 1)
1093: {
1094: TEMP = Z[J + o_z] / DELTA[J + o_delta];
1095: PHI = PHI + Z[J + o_z] * TEMP;
1096: DPHI = DPHI + TEMP * TEMP;
1097: ERRETM = ERRETM + PHI;
1098: }
1099: // *
1100: TEMP = Z[II + o_z] / DELTA[II + o_delta];
1101: DW = DPSI + DPHI + TEMP * TEMP;
1102: TEMP = Z[II + o_z] * TEMP;
1103: W = RHOINV + PHI + PSI + TEMP;
1104: ERRETM = EIGHT * (PHI - PSI) + ERRETM + TWO * RHOINV + THREE * Math.Abs(TEMP) + Math.Abs(TAU) * DW;
1105: if (W * PREW > ZERO && Math.Abs(W) > Math.Abs(PREW) / TEN) SWTCH = !SWTCH;
1106: // *
1107: }
1108: // *
1109: // * Return with INFO = 1, NITER = MAXIT and not converged
1110: // *
1111: INFO = 1;
1112: if (ORGATI)
1113: {
1114: DLAM = D[I + o_d] + TAU;
1115: }
1116: else
1117: {
1118: DLAM = D[IP1 + o_d] + TAU;
1119: }
1120: // *
1121: }
1122: // *
1123: LABEL250:;
1124: // *
1125: return;
1126: // *
1127: // * End of DLAED4
1128: // *
1129:
1130: #endregion
1131:
1132: }
1133: }
1134: }