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: /// DLARFG generates a real elementary reflector H of order n, such
27: /// that
28: ///
29: /// H * ( alpha ) = ( beta ), H' * H = I.
30: /// ( x ) ( 0 )
31: ///
32: /// where alpha and beta are scalars, and x is an (n-1)-element real
33: /// vector. H is represented in the form
34: ///
35: /// H = I - tau * ( 1 ) * ( 1 v' ) ,
36: /// ( v )
37: ///
38: /// where tau is a real scalar and v is a real (n-1)-element
39: /// vector.
40: ///
41: /// If the elements of x are all zero, then tau = 0 and H is taken to be
42: /// the unit matrix.
43: ///
44: /// Otherwise 1 .LE. tau .LE. 2.
45: ///
46: ///</summary>
47: public class DLARFG
48: {
49:
50:
51: #region Dependencies
52:
53: DLAMCH _dlamch; DLAPY2 _dlapy2; DNRM2 _dnrm2; DSCAL _dscal;
54:
55: #endregion
56:
57:
58: #region Fields
59:
60: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; int J = 0; int KNT = 0; double BETA = 0; double RSAFMN = 0;
61: double SAFMIN = 0;double XNORM = 0;
62:
63: #endregion
64:
65: public DLARFG(DLAMCH dlamch, DLAPY2 dlapy2, DNRM2 dnrm2, DSCAL dscal)
66: {
67:
68:
69: #region Set Dependencies
70:
71: this._dlamch = dlamch; this._dlapy2 = dlapy2; this._dnrm2 = dnrm2; this._dscal = dscal;
72:
73: #endregion
74:
75: }
76:
77: public DLARFG()
78: {
79:
80:
81: #region Dependencies (Initialization)
82:
83: LSAME lsame = new LSAME();
84: DLAMC3 dlamc3 = new DLAMC3();
85: DLAPY2 dlapy2 = new DLAPY2();
86: DNRM2 dnrm2 = new DNRM2();
87: DSCAL dscal = new DSCAL();
88: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
89: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
90: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
91: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
92: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
93:
94: #endregion
95:
96:
97: #region Set Dependencies
98:
99: this._dlamch = dlamch; this._dlapy2 = dlapy2; this._dnrm2 = dnrm2; this._dscal = dscal;
100:
101: #endregion
102:
103: }
104: /// <summary>
105: /// Purpose
106: /// =======
107: ///
108: /// DLARFG generates a real elementary reflector H of order n, such
109: /// that
110: ///
111: /// H * ( alpha ) = ( beta ), H' * H = I.
112: /// ( x ) ( 0 )
113: ///
114: /// where alpha and beta are scalars, and x is an (n-1)-element real
115: /// vector. H is represented in the form
116: ///
117: /// H = I - tau * ( 1 ) * ( 1 v' ) ,
118: /// ( v )
119: ///
120: /// where tau is a real scalar and v is a real (n-1)-element
121: /// vector.
122: ///
123: /// If the elements of x are all zero, then tau = 0 and H is taken to be
124: /// the unit matrix.
125: ///
126: /// Otherwise 1 .LE. tau .LE. 2.
127: ///
128: ///</summary>
129: /// <param name="N">
130: /// (input) INTEGER
131: /// The order of the elementary reflector.
132: ///</param>
133: /// <param name="ALPHA">
134: /// (input/output) DOUBLE PRECISION
135: /// On entry, the value alpha.
136: /// On exit, it is overwritten with the value beta.
137: ///</param>
138: /// <param name="X">
139: /// (input/output) DOUBLE PRECISION array, dimension
140: /// (1+(N-2)*abs(INCX))
141: /// On entry, the vector x.
142: /// On exit, it is overwritten with the vector v.
143: ///</param>
144: /// <param name="INCX">
145: /// (input) INTEGER
146: /// The increment between elements of X. INCX .GT. 0.
147: ///</param>
148: /// <param name="TAU">
149: /// (output) DOUBLE PRECISION
150: /// The value tau.
151: ///</param>
152: public void Run(int N, ref double ALPHA, ref double[] X, int offset_x, int INCX, ref double TAU)
153: {
154:
155: #region Array Index Correction
156:
157: int o_x = -1 + offset_x;
158:
159: #endregion
160:
161:
162: #region Prolog
163:
164: // *
165: // * -- LAPACK auxiliary routine (version 3.1) --
166: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
167: // * November 2006
168: // *
169: // * .. Scalar Arguments ..
170: // * ..
171: // * .. Array Arguments ..
172: // * ..
173: // *
174: // * Purpose
175: // * =======
176: // *
177: // * DLARFG generates a real elementary reflector H of order n, such
178: // * that
179: // *
180: // * H * ( alpha ) = ( beta ), H' * H = I.
181: // * ( x ) ( 0 )
182: // *
183: // * where alpha and beta are scalars, and x is an (n-1)-element real
184: // * vector. H is represented in the form
185: // *
186: // * H = I - tau * ( 1 ) * ( 1 v' ) ,
187: // * ( v )
188: // *
189: // * where tau is a real scalar and v is a real (n-1)-element
190: // * vector.
191: // *
192: // * If the elements of x are all zero, then tau = 0 and H is taken to be
193: // * the unit matrix.
194: // *
195: // * Otherwise 1 <= tau <= 2.
196: // *
197: // * Arguments
198: // * =========
199: // *
200: // * N (input) INTEGER
201: // * The order of the elementary reflector.
202: // *
203: // * ALPHA (input/output) DOUBLE PRECISION
204: // * On entry, the value alpha.
205: // * On exit, it is overwritten with the value beta.
206: // *
207: // * X (input/output) DOUBLE PRECISION array, dimension
208: // * (1+(N-2)*abs(INCX))
209: // * On entry, the vector x.
210: // * On exit, it is overwritten with the vector v.
211: // *
212: // * INCX (input) INTEGER
213: // * The increment between elements of X. INCX > 0.
214: // *
215: // * TAU (output) DOUBLE PRECISION
216: // * The value tau.
217: // *
218: // * =====================================================================
219: // *
220: // * .. Parameters ..
221: // * ..
222: // * .. Local Scalars ..
223: // * ..
224: // * .. External Functions ..
225: // * ..
226: // * .. Intrinsic Functions ..
227: // INTRINSIC ABS, SIGN;
228: // * ..
229: // * .. External Subroutines ..
230: // * ..
231: // * .. Executable Statements ..
232: // *
233:
234: #endregion
235:
236:
237: #region Body
238:
239: if (N <= 1)
240: {
241: TAU = ZERO;
242: return;
243: }
244: // *
245: XNORM = this._dnrm2.Run(N - 1, X, offset_x, INCX);
246: // *
247: if (XNORM == ZERO)
248: {
249: // *
250: // * H = I
251: // *
252: TAU = ZERO;
253: }
254: else
255: {
256: // *
257: // * general case
258: // *
259: BETA = - FortranLib.Sign(this._dlapy2.Run(ALPHA, XNORM),ALPHA);
260: SAFMIN = this._dlamch.Run("S") / this._dlamch.Run("E");
261: if (Math.Abs(BETA) < SAFMIN)
262: {
263: // *
264: // * XNORM, BETA may be inaccurate; scale X and recompute them
265: // *
266: RSAFMN = ONE / SAFMIN;
267: KNT = 0;
268: LABEL10:;
269: KNT = KNT + 1;
270: this._dscal.Run(N - 1, RSAFMN, ref X, offset_x, INCX);
271: BETA = BETA * RSAFMN;
272: ALPHA = ALPHA * RSAFMN;
273: if (Math.Abs(BETA) < SAFMIN) goto LABEL10;
274: // *
275: // * New BETA is at most 1, at least SAFMIN
276: // *
277: XNORM = this._dnrm2.Run(N - 1, X, offset_x, INCX);
278: BETA = - FortranLib.Sign(this._dlapy2.Run(ALPHA, XNORM),ALPHA);
279: TAU = (BETA - ALPHA) / BETA;
280: this._dscal.Run(N - 1, ONE / (ALPHA - BETA), ref X, offset_x, INCX);
281: // *
282: // * If ALPHA is subnormal, it may lose relative accuracy
283: // *
284: ALPHA = BETA;
285: for (J = 1; J <= KNT; J++)
286: {
287: ALPHA = ALPHA * SAFMIN;
288: }
289: }
290: else
291: {
292: TAU = (BETA - ALPHA) / BETA;
293: this._dscal.Run(N - 1, ONE / (ALPHA - BETA), ref X, offset_x, INCX);
294: ALPHA = BETA;
295: }
296: }
297: // *
298: return;
299: // *
300: // * End of DLARFG
301: // *
302:
303: #endregion
304:
305: }
306: }
307: }