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: /// DLAZQ4 computes an approximation TAU to the smallest eigenvalue
27: /// using values of d from the previous transform.
28: ///
29: /// I0 (input) INTEGER
30: /// First index.
31: ///
32: /// N0 (input) INTEGER
33: /// Last index.
34: ///
35: /// Z (input) DOUBLE PRECISION array, dimension ( 4*N )
36: /// Z holds the qd array.
37: ///
38: /// PP (input) INTEGER
39: /// PP=0 for ping, PP=1 for pong.
40: ///
41: /// N0IN (input) INTEGER
42: /// The value of N0 at start of EIGTEST.
43: ///
44: /// DMIN (input) DOUBLE PRECISION
45: /// Minimum value of d.
46: ///
47: /// DMIN1 (input) DOUBLE PRECISION
48: /// Minimum value of d, excluding D( N0 ).
49: ///
50: /// DMIN2 (input) DOUBLE PRECISION
51: /// Minimum value of d, excluding D( N0 ) and D( N0-1 ).
52: ///
53: /// DN (input) DOUBLE PRECISION
54: /// d(N)
55: ///
56: /// DN1 (input) DOUBLE PRECISION
57: /// d(N-1)
58: ///
59: /// DN2 (input) DOUBLE PRECISION
60: /// d(N-2)
61: ///
62: /// TAU (output) DOUBLE PRECISION
63: /// This is the shift.
64: ///
65: /// TTYPE (output) INTEGER
66: /// Shift type.
67: ///
68: /// G (input/output) DOUBLE PRECISION
69: /// G is passed as an argument in order to save its value between
70: /// calls to DLAZQ4
71: ///
72: ///</summary>
73: public class DLAZQ4
74: {
75:
76:
77: #region Fields
78:
79: const double CNST1 = 0.5630E0; const double CNST2 = 1.010E0; const double CNST3 = 1.050E0; const double QURTR = 0.250E0;
80: const double THIRD = 0.3330E0;const double HALF = 0.50E0; const double ZERO = 0.0E0; const double ONE = 1.0E0;
81: const double TWO = 2.0E0;const double HUNDRD = 100.0E0; int I4 = 0; int NN = 0; int NP = 0; double A2 = 0; double B1 = 0;
82: double B2 = 0;double GAM = 0; double GAP1 = 0; double GAP2 = 0; double S = 0;
83:
84: #endregion
85:
86: public DLAZQ4()
87: {
88:
89: }
90:
91: /// <summary>
92: /// Purpose
93: /// =======
94: ///
95: /// DLAZQ4 computes an approximation TAU to the smallest eigenvalue
96: /// using values of d from the previous transform.
97: ///
98: /// I0 (input) INTEGER
99: /// First index.
100: ///
101: /// N0 (input) INTEGER
102: /// Last index.
103: ///
104: /// Z (input) DOUBLE PRECISION array, dimension ( 4*N )
105: /// Z holds the qd array.
106: ///
107: /// PP (input) INTEGER
108: /// PP=0 for ping, PP=1 for pong.
109: ///
110: /// N0IN (input) INTEGER
111: /// The value of N0 at start of EIGTEST.
112: ///
113: /// DMIN (input) DOUBLE PRECISION
114: /// Minimum value of d.
115: ///
116: /// DMIN1 (input) DOUBLE PRECISION
117: /// Minimum value of d, excluding D( N0 ).
118: ///
119: /// DMIN2 (input) DOUBLE PRECISION
120: /// Minimum value of d, excluding D( N0 ) and D( N0-1 ).
121: ///
122: /// DN (input) DOUBLE PRECISION
123: /// d(N)
124: ///
125: /// DN1 (input) DOUBLE PRECISION
126: /// d(N-1)
127: ///
128: /// DN2 (input) DOUBLE PRECISION
129: /// d(N-2)
130: ///
131: /// TAU (output) DOUBLE PRECISION
132: /// This is the shift.
133: ///
134: /// TTYPE (output) INTEGER
135: /// Shift type.
136: ///
137: /// G (input/output) DOUBLE PRECISION
138: /// G is passed as an argument in order to save its value between
139: /// calls to DLAZQ4
140: ///
141: ///</summary>
142: /// <param name="I0">
143: /// (input) INTEGER
144: /// First index.
145: ///</param>
146: /// <param name="N0">
147: /// (input) INTEGER
148: /// Last index.
149: ///</param>
150: /// <param name="Z">
151: /// (input) DOUBLE PRECISION array, dimension ( 4*N )
152: /// Z holds the qd array.
153: ///</param>
154: /// <param name="PP">
155: /// (input) INTEGER
156: /// PP=0 for ping, PP=1 for pong.
157: ///</param>
158: /// <param name="N0IN">
159: /// (input) INTEGER
160: /// The value of N0 at start of EIGTEST.
161: ///</param>
162: /// <param name="DMIN">
163: /// (input) DOUBLE PRECISION
164: /// Minimum value of d.
165: ///</param>
166: /// <param name="DMIN1">
167: /// (input) DOUBLE PRECISION
168: /// Minimum value of d, excluding D( N0 ).
169: ///</param>
170: /// <param name="DMIN2">
171: /// (input) DOUBLE PRECISION
172: /// Minimum value of d, excluding D( N0 ) and D( N0-1 ).
173: ///</param>
174: /// <param name="DN">
175: /// (input) DOUBLE PRECISION
176: /// d(N)
177: ///</param>
178: /// <param name="DN1">
179: /// (input) DOUBLE PRECISION
180: /// d(N-1)
181: ///</param>
182: /// <param name="DN2">
183: /// (input) DOUBLE PRECISION
184: /// d(N-2)
185: ///</param>
186: /// <param name="TAU">
187: /// (output) DOUBLE PRECISION
188: /// This is the shift.
189: ///</param>
190: /// <param name="TTYPE">
191: /// (output) INTEGER
192: /// Shift type.
193: ///</param>
194: /// <param name="G">
195: /// (input/output) DOUBLE PRECISION
196: /// G is passed as an argument in order to save its value between
197: /// calls to DLAZQ4
198: ///</param>
199: public void Run(int I0, int N0, double[] Z, int offset_z, int PP, int N0IN, double DMIN
200: , double DMIN1, double DMIN2, double DN, double DN1, double DN2, ref double TAU
201: , ref int TTYPE, ref double G)
202: {
203:
204: #region Array Index Correction
205:
206: int o_z = -1 + offset_z;
207:
208: #endregion
209:
210:
211: #region Prolog
212:
213: // *
214: // * -- LAPACK auxiliary routine (version 3.1) --
215: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
216: // * November 2006
217: // *
218: // * .. Scalar Arguments ..
219: // * ..
220: // * .. Array Arguments ..
221: // * ..
222: // *
223: // * Purpose
224: // * =======
225: // *
226: // * DLAZQ4 computes an approximation TAU to the smallest eigenvalue
227: // * using values of d from the previous transform.
228: // *
229: // * I0 (input) INTEGER
230: // * First index.
231: // *
232: // * N0 (input) INTEGER
233: // * Last index.
234: // *
235: // * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
236: // * Z holds the qd array.
237: // *
238: // * PP (input) INTEGER
239: // * PP=0 for ping, PP=1 for pong.
240: // *
241: // * N0IN (input) INTEGER
242: // * The value of N0 at start of EIGTEST.
243: // *
244: // * DMIN (input) DOUBLE PRECISION
245: // * Minimum value of d.
246: // *
247: // * DMIN1 (input) DOUBLE PRECISION
248: // * Minimum value of d, excluding D( N0 ).
249: // *
250: // * DMIN2 (input) DOUBLE PRECISION
251: // * Minimum value of d, excluding D( N0 ) and D( N0-1 ).
252: // *
253: // * DN (input) DOUBLE PRECISION
254: // * d(N)
255: // *
256: // * DN1 (input) DOUBLE PRECISION
257: // * d(N-1)
258: // *
259: // * DN2 (input) DOUBLE PRECISION
260: // * d(N-2)
261: // *
262: // * TAU (output) DOUBLE PRECISION
263: // * This is the shift.
264: // *
265: // * TTYPE (output) INTEGER
266: // * Shift type.
267: // *
268: // * G (input/output) DOUBLE PRECISION
269: // * G is passed as an argument in order to save its value between
270: // * calls to DLAZQ4
271: // *
272: // * Further Details
273: // * ===============
274: // * CNST1 = 9/16
275: // *
276: // * This is a thread safe version of DLASQ4, which passes G through the
277: // * argument list in place of declaring G in a SAVE statment.
278: // *
279: // * =====================================================================
280: // *
281: // * .. Parameters ..
282: // * ..
283: // * .. Local Scalars ..
284: // * ..
285: // * .. Intrinsic Functions ..
286: // INTRINSIC MAX, MIN, SQRT;
287: // * ..
288: // * .. Executable Statements ..
289: // *
290: // * A negative DMIN forces the shift to take that absolute value
291: // * TTYPE records the type of shift.
292: // *
293:
294: #endregion
295:
296:
297: #region Body
298:
299: if (DMIN <= ZERO)
300: {
301: TAU = - DMIN;
302: TTYPE = - 1;
303: return;
304: }
305: // *
306: NN = 4 * N0 + PP;
307: if (N0IN == N0)
308: {
309: // *
310: // * No eigenvalues deflated.
311: // *
312: if (DMIN == DN || DMIN == DN1)
313: {
314: // *
315: B1 = Math.Sqrt(Z[NN - 3 + o_z]) * Math.Sqrt(Z[NN - 5 + o_z]);
316: B2 = Math.Sqrt(Z[NN - 7 + o_z]) * Math.Sqrt(Z[NN - 9 + o_z]);
317: A2 = Z[NN - 7 + o_z] + Z[NN - 5 + o_z];
318: // *
319: // * Cases 2 and 3.
320: // *
321: if (DMIN == DN && DMIN1 == DN1)
322: {
323: GAP2 = DMIN2 - A2 - DMIN2 * QURTR;
324: if (GAP2 > ZERO && GAP2 > B2)
325: {
326: GAP1 = A2 - DN - (B2 / GAP2) * B2;
327: }
328: else
329: {
330: GAP1 = A2 - DN - (B1 + B2);
331: }
332: if (GAP1 > ZERO && GAP1 > B1)
333: {
334: S = Math.Max(DN - (B1 / GAP1) * B1, HALF * DMIN);
335: TTYPE = - 2;
336: }
337: else
338: {
339: S = ZERO;
340: if (DN > B1) S = DN - B1;
341: if (A2 > (B1 + B2)) S = Math.Min(S, A2 - (B1 + B2));
342: S = Math.Max(S, THIRD * DMIN);
343: TTYPE = - 3;
344: }
345: }
346: else
347: {
348: // *
349: // * Case 4.
350: // *
351: TTYPE = - 4;
352: S = QURTR * DMIN;
353: if (DMIN == DN)
354: {
355: GAM = DN;
356: A2 = ZERO;
357: if (Z[NN - 5 + o_z] > Z[NN - 7 + o_z]) return;
358: B2 = Z[NN - 5 + o_z] / Z[NN - 7 + o_z];
359: NP = NN - 9;
360: }
361: else
362: {
363: NP = NN - 2 * PP;
364: B2 = Z[NP - 2 + o_z];
365: GAM = DN1;
366: if (Z[NP - 4 + o_z] > Z[NP - 2 + o_z]) return;
367: A2 = Z[NP - 4 + o_z] / Z[NP - 2 + o_z];
368: if (Z[NN - 9 + o_z] > Z[NN - 11 + o_z]) return;
369: B2 = Z[NN - 9 + o_z] / Z[NN - 11 + o_z];
370: NP = NN - 13;
371: }
372: // *
373: // * Approximate contribution to norm squared from I < NN-1.
374: // *
375: A2 = A2 + B2;
376: for (I4 = NP; I4 >= 4 * I0 - 1 + PP; I4 += - 4)
377: {
378: if (B2 == ZERO) goto LABEL20;
379: B1 = B2;
380: if (Z[I4 + o_z] > Z[I4 - 2 + o_z]) return;
381: B2 = B2 * (Z[I4 + o_z] / Z[I4 - 2 + o_z]);
382: A2 = A2 + B2;
383: if (HUNDRD * Math.Max(B2, B1) < A2 || CNST1 < A2) goto LABEL20;
384: }
385: LABEL20:;
386: A2 = CNST3 * A2;
387: // *
388: // * Rayleigh quotient residual bound.
389: // *
390: if (A2 < CNST1) S = GAM * (ONE - Math.Sqrt(A2)) / (ONE + A2);
391: }
392: }
393: else
394: {
395: if (DMIN == DN2)
396: {
397: // *
398: // * Case 5.
399: // *
400: TTYPE = - 5;
401: S = QURTR * DMIN;
402: // *
403: // * Compute contribution to norm squared from I > NN-2.
404: // *
405: NP = NN - 2 * PP;
406: B1 = Z[NP - 2 + o_z];
407: B2 = Z[NP - 6 + o_z];
408: GAM = DN2;
409: if (Z[NP - 8 + o_z] > B2 || Z[NP - 4 + o_z] > B1) return;
410: A2 = (Z[NP - 8 + o_z] / B2) * (ONE + Z[NP - 4 + o_z] / B1);
411: // *
412: // * Approximate contribution to norm squared from I < NN-2.
413: // *
414: if (N0 - I0 > 2)
415: {
416: B2 = Z[NN - 13 + o_z] / Z[NN - 15 + o_z];
417: A2 = A2 + B2;
418: for (I4 = NN - 17; I4 >= 4 * I0 - 1 + PP; I4 += - 4)
419: {
420: if (B2 == ZERO) goto LABEL40;
421: B1 = B2;
422: if (Z[I4 + o_z] > Z[I4 - 2 + o_z]) return;
423: B2 = B2 * (Z[I4 + o_z] / Z[I4 - 2 + o_z]);
424: A2 = A2 + B2;
425: if (HUNDRD * Math.Max(B2, B1) < A2 || CNST1 < A2) goto LABEL40;
426: }
427: LABEL40:;
428: A2 = CNST3 * A2;
429: }
430: // *
431: if (A2 < CNST1) S = GAM * (ONE - Math.Sqrt(A2)) / (ONE + A2);
432: }
433: else
434: {
435: // *
436: // * Case 6, no information to guide us.
437: // *
438: if (TTYPE == - 6)
439: {
440: G = G + THIRD * (ONE - G);
441: }
442: else
443: {
444: if (TTYPE == - 18)
445: {
446: G = QURTR * THIRD;
447: }
448: else
449: {
450: G = QURTR;
451: }
452: }
453: S = G * DMIN;
454: TTYPE = - 6;
455: }
456: }
457: // *
458: }
459: else
460: {
461: if (N0IN == (N0 + 1))
462: {
463: // *
464: // * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
465: // *
466: if (DMIN1 == DN1 && DMIN2 == DN2)
467: {
468: // *
469: // * Cases 7 and 8.
470: // *
471: TTYPE = - 7;
472: S = THIRD * DMIN1;
473: if (Z[NN - 5 + o_z] > Z[NN - 7 + o_z]) return;
474: B1 = Z[NN - 5 + o_z] / Z[NN - 7 + o_z];
475: B2 = B1;
476: if (B2 == ZERO) goto LABEL60;
477: for (I4 = 4 * N0 - 9 + PP; I4 >= 4 * I0 - 1 + PP; I4 += - 4)
478: {
479: A2 = B1;
480: if (Z[I4 + o_z] > Z[I4 - 2 + o_z]) return;
481: B1 = B1 * (Z[I4 + o_z] / Z[I4 - 2 + o_z]);
482: B2 = B2 + B1;
483: if (HUNDRD * Math.Max(B1, A2) < B2) goto LABEL60;
484: }
485: LABEL60:;
486: B2 = Math.Sqrt(CNST3 * B2);
487: A2 = DMIN1 / (ONE + Math.Pow(B2,2));
488: GAP2 = HALF * DMIN2 - A2;
489: if (GAP2 > ZERO && GAP2 > B2 * A2)
490: {
491: S = Math.Max(S, A2 * (ONE - CNST2 * A2 * (B2 / GAP2) * B2));
492: }
493: else
494: {
495: S = Math.Max(S, A2 * (ONE - CNST2 * B2));
496: TTYPE = - 8;
497: }
498: }
499: else
500: {
501: // *
502: // * Case 9.
503: // *
504: S = QURTR * DMIN1;
505: if (DMIN1 == DN1) S = HALF * DMIN1;
506: TTYPE = - 9;
507: }
508: // *
509: }
510: else
511: {
512: if (N0IN == (N0 + 2))
513: {
514: // *
515: // * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
516: // *
517: // * Cases 10 and 11.
518: // *
519: if (DMIN2 == DN2 && TWO * Z[NN - 5 + o_z] < Z[NN - 7 + o_z])
520: {
521: TTYPE = - 10;
522: S = THIRD * DMIN2;
523: if (Z[NN - 5 + o_z] > Z[NN - 7 + o_z]) return;
524: B1 = Z[NN - 5 + o_z] / Z[NN - 7 + o_z];
525: B2 = B1;
526: if (B2 == ZERO) goto LABEL80;
527: for (I4 = 4 * N0 - 9 + PP; I4 >= 4 * I0 - 1 + PP; I4 += - 4)
528: {
529: if (Z[I4 + o_z] > Z[I4 - 2 + o_z]) return;
530: B1 = B1 * (Z[I4 + o_z] / Z[I4 - 2 + o_z]);
531: B2 = B2 + B1;
532: if (HUNDRD * B1 < B2) goto LABEL80;
533: }
534: LABEL80:;
535: B2 = Math.Sqrt(CNST3 * B2);
536: A2 = DMIN2 / (ONE + Math.Pow(B2,2));
537: GAP2 = Z[NN - 7 + o_z] + Z[NN - 9 + o_z] - Math.Sqrt(Z[NN - 11 + o_z]) * Math.Sqrt(Z[NN - 9 + o_z]) - A2;
538: if (GAP2 > ZERO && GAP2 > B2 * A2)
539: {
540: S = Math.Max(S, A2 * (ONE - CNST2 * A2 * (B2 / GAP2) * B2));
541: }
542: else
543: {
544: S = Math.Max(S, A2 * (ONE - CNST2 * B2));
545: }
546: }
547: else
548: {
549: S = QURTR * DMIN2;
550: TTYPE = - 11;
551: }
552: }
553: else
554: {
555: if (N0IN > (N0 + 2))
556: {
557: // *
558: // * Case 12, more than two eigenvalues deflated. No information.
559: // *
560: S = ZERO;
561: TTYPE = - 12;
562: }
563: }
564: }
565: }
566: // *
567: TAU = S;
568: return;
569: // *
570: // * End of DLAZQ4
571: // *
572:
573: #endregion
574:
575: }
576: }
577: }