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: /// DLAZQ3 checks for deflation, computes a shift (TAU) and calls dqds.
27: /// In case of failure it changes shifts, and tries again until output
28: /// is positive.
29: ///
30: ///</summary>
31: public class DLAZQ3
32: {
33:
34:
35: #region Dependencies
36:
37: DLASQ5 _dlasq5; DLASQ6 _dlasq6; DLAZQ4 _dlazq4; DLAMCH _dlamch;
38:
39: #endregion
40:
41:
42: #region Fields
43:
44: const double CBIAS = 1.50E0; const double ZERO = 0.0E0; const double QURTR = 0.250E0; const double HALF = 0.5E0;
45: const double ONE = 1.0E0;const double TWO = 2.0E0; const double HUNDRD = 100.0E0; int IPN4 = 0; int J4 = 0; int N0IN = 0;
46: int NN = 0;double EPS = 0; double G = 0; double S = 0; double SAFMIN = 0; double T = 0; double TEMP = 0; double TOL = 0;
47: double TOL2 = 0;
48:
49: #endregion
50:
51: public DLAZQ3(DLASQ5 dlasq5, DLASQ6 dlasq6, DLAZQ4 dlazq4, DLAMCH dlamch)
52: {
53:
54:
55: #region Set Dependencies
56:
57: this._dlasq5 = dlasq5; this._dlasq6 = dlasq6; this._dlazq4 = dlazq4; this._dlamch = dlamch;
58:
59: #endregion
60:
61: }
62:
63: public DLAZQ3()
64: {
65:
66:
67: #region Dependencies (Initialization)
68:
69: DLASQ5 dlasq5 = new DLASQ5();
70: LSAME lsame = new LSAME();
71: DLAMC3 dlamc3 = new DLAMC3();
72: DLAZQ4 dlazq4 = new DLAZQ4();
73: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
74: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
75: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
76: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
77: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
78: DLASQ6 dlasq6 = new DLASQ6(dlamch);
79:
80: #endregion
81:
82:
83: #region Set Dependencies
84:
85: this._dlasq5 = dlasq5; this._dlasq6 = dlasq6; this._dlazq4 = dlazq4; this._dlamch = dlamch;
86:
87: #endregion
88:
89: }
90: /// <summary>
91: /// Purpose
92: /// =======
93: ///
94: /// DLAZQ3 checks for deflation, computes a shift (TAU) and calls dqds.
95: /// In case of failure it changes shifts, and tries again until output
96: /// is positive.
97: ///
98: ///</summary>
99: /// <param name="I0">
100: /// (input) INTEGER
101: /// First index.
102: ///</param>
103: /// <param name="N0">
104: /// (input) INTEGER
105: /// Last index.
106: ///</param>
107: /// <param name="Z">
108: /// (input) DOUBLE PRECISION array, dimension ( 4*N )
109: /// Z holds the qd array.
110: ///</param>
111: /// <param name="PP">
112: /// (input) INTEGER
113: /// PP=0 for ping, PP=1 for pong.
114: ///</param>
115: /// <param name="DMIN">
116: /// (output) DOUBLE PRECISION
117: /// Minimum value of d.
118: ///</param>
119: /// <param name="SIGMA">
120: /// (output) DOUBLE PRECISION
121: /// Sum of shifts used in current segment.
122: ///</param>
123: /// <param name="DESIG">
124: /// (input/output) DOUBLE PRECISION
125: /// Lower order part of SIGMA
126: ///</param>
127: /// <param name="QMAX">
128: /// (input) DOUBLE PRECISION
129: /// Maximum value of q.
130: ///</param>
131: /// <param name="NFAIL">
132: /// (output) INTEGER
133: /// Number of times shift was too big.
134: ///</param>
135: /// <param name="ITER">
136: /// (output) INTEGER
137: /// Number of iterations.
138: ///</param>
139: /// <param name="NDIV">
140: /// (output) INTEGER
141: /// Number of divisions.
142: ///</param>
143: /// <param name="IEEE">
144: /// (input) LOGICAL
145: /// Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
146: ///</param>
147: /// <param name="TTYPE">
148: /// (input/output) INTEGER
149: /// Shift type. TTYPE is passed as an argument in order to save
150: /// its value between calls to DLAZQ3
151: ///</param>
152: /// <param name="DMIN1">
153: /// (input/output) REAL
154: ///</param>
155: /// <param name="DMIN2">
156: /// (input/output) REAL
157: ///</param>
158: /// <param name="DN">
159: /// (input/output) REAL
160: ///</param>
161: /// <param name="DN1">
162: /// (input/output) REAL
163: ///</param>
164: /// <param name="DN2">
165: /// (input/output) REAL
166: ///</param>
167: /// <param name="TAU">
168: /// (input/output) REAL
169: /// These are passed as arguments in order to save their values
170: /// between calls to DLAZQ3
171: ///</param>
172: public void Run(int I0, ref int N0, ref double[] Z, int offset_z, int PP, ref double DMIN, ref double SIGMA
173: , ref double DESIG, ref double QMAX, ref int NFAIL, ref int ITER, ref int NDIV, bool IEEE
174: , ref int TTYPE, ref double DMIN1, ref double DMIN2, ref double DN, ref double DN1, ref double DN2
175: , ref double TAU)
176: {
177:
178: #region Array Index Correction
179:
180: int o_z = -1 + offset_z;
181:
182: #endregion
183:
184:
185: #region Prolog
186:
187: // *
188: // * -- LAPACK auxiliary routine (version 3.1) --
189: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
190: // * November 2006
191: // *
192: // * .. Scalar Arguments ..
193: // * ..
194: // * .. Array Arguments ..
195: // * ..
196: // *
197: // * Purpose
198: // * =======
199: // *
200: // * DLAZQ3 checks for deflation, computes a shift (TAU) and calls dqds.
201: // * In case of failure it changes shifts, and tries again until output
202: // * is positive.
203: // *
204: // * Arguments
205: // * =========
206: // *
207: // * I0 (input) INTEGER
208: // * First index.
209: // *
210: // * N0 (input) INTEGER
211: // * Last index.
212: // *
213: // * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
214: // * Z holds the qd array.
215: // *
216: // * PP (input) INTEGER
217: // * PP=0 for ping, PP=1 for pong.
218: // *
219: // * DMIN (output) DOUBLE PRECISION
220: // * Minimum value of d.
221: // *
222: // * SIGMA (output) DOUBLE PRECISION
223: // * Sum of shifts used in current segment.
224: // *
225: // * DESIG (input/output) DOUBLE PRECISION
226: // * Lower order part of SIGMA
227: // *
228: // * QMAX (input) DOUBLE PRECISION
229: // * Maximum value of q.
230: // *
231: // * NFAIL (output) INTEGER
232: // * Number of times shift was too big.
233: // *
234: // * ITER (output) INTEGER
235: // * Number of iterations.
236: // *
237: // * NDIV (output) INTEGER
238: // * Number of divisions.
239: // *
240: // * IEEE (input) LOGICAL
241: // * Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
242: // *
243: // * TTYPE (input/output) INTEGER
244: // * Shift type. TTYPE is passed as an argument in order to save
245: // * its value between calls to DLAZQ3
246: // *
247: // * DMIN1 (input/output) REAL
248: // * DMIN2 (input/output) REAL
249: // * DN (input/output) REAL
250: // * DN1 (input/output) REAL
251: // * DN2 (input/output) REAL
252: // * TAU (input/output) REAL
253: // * These are passed as arguments in order to save their values
254: // * between calls to DLAZQ3
255: // *
256: // * This is a thread safe version of DLASQ3, which passes TTYPE, DMIN1,
257: // * DMIN2, DN, DN1. DN2 and TAU through the argument list in place of
258: // * declaring them in a SAVE statment.
259: // *
260: // * =====================================================================
261: // *
262: // * .. Parameters ..
263: // * ..
264: // * .. Local Scalars ..
265: // * ..
266: // * .. External Subroutines ..
267: // * ..
268: // * .. External Function ..
269: // * ..
270: // * .. Intrinsic Functions ..
271: // INTRINSIC ABS, MIN, SQRT;
272: // * ..
273: // * .. Executable Statements ..
274: // *
275:
276: #endregion
277:
278:
279: #region Body
280:
281: N0IN = N0;
282: EPS = this._dlamch.Run("Precision");
283: SAFMIN = this._dlamch.Run("Safe minimum");
284: TOL = EPS * HUNDRD;
285: TOL2 = Math.Pow(TOL,2);
286: G = ZERO;
287: // *
288: // * Check for deflation.
289: // *
290: LABEL10:;
291: // *
292: if (N0 < I0) return;
293: if (N0 == I0) goto LABEL20;
294: NN = 4 * N0 + PP;
295: if (N0 == (I0 + 1)) goto LABEL40;
296: // *
297: // * Check whether E(N0-1) is negligible, 1 eigenvalue.
298: // *
299: if (Z[NN - 5 + o_z] > TOL2 * (SIGMA + Z[NN - 3 + o_z]) && Z[NN - 2 * PP - 4 + o_z] > TOL2 * Z[NN - 7 + o_z]) goto LABEL30;
300: // *
301: LABEL20:;
302: // *
303: Z[4 * N0 - 3 + o_z] = Z[4 * N0 + PP - 3 + o_z] + SIGMA;
304: N0 = N0 - 1;
305: goto LABEL10;
306: // *
307: // * Check whether E(N0-2) is negligible, 2 eigenvalues.
308: // *
309: LABEL30:;
310: // *
311: if (Z[NN - 9 + o_z] > TOL2 * SIGMA && Z[NN - 2 * PP - 8 + o_z] > TOL2 * Z[NN - 11 + o_z]) goto LABEL50;
312: // *
313: LABEL40:;
314: // *
315: if (Z[NN - 3 + o_z] > Z[NN - 7 + o_z])
316: {
317: S = Z[NN - 3 + o_z];
318: Z[NN - 3 + o_z] = Z[NN - 7 + o_z];
319: Z[NN - 7 + o_z] = S;
320: }
321: if (Z[NN - 5 + o_z] > Z[NN - 3 + o_z] * TOL2)
322: {
323: T = HALF * ((Z[NN - 7 + o_z] - Z[NN - 3 + o_z]) + Z[NN - 5 + o_z]);
324: S = Z[NN - 3 + o_z] * (Z[NN - 5 + o_z] / T);
325: if (S <= T)
326: {
327: S = Z[NN - 3 + o_z] * (Z[NN - 5 + o_z] / (T * (ONE + Math.Sqrt(ONE + S / T))));
328: }
329: else
330: {
331: S = Z[NN - 3 + o_z] * (Z[NN - 5 + o_z] / (T + Math.Sqrt(T) * Math.Sqrt(T + S)));
332: }
333: T = Z[NN - 7 + o_z] + (S + Z[NN - 5 + o_z]);
334: Z[NN - 3 + o_z] = Z[NN - 3 + o_z] * (Z[NN - 7 + o_z] / T);
335: Z[NN - 7 + o_z] = T;
336: }
337: Z[4 * N0 - 7 + o_z] = Z[NN - 7 + o_z] + SIGMA;
338: Z[4 * N0 - 3 + o_z] = Z[NN - 3 + o_z] + SIGMA;
339: N0 = N0 - 2;
340: goto LABEL10;
341: // *
342: LABEL50:;
343: // *
344: // * Reverse the qd-array, if warranted.
345: // *
346: if (DMIN <= ZERO || N0 < N0IN)
347: {
348: if (CBIAS * Z[4 * I0 + PP - 3 + o_z] < Z[4 * N0 + PP - 3 + o_z])
349: {
350: IPN4 = 4 * (I0 + N0);
351: for (J4 = 4 * I0; J4 <= 2 * (I0 + N0 - 1); J4 += 4)
352: {
353: TEMP = Z[J4 - 3 + o_z];
354: Z[J4 - 3 + o_z] = Z[IPN4 - J4 - 3 + o_z];
355: Z[IPN4 - J4 - 3 + o_z] = TEMP;
356: TEMP = Z[J4 - 2 + o_z];
357: Z[J4 - 2 + o_z] = Z[IPN4 - J4 - 2 + o_z];
358: Z[IPN4 - J4 - 2 + o_z] = TEMP;
359: TEMP = Z[J4 - 1 + o_z];
360: Z[J4 - 1 + o_z] = Z[IPN4 - J4 - 5 + o_z];
361: Z[IPN4 - J4 - 5 + o_z] = TEMP;
362: TEMP = Z[J4 + o_z];
363: Z[J4 + o_z] = Z[IPN4 - J4 - 4 + o_z];
364: Z[IPN4 - J4 - 4 + o_z] = TEMP;
365: }
366: if (N0 - I0 <= 4)
367: {
368: Z[4 * N0 + PP - 1 + o_z] = Z[4 * I0 + PP - 1 + o_z];
369: Z[4 * N0 - PP + o_z] = Z[4 * I0 - PP + o_z];
370: }
371: DMIN2 = Math.Min(DMIN2, Z[4 * N0 + PP - 1 + o_z]);
372: Z[4 * N0 + PP - 1 + o_z] = Math.Min(Z[4 * N0 + PP - 1 + o_z], Math.Min(Z[4 * I0 + PP - 1 + o_z], Z[4 * I0 + PP + 3 + o_z]));
373: Z[4 * N0 - PP + o_z] = Math.Min(Z[4 * N0 - PP + o_z], Math.Min(Z[4 * I0 - PP + o_z], Z[4 * I0 - PP + 4 + o_z]));
374: QMAX = Math.Max(QMAX, Math.Max(Z[4 * I0 + PP - 3 + o_z], Z[4 * I0 + PP + 1 + o_z]));
375: DMIN = - ZERO;
376: }
377: }
378: // *
379: if (DMIN < ZERO || SAFMIN * QMAX < Math.Min(Z[4 * N0 + PP - 1 + o_z], Math.Min(Z[4 * N0 + PP - 9 + o_z], DMIN2 + Z[4 * N0 - PP + o_z])))
380: {
381: // *
382: // * Choose a shift.
383: // *
384: this._dlazq4.Run(I0, N0, Z, offset_z, PP, N0IN, DMIN
385: , DMIN1, DMIN2, DN, DN1, DN2, ref TAU
386: , ref TTYPE, ref G);
387: // *
388: // * Call dqds until DMIN > 0.
389: // *
390: LABEL80:;
391: // *
392: this._dlasq5.Run(I0, N0, ref Z, offset_z, PP, TAU, ref DMIN
393: , ref DMIN1, ref DMIN2, ref DN, ref DN1, ref DN2, IEEE);
394: // *
395: NDIV = NDIV + (N0 - I0 + 2);
396: ITER = ITER + 1;
397: // *
398: // * Check status.
399: // *
400: if (DMIN >= ZERO && DMIN1 > ZERO)
401: {
402: // *
403: // * Success.
404: // *
405: goto LABEL100;
406: // *
407: }
408: else
409: {
410: if (DMIN < ZERO && DMIN1 > ZERO && Z[4 * (N0 - 1) - PP + o_z] < TOL * (SIGMA + DN1) && Math.Abs(DN) < TOL * SIGMA)
411: {
412: // *
413: // * Convergence hidden by negative DN.
414: // *
415: Z[4 * (N0 - 1) - PP + 2 + o_z] = ZERO;
416: DMIN = ZERO;
417: goto LABEL100;
418: }
419: else
420: {
421: if (DMIN < ZERO)
422: {
423: // *
424: // * TAU too big. Select new TAU and try again.
425: // *
426: NFAIL = NFAIL + 1;
427: if (TTYPE < - 22)
428: {
429: // *
430: // * Failed twice. Play it safe.
431: // *
432: TAU = ZERO;
433: }
434: else
435: {
436: if (DMIN1 > ZERO)
437: {
438: // *
439: // * Late failure. Gives excellent shift.
440: // *
441: TAU = (TAU + DMIN) * (ONE - TWO * EPS);
442: TTYPE = TTYPE - 11;
443: }
444: else
445: {
446: // *
447: // * Early failure. Divide by 4.
448: // *
449: TAU = QURTR * TAU;
450: TTYPE = TTYPE - 12;
451: }
452: }
453: goto LABEL80;
454: }
455: else
456: {
457: if (DMIN != DMIN)
458: {
459: // *
460: // * NaN.
461: // *
462: TAU = ZERO;
463: goto LABEL80;
464: }
465: else
466: {
467: // *
468: // * Possible underflow. Play it safe.
469: // *
470: goto LABEL90;
471: }
472: }
473: }
474: }
475: }
476: // *
477: // * Risk of underflow.
478: // *
479: LABEL90:;
480: this._dlasq6.Run(I0, N0, ref Z, offset_z, PP, ref DMIN, ref DMIN1
481: , ref DMIN2, ref DN, ref DN1, ref DN2);
482: NDIV = NDIV + (N0 - I0 + 2);
483: ITER = ITER + 1;
484: TAU = ZERO;
485: // *
486: LABEL100:;
487: if (TAU < SIGMA)
488: {
489: DESIG = DESIG + TAU;
490: T = SIGMA + DESIG;
491: DESIG = DESIG - (T - SIGMA);
492: }
493: else
494: {
495: T = SIGMA + TAU;
496: DESIG = SIGMA - (T - TAU) + DESIG;
497: }
498: SIGMA = T;
499: // *
500: return;
501: // *
502: // * End of DLAZQ3
503: // *
504:
505: #endregion
506:
507: }
508: }
509: }