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: /// DLASV2 computes the singular value decomposition of a 2-by-2
27: /// triangular matrix
28: /// [ F G ]
29: /// [ 0 H ].
30: /// On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
31: /// smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
32: /// right singular vectors for abs(SSMAX), giving the decomposition
33: ///
34: /// [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ]
35: /// [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ].
36: ///
37: ///</summary>
38: public class DLASV2
39: {
40:
41:
42: #region Dependencies
43:
44: DLAMCH _dlamch;
45:
46: #endregion
47:
48:
49: #region Fields
50:
51: const double ZERO = 0.0E0; const double HALF = 0.5E0; const double ONE = 1.0E0; const double TWO = 2.0E0;
52: const double FOUR = 4.0E0;bool GASMAL = false; bool SWAP = false; int PMAX = 0; double A = 0; double CLT = 0;
53: double CRT = 0;double D = 0; double FA = 0; double FT = 0; double GA = 0; double GT = 0; double HA = 0; double HT = 0;
54: double L = 0;double M = 0; double MM = 0; double R = 0; double S = 0; double SLT = 0; double SRT = 0; double T = 0;
55: double TEMP = 0;double TSIGN = 0; double TT = 0;
56:
57: #endregion
58:
59: public DLASV2(DLAMCH dlamch)
60: {
61:
62:
63: #region Set Dependencies
64:
65: this._dlamch = dlamch;
66:
67: #endregion
68:
69: }
70:
71: public DLASV2()
72: {
73:
74:
75: #region Dependencies (Initialization)
76:
77: LSAME lsame = new LSAME();
78: DLAMC3 dlamc3 = new DLAMC3();
79: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
80: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
81: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
82: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
83: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
84:
85: #endregion
86:
87:
88: #region Set Dependencies
89:
90: this._dlamch = dlamch;
91:
92: #endregion
93:
94: }
95: /// <summary>
96: /// Purpose
97: /// =======
98: ///
99: /// DLASV2 computes the singular value decomposition of a 2-by-2
100: /// triangular matrix
101: /// [ F G ]
102: /// [ 0 H ].
103: /// On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
104: /// smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
105: /// right singular vectors for abs(SSMAX), giving the decomposition
106: ///
107: /// [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ]
108: /// [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ].
109: ///
110: ///</summary>
111: /// <param name="F">
112: /// (input) DOUBLE PRECISION
113: /// The (1,1) element of the 2-by-2 matrix.
114: ///</param>
115: /// <param name="G">
116: /// (input) DOUBLE PRECISION
117: /// The (1,2) element of the 2-by-2 matrix.
118: ///</param>
119: /// <param name="H">
120: /// (input) DOUBLE PRECISION
121: /// The (2,2) element of the 2-by-2 matrix.
122: ///</param>
123: /// <param name="SSMIN">
124: /// (output) DOUBLE PRECISION
125: /// abs(SSMIN) is the smaller singular value.
126: ///</param>
127: /// <param name="SSMAX">
128: /// (output) DOUBLE PRECISION
129: /// abs(SSMAX) is the larger singular value.
130: ///</param>
131: /// <param name="SNR">
132: /// (output) DOUBLE PRECISION
133: ///</param>
134: /// <param name="CSR">
135: /// (output) DOUBLE PRECISION
136: /// The vector (CSR, SNR) is a unit right singular vector for the
137: /// singular value abs(SSMAX).
138: ///</param>
139: /// <param name="SNL">
140: /// (output) DOUBLE PRECISION
141: ///</param>
142: /// <param name="CSL">
143: /// (output) DOUBLE PRECISION
144: /// The vector (CSL, SNL) is a unit left singular vector for the
145: /// singular value abs(SSMAX).
146: ///</param>
147: public void Run(double F, double G, double H, ref double SSMIN, ref double SSMAX, ref double SNR
148: , ref double CSR, ref double SNL, ref double CSL)
149: {
150:
151: #region Prolog
152:
153: // *
154: // * -- LAPACK auxiliary routine (version 3.1) --
155: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
156: // * November 2006
157: // *
158: // * .. Scalar Arguments ..
159: // * ..
160: // *
161: // * Purpose
162: // * =======
163: // *
164: // * DLASV2 computes the singular value decomposition of a 2-by-2
165: // * triangular matrix
166: // * [ F G ]
167: // * [ 0 H ].
168: // * On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
169: // * smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
170: // * right singular vectors for abs(SSMAX), giving the decomposition
171: // *
172: // * [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ]
173: // * [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ].
174: // *
175: // * Arguments
176: // * =========
177: // *
178: // * F (input) DOUBLE PRECISION
179: // * The (1,1) element of the 2-by-2 matrix.
180: // *
181: // * G (input) DOUBLE PRECISION
182: // * The (1,2) element of the 2-by-2 matrix.
183: // *
184: // * H (input) DOUBLE PRECISION
185: // * The (2,2) element of the 2-by-2 matrix.
186: // *
187: // * SSMIN (output) DOUBLE PRECISION
188: // * abs(SSMIN) is the smaller singular value.
189: // *
190: // * SSMAX (output) DOUBLE PRECISION
191: // * abs(SSMAX) is the larger singular value.
192: // *
193: // * SNL (output) DOUBLE PRECISION
194: // * CSL (output) DOUBLE PRECISION
195: // * The vector (CSL, SNL) is a unit left singular vector for the
196: // * singular value abs(SSMAX).
197: // *
198: // * SNR (output) DOUBLE PRECISION
199: // * CSR (output) DOUBLE PRECISION
200: // * The vector (CSR, SNR) is a unit right singular vector for the
201: // * singular value abs(SSMAX).
202: // *
203: // * Further Details
204: // * ===============
205: // *
206: // * Any input parameter may be aliased with any output parameter.
207: // *
208: // * Barring over/underflow and assuming a guard digit in subtraction, all
209: // * output quantities are correct to within a few units in the last
210: // * place (ulps).
211: // *
212: // * In IEEE arithmetic, the code works correctly if one matrix element is
213: // * infinite.
214: // *
215: // * Overflow will not occur unless the largest singular value itself
216: // * overflows or is within a few ulps of overflow. (On machines with
217: // * partial overflow, like the Cray, overflow may occur if the largest
218: // * singular value is within a factor of 2 of overflow.)
219: // *
220: // * Underflow is harmless if underflow is gradual. Otherwise, results
221: // * may correspond to a matrix modified by perturbations of size near
222: // * the underflow threshold.
223: // *
224: // * =====================================================================
225: // *
226: // * .. Parameters ..
227: // * ..
228: // * .. Local Scalars ..
229: // * ..
230: // * .. Intrinsic Functions ..
231: // INTRINSIC ABS, SIGN, SQRT;
232: // * ..
233: // * .. External Functions ..
234: // * ..
235: // * .. Executable Statements ..
236: // *
237:
238: #endregion
239:
240:
241: #region Body
242:
243: FT = F;
244: FA = Math.Abs(FT);
245: HT = H;
246: HA = Math.Abs(H);
247: // *
248: // * PMAX points to the maximum absolute element of matrix
249: // * PMAX = 1 if F largest in absolute values
250: // * PMAX = 2 if G largest in absolute values
251: // * PMAX = 3 if H largest in absolute values
252: // *
253: PMAX = 1;
254: SWAP = (HA > FA);
255: if (SWAP)
256: {
257: PMAX = 3;
258: TEMP = FT;
259: FT = HT;
260: HT = TEMP;
261: TEMP = FA;
262: FA = HA;
263: HA = TEMP;
264: // *
265: // * Now FA .ge. HA
266: // *
267: }
268: GT = G;
269: GA = Math.Abs(GT);
270: if (GA == ZERO)
271: {
272: // *
273: // * Diagonal matrix
274: // *
275: SSMIN = HA;
276: SSMAX = FA;
277: CLT = ONE;
278: CRT = ONE;
279: SLT = ZERO;
280: SRT = ZERO;
281: }
282: else
283: {
284: GASMAL = true;
285: if (GA > FA)
286: {
287: PMAX = 2;
288: if ((FA / GA) < this._dlamch.Run("EPS"))
289: {
290: // *
291: // * Case of very large GA
292: // *
293: GASMAL = false;
294: SSMAX = GA;
295: if (HA > ONE)
296: {
297: SSMIN = FA / (GA / HA);
298: }
299: else
300: {
301: SSMIN = (FA / GA) * HA;
302: }
303: CLT = ONE;
304: SLT = HT / GT;
305: SRT = ONE;
306: CRT = FT / GT;
307: }
308: }
309: if (GASMAL)
310: {
311: // *
312: // * Normal case
313: // *
314: D = FA - HA;
315: if (D == FA)
316: {
317: // *
318: // * Copes with infinite F or H
319: // *
320: L = ONE;
321: }
322: else
323: {
324: L = D / FA;
325: }
326: // *
327: // * Note that 0 .le. L .le. 1
328: // *
329: M = GT / FT;
330: // *
331: // * Note that abs(M) .le. 1/macheps
332: // *
333: T = TWO - L;
334: // *
335: // * Note that T .ge. 1
336: // *
337: MM = M * M;
338: TT = T * T;
339: S = Math.Sqrt(TT + MM);
340: // *
341: // * Note that 1 .le. S .le. 1 + 1/macheps
342: // *
343: if (L == ZERO)
344: {
345: R = Math.Abs(M);
346: }
347: else
348: {
349: R = Math.Sqrt(L * L + MM);
350: }
351: // *
352: // * Note that 0 .le. R .le. 1 + 1/macheps
353: // *
354: A = HALF * (S + R);
355: // *
356: // * Note that 1 .le. A .le. 1 + abs(M)
357: // *
358: SSMIN = HA / A;
359: SSMAX = FA * A;
360: if (MM == ZERO)
361: {
362: // *
363: // * Note that M is very tiny
364: // *
365: if (L == ZERO)
366: {
367: T = FortranLib.Sign(TWO,FT) * FortranLib.Sign(ONE,GT);
368: }
369: else
370: {
371: T = GT / FortranLib.Sign(D,FT) + M / T;
372: }
373: }
374: else
375: {
376: T = (M / (S + T) + M / (R + L)) * (ONE + A);
377: }
378: L = Math.Sqrt(T * T + FOUR);
379: CRT = TWO / L;
380: SRT = T / L;
381: CLT = (CRT + SRT * M) / A;
382: SLT = (HT / FT) * SRT / A;
383: }
384: }
385: if (SWAP)
386: {
387: CSL = SRT;
388: SNL = CRT;
389: CSR = SLT;
390: SNR = CLT;
391: }
392: else
393: {
394: CSL = CLT;
395: SNL = SLT;
396: CSR = CRT;
397: SNR = SRT;
398: }
399: // *
400: // * Correct signs of SSMAX and SSMIN
401: // *
402: if (PMAX == 1) TSIGN = FortranLib.Sign(ONE,CSR) * FortranLib.Sign(ONE,CSL) * FortranLib.Sign(ONE,F);
403: if (PMAX == 2) TSIGN = FortranLib.Sign(ONE,SNR) * FortranLib.Sign(ONE,CSL) * FortranLib.Sign(ONE,G);
404: if (PMAX == 3) TSIGN = FortranLib.Sign(ONE,SNR) * FortranLib.Sign(ONE,SNL) * FortranLib.Sign(ONE,H);
405: SSMAX = FortranLib.Sign(SSMAX,TSIGN);
406: SSMIN = FortranLib.Sign(SSMIN,TSIGN * FortranLib.Sign(ONE,F) * FortranLib.Sign(ONE,H));
407: return;
408: // *
409: // * End of DLASV2
410: // *
411:
412: #endregion
413:
414: }
415: }
416: }