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