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.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// February 2007
23: /// Purpose
24: /// =======
25: ///
26: /// DLAED6 computes the positive or negative root (closest to the origin)
27: /// of
28: /// z(1) z(2) z(3)
29: /// f(x) = rho + --------- + ---------- + ---------
30: /// d(1)-x d(2)-x d(3)-x
31: ///
32: /// It is assumed that
33: ///
34: /// if ORGATI = .true. the root is between d(2) and d(3);
35: /// otherwise it is between d(1) and d(2)
36: ///
37: /// This routine will be called by DLAED4 when necessary. In most cases,
38: /// the root sought is the smallest in magnitude, though it might not be
39: /// in some extremely rare situations.
40: ///
41: ///</summary>
42: public class DLAED6
43: {
44:
45:
46: #region Dependencies
47:
48: DLAMCH _dlamch;
49:
50: #endregion
51:
52:
53: #region Fields
54:
55: const int MAXIT = 40; const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0;
56: const double THREE = 3.0E0;const double FOUR = 4.0E0; const double EIGHT = 8.0E0;
57: double[] DSCALE = new double[3]; int o_dscale = -1;
58: double[] ZSCALE = new double[3]; int o_zscale = -1;bool SCALE = false; int I = 0; int ITER = 0;
59: int NITER = 0;double A = 0; double B = 0; double BASE = 0; double C = 0; double DDF = 0; double DF = 0; double EPS = 0;
60: double ERRETM = 0;double ETA = 0; double F = 0; double FC = 0; double SCLFAC = 0; double SCLINV = 0; double SMALL1 = 0;
61: double SMALL2 = 0;double SMINV1 = 0; double SMINV2 = 0; double TEMP = 0; double TEMP1 = 0; double TEMP2 = 0;
62: double TEMP3 = 0;double TEMP4 = 0; double LBD = 0; double UBD = 0;
63:
64: #endregion
65:
66: public DLAED6(DLAMCH dlamch)
67: {
68:
69:
70: #region Set Dependencies
71:
72: this._dlamch = dlamch;
73:
74: #endregion
75:
76: }
77:
78: public DLAED6()
79: {
80:
81:
82: #region Dependencies (Initialization)
83:
84: LSAME lsame = new LSAME();
85: DLAMC3 dlamc3 = new DLAMC3();
86: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
87: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
88: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
89: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
90: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
91:
92: #endregion
93:
94:
95: #region Set Dependencies
96:
97: this._dlamch = dlamch;
98:
99: #endregion
100:
101: }
102: /// <summary>
103: /// Purpose
104: /// =======
105: ///
106: /// DLAED6 computes the positive or negative root (closest to the origin)
107: /// of
108: /// z(1) z(2) z(3)
109: /// f(x) = rho + --------- + ---------- + ---------
110: /// d(1)-x d(2)-x d(3)-x
111: ///
112: /// It is assumed that
113: ///
114: /// if ORGATI = .true. the root is between d(2) and d(3);
115: /// otherwise it is between d(1) and d(2)
116: ///
117: /// This routine will be called by DLAED4 when necessary. In most cases,
118: /// the root sought is the smallest in magnitude, though it might not be
119: /// in some extremely rare situations.
120: ///
121: ///</summary>
122: /// <param name="KNITER">
123: /// (input) INTEGER
124: /// Refer to DLAED4 for its significance.
125: ///</param>
126: /// <param name="ORGATI">
127: /// (input) LOGICAL
128: /// If ORGATI is true, the needed root is between d(2) and
129: /// d(3); otherwise it is between d(1) and d(2). See
130: /// DLAED4 for further details.
131: ///</param>
132: /// <param name="RHO">
133: /// (input) DOUBLE PRECISION
134: /// Refer to the equation f(x) above.
135: ///</param>
136: /// <param name="D">
137: /// (input) DOUBLE PRECISION array, dimension (3)
138: /// D satisfies d(1) .LT. d(2) .LT. d(3).
139: ///</param>
140: /// <param name="Z">
141: /// (input) DOUBLE PRECISION array, dimension (3)
142: /// Each of the elements in z must be positive.
143: ///</param>
144: /// <param name="FINIT">
145: /// (input) DOUBLE PRECISION
146: /// The value of f at 0. It is more accurate than the one
147: /// evaluated inside this routine (if someone wants to do
148: /// so).
149: ///</param>
150: /// <param name="TAU">
151: /// (output) DOUBLE PRECISION
152: /// The root of the equation f(x).
153: ///</param>
154: /// <param name="INFO">
155: /// (output) INTEGER
156: /// = 0: successful exit
157: /// .GT. 0: if INFO = 1, failure to converge
158: ///</param>
159: public void Run(int KNITER, bool ORGATI, double RHO, double[] D, int offset_d, double[] Z, int offset_z, double FINIT
160: , ref double TAU, ref int INFO)
161: {
162:
163: #region Array Index Correction
164:
165: int o_d = -1 + offset_d; int o_z = -1 + offset_z;
166:
167: #endregion
168:
169:
170: #region Prolog
171:
172: // *
173: // * -- LAPACK routine (version 3.1.1) --
174: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
175: // * February 2007
176: // *
177: // * .. Scalar Arguments ..
178: // * ..
179: // * .. Array Arguments ..
180: // * ..
181: // *
182: // * Purpose
183: // * =======
184: // *
185: // * DLAED6 computes the positive or negative root (closest to the origin)
186: // * of
187: // * z(1) z(2) z(3)
188: // * f(x) = rho + --------- + ---------- + ---------
189: // * d(1)-x d(2)-x d(3)-x
190: // *
191: // * It is assumed that
192: // *
193: // * if ORGATI = .true. the root is between d(2) and d(3);
194: // * otherwise it is between d(1) and d(2)
195: // *
196: // * This routine will be called by DLAED4 when necessary. In most cases,
197: // * the root sought is the smallest in magnitude, though it might not be
198: // * in some extremely rare situations.
199: // *
200: // * Arguments
201: // * =========
202: // *
203: // * KNITER (input) INTEGER
204: // * Refer to DLAED4 for its significance.
205: // *
206: // * ORGATI (input) LOGICAL
207: // * If ORGATI is true, the needed root is between d(2) and
208: // * d(3); otherwise it is between d(1) and d(2). See
209: // * DLAED4 for further details.
210: // *
211: // * RHO (input) DOUBLE PRECISION
212: // * Refer to the equation f(x) above.
213: // *
214: // * D (input) DOUBLE PRECISION array, dimension (3)
215: // * D satisfies d(1) < d(2) < d(3).
216: // *
217: // * Z (input) DOUBLE PRECISION array, dimension (3)
218: // * Each of the elements in z must be positive.
219: // *
220: // * FINIT (input) DOUBLE PRECISION
221: // * The value of f at 0. It is more accurate than the one
222: // * evaluated inside this routine (if someone wants to do
223: // * so).
224: // *
225: // * TAU (output) DOUBLE PRECISION
226: // * The root of the equation f(x).
227: // *
228: // * INFO (output) INTEGER
229: // * = 0: successful exit
230: // * > 0: if INFO = 1, failure to converge
231: // *
232: // * Further Details
233: // * ===============
234: // *
235: // * 30/06/99: Based on contributions by
236: // * Ren-Cang Li, Computer Science Division, University of California
237: // * at Berkeley, USA
238: // *
239: // * 10/02/03: This version has a few statements commented out for thread
240: // * safety (machine parameters are computed on each entry). SJH.
241: // *
242: // * 05/10/06: Modified from a new version of Ren-Cang Li, use
243: // * Gragg-Thornton-Warner cubic convergent scheme for better stability.
244: // *
245: // * =====================================================================
246: // *
247: // * .. Parameters ..
248: // * ..
249: // * .. External Functions ..
250: // * ..
251: // * .. Local Arrays ..
252: // * ..
253: // * .. Local Scalars ..
254: // * ..
255: // * .. Intrinsic Functions ..
256: // INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT;
257: // * ..
258: // * .. Executable Statements ..
259: // *
260:
261: #endregion
262:
263:
264: #region Body
265:
266: INFO = 0;
267: // *
268: if (ORGATI)
269: {
270: LBD = D[2 + o_d];
271: UBD = D[3 + o_d];
272: }
273: else
274: {
275: LBD = D[1 + o_d];
276: UBD = D[2 + o_d];
277: }
278: if (FINIT < ZERO)
279: {
280: LBD = ZERO;
281: }
282: else
283: {
284: UBD = ZERO;
285: }
286: // *
287: NITER = 1;
288: TAU = ZERO;
289: if (KNITER == 2)
290: {
291: if (ORGATI)
292: {
293: TEMP = (D[3 + o_d] - D[2 + o_d]) / TWO;
294: C = RHO + Z[1 + o_z] / ((D[1 + o_d] - D[2 + o_d]) - TEMP);
295: A = C * (D[2 + o_d] + D[3 + o_d]) + Z[2 + o_z] + Z[3 + o_z];
296: B = C * D[2 + o_d] * D[3 + o_d] + Z[2 + o_z] * D[3 + o_d] + Z[3 + o_z] * D[2 + o_d];
297: }
298: else
299: {
300: TEMP = (D[1 + o_d] - D[2 + o_d]) / TWO;
301: C = RHO + Z[3 + o_z] / ((D[3 + o_d] - D[2 + o_d]) - TEMP);
302: A = C * (D[1 + o_d] + D[2 + o_d]) + Z[1 + o_z] + Z[2 + o_z];
303: B = C * D[1 + o_d] * D[2 + o_d] + Z[1 + o_z] * D[2 + o_d] + Z[2 + o_z] * D[1 + o_d];
304: }
305: TEMP = Math.Max(Math.Abs(A), Math.Max(Math.Abs(B), Math.Abs(C)));
306: A = A / TEMP;
307: B = B / TEMP;
308: C = C / TEMP;
309: if (C == ZERO)
310: {
311: TAU = B / A;
312: }
313: else
314: {
315: if (A <= ZERO)
316: {
317: TAU = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);
318: }
319: else
320: {
321: TAU = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));
322: }
323: }
324: if (TAU < LBD || TAU > UBD) TAU = (LBD + UBD) / TWO;
325: if (D[1 + o_d] == TAU || D[2 + o_d] == TAU || D[3 + o_d] == TAU)
326: {
327: TAU = ZERO;
328: }
329: else
330: {
331: TEMP = FINIT + TAU * Z[1 + o_z] / (D[1 + o_d] * (D[1 + o_d] - TAU)) + TAU * Z[2 + o_z] / (D[2 + o_d] * (D[2 + o_d] - TAU)) + TAU * Z[3 + o_z] / (D[3 + o_d] * (D[3 + o_d] - TAU));
332: if (TEMP <= ZERO)
333: {
334: LBD = TAU;
335: }
336: else
337: {
338: UBD = TAU;
339: }
340: if (Math.Abs(FINIT) <= Math.Abs(TEMP)) TAU = ZERO;
341: }
342: }
343: // *
344: // * get machine parameters for possible scaling to avoid overflow
345: // *
346: // * modified by Sven: parameters SMALL1, SMINV1, SMALL2,
347: // * SMINV2, EPS are not SAVEd anymore between one call to the
348: // * others but recomputed at each call
349: // *
350: EPS = this._dlamch.Run("Epsilon");
351: BASE = this._dlamch.Run("Base");
352: SMALL1 = Math.Pow(BASE,Convert.ToInt32(Math.Truncate(Math.Log(this._dlamch.Run("SafMin")) / Math.Log(BASE) / THREE)));
353: SMINV1 = ONE / SMALL1;
354: SMALL2 = SMALL1 * SMALL1;
355: SMINV2 = SMINV1 * SMINV1;
356: // *
357: // * Determine if scaling of inputs necessary to avoid overflow
358: // * when computing 1/TEMP**3
359: // *
360: if (ORGATI)
361: {
362: TEMP = Math.Min(Math.Abs(D[2 + o_d] - TAU), Math.Abs(D[3 + o_d] - TAU));
363: }
364: else
365: {
366: TEMP = Math.Min(Math.Abs(D[1 + o_d] - TAU), Math.Abs(D[2 + o_d] - TAU));
367: }
368: SCALE = false;
369: if (TEMP <= SMALL1)
370: {
371: SCALE = true;
372: if (TEMP <= SMALL2)
373: {
374: // *
375: // * Scale up by power of radix nearest 1/SAFMIN**(2/3)
376: // *
377: SCLFAC = SMINV2;
378: SCLINV = SMALL2;
379: }
380: else
381: {
382: // *
383: // * Scale up by power of radix nearest 1/SAFMIN**(1/3)
384: // *
385: SCLFAC = SMINV1;
386: SCLINV = SMALL1;
387: }
388: // *
389: // * Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
390: // *
391: for (I = 1; I <= 3; I++)
392: {
393: DSCALE[I + o_dscale] = D[I + o_d] * SCLFAC;
394: ZSCALE[I + o_zscale] = Z[I + o_z] * SCLFAC;
395: }
396: TAU = TAU * SCLFAC;
397: LBD = LBD * SCLFAC;
398: UBD = UBD * SCLFAC;
399: }
400: else
401: {
402: // *
403: // * Copy D and Z to DSCALE and ZSCALE
404: // *
405: for (I = 1; I <= 3; I++)
406: {
407: DSCALE[I + o_dscale] = D[I + o_d];
408: ZSCALE[I + o_zscale] = Z[I + o_z];
409: }
410: }
411: // *
412: FC = ZERO;
413: DF = ZERO;
414: DDF = ZERO;
415: for (I = 1; I <= 3; I++)
416: {
417: TEMP = ONE / (DSCALE[I + o_dscale] - TAU);
418: TEMP1 = ZSCALE[I + o_zscale] * TEMP;
419: TEMP2 = TEMP1 * TEMP;
420: TEMP3 = TEMP2 * TEMP;
421: FC = FC + TEMP1 / DSCALE[I + o_dscale];
422: DF = DF + TEMP2;
423: DDF = DDF + TEMP3;
424: }
425: F = FINIT + TAU * FC;
426: // *
427: if (Math.Abs(F) <= ZERO) goto LABEL60;
428: if (F <= ZERO)
429: {
430: LBD = TAU;
431: }
432: else
433: {
434: UBD = TAU;
435: }
436: // *
437: // * Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
438: // * scheme
439: // *
440: // * It is not hard to see that
441: // *
442: // * 1) Iterations will go up monotonically
443: // * if FINIT < 0;
444: // *
445: // * 2) Iterations will go down monotonically
446: // * if FINIT > 0.
447: // *
448: ITER = NITER + 1;
449: // *
450: for (NITER = ITER; NITER <= MAXIT; NITER++)
451: {
452: // *
453: if (ORGATI)
454: {
455: TEMP1 = DSCALE[2 + o_dscale] - TAU;
456: TEMP2 = DSCALE[3 + o_dscale] - TAU;
457: }
458: else
459: {
460: TEMP1 = DSCALE[1 + o_dscale] - TAU;
461: TEMP2 = DSCALE[2 + o_dscale] - TAU;
462: }
463: A = (TEMP1 + TEMP2) * F - TEMP1 * TEMP2 * DF;
464: B = TEMP1 * TEMP2 * F;
465: C = F - (TEMP1 + TEMP2) * DF + TEMP1 * TEMP2 * DDF;
466: TEMP = Math.Max(Math.Abs(A), Math.Max(Math.Abs(B), Math.Abs(C)));
467: A = A / TEMP;
468: B = B / TEMP;
469: C = C / TEMP;
470: if (C == ZERO)
471: {
472: ETA = B / A;
473: }
474: else
475: {
476: if (A <= ZERO)
477: {
478: ETA = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);
479: }
480: else
481: {
482: ETA = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));
483: }
484: }
485: if (F * ETA >= ZERO)
486: {
487: ETA = - F / DF;
488: }
489: // *
490: TAU = TAU + ETA;
491: if (TAU < LBD || TAU > UBD) TAU = (LBD + UBD) / TWO;
492: // *
493: FC = ZERO;
494: ERRETM = ZERO;
495: DF = ZERO;
496: DDF = ZERO;
497: for (I = 1; I <= 3; I++)
498: {
499: TEMP = ONE / (DSCALE[I + o_dscale] - TAU);
500: TEMP1 = ZSCALE[I + o_zscale] * TEMP;
501: TEMP2 = TEMP1 * TEMP;
502: TEMP3 = TEMP2 * TEMP;
503: TEMP4 = TEMP1 / DSCALE[I + o_dscale];
504: FC = FC + TEMP4;
505: ERRETM = ERRETM + Math.Abs(TEMP4);
506: DF = DF + TEMP2;
507: DDF = DDF + TEMP3;
508: }
509: F = FINIT + TAU * FC;
510: ERRETM = EIGHT * (Math.Abs(FINIT) + Math.Abs(TAU) * ERRETM) + Math.Abs(TAU) * DF;
511: if (Math.Abs(F) <= EPS * ERRETM) goto LABEL60;
512: if (F <= ZERO)
513: {
514: LBD = TAU;
515: }
516: else
517: {
518: UBD = TAU;
519: }
520: }
521: INFO = 1;
522: LABEL60:;
523: // *
524: // * Undo scaling
525: // *
526: if (SCALE) TAU = TAU * SCLINV;
527: return;
528: // *
529: // * End of DLAED6
530: // *
531:
532: #endregion
533:
534: }
535: }
536: }