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