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: /// DLASQ5 computes one dqds transform in ping-pong form, one
27: /// version for IEEE machines another for non IEEE machines.
28: ///
29: ///</summary>
30: public class DLASQ5
31: {
32:
33:
34: #region Fields
35:
36: const double ZERO = 0.0E0; int J4 = 0; int J4P2 = 0; double D = 0; double EMIN = 0; double TEMP = 0;
37:
38: #endregion
39:
40: public DLASQ5()
41: {
42:
43: }
44:
45: /// <summary>
46: /// Purpose
47: /// =======
48: ///
49: /// DLASQ5 computes one dqds transform in ping-pong form, one
50: /// version for IEEE machines another for non IEEE machines.
51: ///
52: ///</summary>
53: /// <param name="I0">
54: /// (input) INTEGER
55: /// First index.
56: ///</param>
57: /// <param name="N0">
58: /// (input) INTEGER
59: /// Last index.
60: ///</param>
61: /// <param name="Z">
62: /// (input) DOUBLE PRECISION array, dimension ( 4*N )
63: /// Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
64: /// an extra argument.
65: ///</param>
66: /// <param name="PP">
67: /// (input) INTEGER
68: /// PP=0 for ping, PP=1 for pong.
69: ///</param>
70: /// <param name="TAU">
71: /// (input) DOUBLE PRECISION
72: /// This is the shift.
73: ///</param>
74: /// <param name="DMIN">
75: /// (output) DOUBLE PRECISION
76: /// Minimum value of d.
77: ///</param>
78: /// <param name="DMIN1">
79: /// (output) DOUBLE PRECISION
80: /// Minimum value of d, excluding D( N0 ).
81: ///</param>
82: /// <param name="DMIN2">
83: /// (output) DOUBLE PRECISION
84: /// Minimum value of d, excluding D( N0 ) and D( N0-1 ).
85: ///</param>
86: /// <param name="DN">
87: /// (output) DOUBLE PRECISION
88: /// d(N0), the last value of d.
89: ///</param>
90: /// <param name="DNM1">
91: /// (output) DOUBLE PRECISION
92: /// d(N0-1).
93: ///</param>
94: /// <param name="DNM2">
95: /// (output) DOUBLE PRECISION
96: /// d(N0-2).
97: ///</param>
98: /// <param name="IEEE">
99: /// (input) LOGICAL
100: /// Flag for IEEE or non IEEE arithmetic.
101: ///</param>
102: public void Run(int I0, int N0, ref double[] Z, int offset_z, int PP, double TAU, ref double DMIN
103: , ref double DMIN1, ref double DMIN2, ref double DN, ref double DNM1, ref double DNM2, bool IEEE)
104: {
105:
106: #region Array Index Correction
107:
108: int o_z = -1 + offset_z;
109:
110: #endregion
111:
112:
113: #region Prolog
114:
115: // *
116: // * -- LAPACK auxiliary routine (version 3.1) --
117: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
118: // * November 2006
119: // *
120: // * .. Scalar Arguments ..
121: // * ..
122: // * .. Array Arguments ..
123: // * ..
124: // *
125: // * Purpose
126: // * =======
127: // *
128: // * DLASQ5 computes one dqds transform in ping-pong form, one
129: // * version for IEEE machines another for non IEEE machines.
130: // *
131: // * Arguments
132: // * =========
133: // *
134: // * I0 (input) INTEGER
135: // * First index.
136: // *
137: // * N0 (input) INTEGER
138: // * Last index.
139: // *
140: // * Z (input) DOUBLE PRECISION array, dimension ( 4*N )
141: // * Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
142: // * an extra argument.
143: // *
144: // * PP (input) INTEGER
145: // * PP=0 for ping, PP=1 for pong.
146: // *
147: // * TAU (input) DOUBLE PRECISION
148: // * This is the shift.
149: // *
150: // * DMIN (output) DOUBLE PRECISION
151: // * Minimum value of d.
152: // *
153: // * DMIN1 (output) DOUBLE PRECISION
154: // * Minimum value of d, excluding D( N0 ).
155: // *
156: // * DMIN2 (output) DOUBLE PRECISION
157: // * Minimum value of d, excluding D( N0 ) and D( N0-1 ).
158: // *
159: // * DN (output) DOUBLE PRECISION
160: // * d(N0), the last value of d.
161: // *
162: // * DNM1 (output) DOUBLE PRECISION
163: // * d(N0-1).
164: // *
165: // * DNM2 (output) DOUBLE PRECISION
166: // * d(N0-2).
167: // *
168: // * IEEE (input) LOGICAL
169: // * Flag for IEEE or non IEEE arithmetic.
170: // *
171: // * =====================================================================
172: // *
173: // * .. Parameter ..
174: // * ..
175: // * .. Local Scalars ..
176: // * ..
177: // * .. Intrinsic Functions ..
178: // INTRINSIC MIN;
179: // * ..
180: // * .. Executable Statements ..
181: // *
182:
183: #endregion
184:
185:
186: #region Body
187:
188: if ((N0 - I0 - 1) <= 0) return;
189: // *
190: J4 = 4 * I0 + PP - 3;
191: EMIN = Z[J4 + 4 + o_z];
192: D = Z[J4 + o_z] - TAU;
193: DMIN = D;
194: DMIN1 = - Z[J4 + o_z];
195: // *
196: if (IEEE)
197: {
198: // *
199: // * Code for IEEE arithmetic.
200: // *
201: if (PP == 0)
202: {
203: for (J4 = 4 * I0; J4 <= 4 * (N0 - 3); J4 += 4)
204: {
205: Z[J4 - 2 + o_z] = D + Z[J4 - 1 + o_z];
206: TEMP = Z[J4 + 1 + o_z] / Z[J4 - 2 + o_z];
207: D = D * TEMP - TAU;
208: DMIN = Math.Min(DMIN, D);
209: Z[J4 + o_z] = Z[J4 - 1 + o_z] * TEMP;
210: EMIN = Math.Min(Z[J4 + o_z], EMIN);
211: }
212: }
213: else
214: {
215: for (J4 = 4 * I0; J4 <= 4 * (N0 - 3); J4 += 4)
216: {
217: Z[J4 - 3 + o_z] = D + Z[J4 + o_z];
218: TEMP = Z[J4 + 2 + o_z] / Z[J4 - 3 + o_z];
219: D = D * TEMP - TAU;
220: DMIN = Math.Min(DMIN, D);
221: Z[J4 - 1 + o_z] = Z[J4 + o_z] * TEMP;
222: EMIN = Math.Min(Z[J4 - 1 + o_z], EMIN);
223: }
224: }
225: // *
226: // * Unroll last two steps.
227: // *
228: DNM2 = D;
229: DMIN2 = DMIN;
230: J4 = 4 * (N0 - 2) - PP;
231: J4P2 = J4 + 2 * PP - 1;
232: Z[J4 - 2 + o_z] = DNM2 + Z[J4P2 + o_z];
233: Z[J4 + o_z] = Z[J4P2 + 2 + o_z] * (Z[J4P2 + o_z] / Z[J4 - 2 + o_z]);
234: DNM1 = Z[J4P2 + 2 + o_z] * (DNM2 / Z[J4 - 2 + o_z]) - TAU;
235: DMIN = Math.Min(DMIN, DNM1);
236: // *
237: DMIN1 = DMIN;
238: J4 = J4 + 4;
239: J4P2 = J4 + 2 * PP - 1;
240: Z[J4 - 2 + o_z] = DNM1 + Z[J4P2 + o_z];
241: Z[J4 + o_z] = Z[J4P2 + 2 + o_z] * (Z[J4P2 + o_z] / Z[J4 - 2 + o_z]);
242: DN = Z[J4P2 + 2 + o_z] * (DNM1 / Z[J4 - 2 + o_z]) - TAU;
243: DMIN = Math.Min(DMIN, DN);
244: // *
245: }
246: else
247: {
248: // *
249: // * Code for non IEEE arithmetic.
250: // *
251: if (PP == 0)
252: {
253: for (J4 = 4 * I0; J4 <= 4 * (N0 - 3); J4 += 4)
254: {
255: Z[J4 - 2 + o_z] = D + Z[J4 - 1 + o_z];
256: if (D < ZERO)
257: {
258: return;
259: }
260: else
261: {
262: Z[J4 + o_z] = Z[J4 + 1 + o_z] * (Z[J4 - 1 + o_z] / Z[J4 - 2 + o_z]);
263: D = Z[J4 + 1 + o_z] * (D / Z[J4 - 2 + o_z]) - TAU;
264: }
265: DMIN = Math.Min(DMIN, D);
266: EMIN = Math.Min(EMIN, Z[J4 + o_z]);
267: }
268: }
269: else
270: {
271: for (J4 = 4 * I0; J4 <= 4 * (N0 - 3); J4 += 4)
272: {
273: Z[J4 - 3 + o_z] = D + Z[J4 + o_z];
274: if (D < ZERO)
275: {
276: return;
277: }
278: else
279: {
280: Z[J4 - 1 + o_z] = Z[J4 + 2 + o_z] * (Z[J4 + o_z] / Z[J4 - 3 + o_z]);
281: D = Z[J4 + 2 + o_z] * (D / Z[J4 - 3 + o_z]) - TAU;
282: }
283: DMIN = Math.Min(DMIN, D);
284: EMIN = Math.Min(EMIN, Z[J4 - 1 + o_z]);
285: }
286: }
287: // *
288: // * Unroll last two steps.
289: // *
290: DNM2 = D;
291: DMIN2 = DMIN;
292: J4 = 4 * (N0 - 2) - PP;
293: J4P2 = J4 + 2 * PP - 1;
294: Z[J4 - 2 + o_z] = DNM2 + Z[J4P2 + o_z];
295: if (DNM2 < ZERO)
296: {
297: return;
298: }
299: else
300: {
301: Z[J4 + o_z] = Z[J4P2 + 2 + o_z] * (Z[J4P2 + o_z] / Z[J4 - 2 + o_z]);
302: DNM1 = Z[J4P2 + 2 + o_z] * (DNM2 / Z[J4 - 2 + o_z]) - TAU;
303: }
304: DMIN = Math.Min(DMIN, DNM1);
305: // *
306: DMIN1 = DMIN;
307: J4 = J4 + 4;
308: J4P2 = J4 + 2 * PP - 1;
309: Z[J4 - 2 + o_z] = DNM1 + Z[J4P2 + o_z];
310: if (DNM1 < ZERO)
311: {
312: return;
313: }
314: else
315: {
316: Z[J4 + o_z] = Z[J4P2 + 2 + o_z] * (Z[J4P2 + o_z] / Z[J4 - 2 + o_z]);
317: DN = Z[J4P2 + 2 + o_z] * (DNM1 / Z[J4 - 2 + o_z]) - TAU;
318: }
319: DMIN = Math.Min(DMIN, DN);
320: // *
321: }
322: // *
323: Z[J4 + 2 + o_z] = DN;
324: Z[4 * N0 - PP + o_z] = EMIN;
325: return;
326: // *
327: // * End of DLASQ5
328: // *
329:
330: #endregion
331:
332: }
333: }
334: }