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: /// DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
27: /// that if ( UPPER ) then
28: ///
29: /// U'*A*Q = U'*( A1 A2 )*Q = ( x 0 )
30: /// ( 0 A3 ) ( x x )
31: /// and
32: /// V'*B*Q = V'*( B1 B2 )*Q = ( x 0 )
33: /// ( 0 B3 ) ( x x )
34: ///
35: /// or if ( .NOT.UPPER ) then
36: ///
37: /// U'*A*Q = U'*( A1 0 )*Q = ( x x )
38: /// ( A2 A3 ) ( 0 x )
39: /// and
40: /// V'*B*Q = V'*( B1 0 )*Q = ( x x )
41: /// ( B2 B3 ) ( 0 x )
42: ///
43: /// The rows of the transformed A and B are parallel, where
44: ///
45: /// U = ( CSU SNU ), V = ( CSV SNV ), Q = ( CSQ SNQ )
46: /// ( -SNU CSU ) ( -SNV CSV ) ( -SNQ CSQ )
47: ///
48: /// Z' denotes the transpose of Z.
49: ///</summary>
50: public class DLAGS2
51: {
52:
53:
54: #region Dependencies
55:
56: DLARTG _dlartg; DLASV2 _dlasv2;
57:
58: #endregion
59:
60:
61: #region Fields
62:
63: const double ZERO = 0.0E+0; double A = 0; double AUA11 = 0; double AUA12 = 0; double AUA21 = 0; double AUA22 = 0;
64: double AVB11 = 0;double AVB12 = 0; double AVB21 = 0; double AVB22 = 0; double B = 0; double C = 0; double CSL = 0;
65: double CSR = 0;double D = 0; double R = 0; double S1 = 0; double S2 = 0; double SNL = 0; double SNR = 0; double UA11 = 0;
66: double UA11R = 0;double UA12 = 0; double UA21 = 0; double UA22 = 0; double UA22R = 0; double VB11 = 0; double VB11R = 0;
67: double VB12 = 0;double VB21 = 0; double VB22 = 0; double VB22R = 0;
68:
69: #endregion
70:
71: public DLAGS2(DLARTG dlartg, DLASV2 dlasv2)
72: {
73:
74:
75: #region Set Dependencies
76:
77: this._dlartg = dlartg; this._dlasv2 = dlasv2;
78:
79: #endregion
80:
81: }
82:
83: public DLAGS2()
84: {
85:
86:
87: #region Dependencies (Initialization)
88:
89: LSAME lsame = new LSAME();
90: DLAMC3 dlamc3 = new DLAMC3();
91: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
92: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
93: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
94: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
95: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
96: DLARTG dlartg = new DLARTG(dlamch);
97: DLASV2 dlasv2 = new DLASV2(dlamch);
98:
99: #endregion
100:
101:
102: #region Set Dependencies
103:
104: this._dlartg = dlartg; this._dlasv2 = dlasv2;
105:
106: #endregion
107:
108: }
109: /// <summary>
110: /// Purpose
111: /// =======
112: ///
113: /// DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
114: /// that if ( UPPER ) then
115: ///
116: /// U'*A*Q = U'*( A1 A2 )*Q = ( x 0 )
117: /// ( 0 A3 ) ( x x )
118: /// and
119: /// V'*B*Q = V'*( B1 B2 )*Q = ( x 0 )
120: /// ( 0 B3 ) ( x x )
121: ///
122: /// or if ( .NOT.UPPER ) then
123: ///
124: /// U'*A*Q = U'*( A1 0 )*Q = ( x x )
125: /// ( A2 A3 ) ( 0 x )
126: /// and
127: /// V'*B*Q = V'*( B1 0 )*Q = ( x x )
128: /// ( B2 B3 ) ( 0 x )
129: ///
130: /// The rows of the transformed A and B are parallel, where
131: ///
132: /// U = ( CSU SNU ), V = ( CSV SNV ), Q = ( CSQ SNQ )
133: /// ( -SNU CSU ) ( -SNV CSV ) ( -SNQ CSQ )
134: ///
135: /// Z' denotes the transpose of Z.
136: ///</summary>
137: /// <param name="UPPER">
138: /// (input) LOGICAL
139: /// = .TRUE.: the input matrices A and B are upper triangular.
140: /// = .FALSE.: the input matrices A and B are lower triangular.
141: ///</param>
142: /// <param name="A1">
143: /// (input) DOUBLE PRECISION
144: ///</param>
145: /// <param name="A2">
146: /// (input) DOUBLE PRECISION
147: ///</param>
148: /// <param name="A3">
149: /// (input) DOUBLE PRECISION
150: /// On entry, A1, A2 and A3 are elements of the input 2-by-2
151: /// upper (lower) triangular matrix A.
152: ///</param>
153: /// <param name="B1">
154: /// (input) DOUBLE PRECISION
155: ///</param>
156: /// <param name="B2">
157: /// (input) DOUBLE PRECISION
158: ///</param>
159: /// <param name="B3">
160: /// (input) DOUBLE PRECISION
161: /// On entry, B1, B2 and B3 are elements of the input 2-by-2
162: /// upper (lower) triangular matrix B.
163: ///</param>
164: /// <param name="CSU">
165: /// (output) DOUBLE PRECISION
166: ///</param>
167: /// <param name="SNU">
168: /// (output) DOUBLE PRECISION
169: /// The desired orthogonal matrix U.
170: ///</param>
171: /// <param name="CSV">
172: /// (output) DOUBLE PRECISION
173: ///</param>
174: /// <param name="SNV">
175: /// (output) DOUBLE PRECISION
176: /// The desired orthogonal matrix V.
177: ///</param>
178: /// <param name="CSQ">
179: /// (output) DOUBLE PRECISION
180: ///</param>
181: /// <param name="SNQ">
182: /// (output) DOUBLE PRECISION
183: /// The desired orthogonal matrix Q.
184: ///</param>
185: public void Run(bool UPPER, double A1, double A2, double A3, double B1, double B2
186: , double B3, ref double CSU, ref double SNU, ref double CSV, ref double SNV, ref double CSQ
187: , ref double SNQ)
188: {
189:
190: #region Prolog
191:
192: // *
193: // * -- LAPACK auxiliary routine (version 3.1) --
194: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
195: // * November 2006
196: // *
197: // * .. Scalar Arguments ..
198: // * ..
199: // *
200: // * Purpose
201: // * =======
202: // *
203: // * DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
204: // * that if ( UPPER ) then
205: // *
206: // * U'*A*Q = U'*( A1 A2 )*Q = ( x 0 )
207: // * ( 0 A3 ) ( x x )
208: // * and
209: // * V'*B*Q = V'*( B1 B2 )*Q = ( x 0 )
210: // * ( 0 B3 ) ( x x )
211: // *
212: // * or if ( .NOT.UPPER ) then
213: // *
214: // * U'*A*Q = U'*( A1 0 )*Q = ( x x )
215: // * ( A2 A3 ) ( 0 x )
216: // * and
217: // * V'*B*Q = V'*( B1 0 )*Q = ( x x )
218: // * ( B2 B3 ) ( 0 x )
219: // *
220: // * The rows of the transformed A and B are parallel, where
221: // *
222: // * U = ( CSU SNU ), V = ( CSV SNV ), Q = ( CSQ SNQ )
223: // * ( -SNU CSU ) ( -SNV CSV ) ( -SNQ CSQ )
224: // *
225: // * Z' denotes the transpose of Z.
226: // *
227: // *
228: // * Arguments
229: // * =========
230: // *
231: // * UPPER (input) LOGICAL
232: // * = .TRUE.: the input matrices A and B are upper triangular.
233: // * = .FALSE.: the input matrices A and B are lower triangular.
234: // *
235: // * A1 (input) DOUBLE PRECISION
236: // * A2 (input) DOUBLE PRECISION
237: // * A3 (input) DOUBLE PRECISION
238: // * On entry, A1, A2 and A3 are elements of the input 2-by-2
239: // * upper (lower) triangular matrix A.
240: // *
241: // * B1 (input) DOUBLE PRECISION
242: // * B2 (input) DOUBLE PRECISION
243: // * B3 (input) DOUBLE PRECISION
244: // * On entry, B1, B2 and B3 are elements of the input 2-by-2
245: // * upper (lower) triangular matrix B.
246: // *
247: // * CSU (output) DOUBLE PRECISION
248: // * SNU (output) DOUBLE PRECISION
249: // * The desired orthogonal matrix U.
250: // *
251: // * CSV (output) DOUBLE PRECISION
252: // * SNV (output) DOUBLE PRECISION
253: // * The desired orthogonal matrix V.
254: // *
255: // * CSQ (output) DOUBLE PRECISION
256: // * SNQ (output) DOUBLE PRECISION
257: // * The desired orthogonal matrix Q.
258: // *
259: // * =====================================================================
260: // *
261: // * .. Parameters ..
262: // * ..
263: // * .. Local Scalars ..
264: // * ..
265: // * .. External Subroutines ..
266: // * ..
267: // * .. Intrinsic Functions ..
268: // INTRINSIC ABS;
269: // * ..
270: // * .. Executable Statements ..
271: // *
272:
273: #endregion
274:
275:
276: #region Body
277:
278: if (UPPER)
279: {
280: // *
281: // * Input matrices A and B are upper triangular matrices
282: // *
283: // * Form matrix C = A*adj(B) = ( a b )
284: // * ( 0 d )
285: // *
286: A = A1 * B3;
287: D = A3 * B1;
288: B = A2 * B1 - A1 * B2;
289: // *
290: // * The SVD of real 2-by-2 triangular C
291: // *
292: // * ( CSL -SNL )*( A B )*( CSR SNR ) = ( R 0 )
293: // * ( SNL CSL ) ( 0 D ) ( -SNR CSR ) ( 0 T )
294: // *
295: this._dlasv2.Run(A, B, D, ref S1, ref S2, ref SNR
296: , ref CSR, ref SNL, ref CSL);
297: // *
298: if (Math.Abs(CSL) >= Math.Abs(SNL) || Math.Abs(CSR) >= Math.Abs(SNR))
299: {
300: // *
301: // * Compute the (1,1) and (1,2) elements of U'*A and V'*B,
302: // * and (1,2) element of |U|'*|A| and |V|'*|B|.
303: // *
304: UA11R = CSL * A1;
305: UA12 = CSL * A2 + SNL * A3;
306: // *
307: VB11R = CSR * B1;
308: VB12 = CSR * B2 + SNR * B3;
309: // *
310: AUA12 = Math.Abs(CSL) * Math.Abs(A2) + Math.Abs(SNL) * Math.Abs(A3);
311: AVB12 = Math.Abs(CSR) * Math.Abs(B2) + Math.Abs(SNR) * Math.Abs(B3);
312: // *
313: // * zero (1,2) elements of U'*A and V'*B
314: // *
315: if ((Math.Abs(UA11R) + Math.Abs(UA12)) != ZERO)
316: {
317: if (AUA12 / (Math.Abs(UA11R) + Math.Abs(UA12)) <= AVB12 / (Math.Abs(VB11R) + Math.Abs(VB12)))
318: {
319: this._dlartg.Run( - UA11R, UA12, ref CSQ, ref SNQ, ref R);
320: }
321: else
322: {
323: this._dlartg.Run( - VB11R, VB12, ref CSQ, ref SNQ, ref R);
324: }
325: }
326: else
327: {
328: this._dlartg.Run( - VB11R, VB12, ref CSQ, ref SNQ, ref R);
329: }
330: // *
331: CSU = CSL;
332: SNU = - SNL;
333: CSV = CSR;
334: SNV = - SNR;
335: // *
336: }
337: else
338: {
339: // *
340: // * Compute the (2,1) and (2,2) elements of U'*A and V'*B,
341: // * and (2,2) element of |U|'*|A| and |V|'*|B|.
342: // *
343: UA21 = - SNL * A1;
344: UA22 = - SNL * A2 + CSL * A3;
345: // *
346: VB21 = - SNR * B1;
347: VB22 = - SNR * B2 + CSR * B3;
348: // *
349: AUA22 = Math.Abs(SNL) * Math.Abs(A2) + Math.Abs(CSL) * Math.Abs(A3);
350: AVB22 = Math.Abs(SNR) * Math.Abs(B2) + Math.Abs(CSR) * Math.Abs(B3);
351: // *
352: // * zero (2,2) elements of U'*A and V'*B, and then swap.
353: // *
354: if ((Math.Abs(UA21) + Math.Abs(UA22)) != ZERO)
355: {
356: if (AUA22 / (Math.Abs(UA21) + Math.Abs(UA22)) <= AVB22 / (Math.Abs(VB21) + Math.Abs(VB22)))
357: {
358: this._dlartg.Run( - UA21, UA22, ref CSQ, ref SNQ, ref R);
359: }
360: else
361: {
362: this._dlartg.Run( - VB21, VB22, ref CSQ, ref SNQ, ref R);
363: }
364: }
365: else
366: {
367: this._dlartg.Run( - VB21, VB22, ref CSQ, ref SNQ, ref R);
368: }
369: // *
370: CSU = SNL;
371: SNU = CSL;
372: CSV = SNR;
373: SNV = CSR;
374: // *
375: }
376: // *
377: }
378: else
379: {
380: // *
381: // * Input matrices A and B are lower triangular matrices
382: // *
383: // * Form matrix C = A*adj(B) = ( a 0 )
384: // * ( c d )
385: // *
386: A = A1 * B3;
387: D = A3 * B1;
388: C = A2 * B3 - A3 * B2;
389: // *
390: // * The SVD of real 2-by-2 triangular C
391: // *
392: // * ( CSL -SNL )*( A 0 )*( CSR SNR ) = ( R 0 )
393: // * ( SNL CSL ) ( C D ) ( -SNR CSR ) ( 0 T )
394: // *
395: this._dlasv2.Run(A, C, D, ref S1, ref S2, ref SNR
396: , ref CSR, ref SNL, ref CSL);
397: // *
398: if (Math.Abs(CSR) >= Math.Abs(SNR) || Math.Abs(CSL) >= Math.Abs(SNL))
399: {
400: // *
401: // * Compute the (2,1) and (2,2) elements of U'*A and V'*B,
402: // * and (2,1) element of |U|'*|A| and |V|'*|B|.
403: // *
404: UA21 = - SNR * A1 + CSR * A2;
405: UA22R = CSR * A3;
406: // *
407: VB21 = - SNL * B1 + CSL * B2;
408: VB22R = CSL * B3;
409: // *
410: AUA21 = Math.Abs(SNR) * Math.Abs(A1) + Math.Abs(CSR) * Math.Abs(A2);
411: AVB21 = Math.Abs(SNL) * Math.Abs(B1) + Math.Abs(CSL) * Math.Abs(B2);
412: // *
413: // * zero (2,1) elements of U'*A and V'*B.
414: // *
415: if ((Math.Abs(UA21) + Math.Abs(UA22R)) != ZERO)
416: {
417: if (AUA21 / (Math.Abs(UA21) + Math.Abs(UA22R)) <= AVB21 / (Math.Abs(VB21) + Math.Abs(VB22R)))
418: {
419: this._dlartg.Run(UA22R, UA21, ref CSQ, ref SNQ, ref R);
420: }
421: else
422: {
423: this._dlartg.Run(VB22R, VB21, ref CSQ, ref SNQ, ref R);
424: }
425: }
426: else
427: {
428: this._dlartg.Run(VB22R, VB21, ref CSQ, ref SNQ, ref R);
429: }
430: // *
431: CSU = CSR;
432: SNU = - SNR;
433: CSV = CSL;
434: SNV = - SNL;
435: // *
436: }
437: else
438: {
439: // *
440: // * Compute the (1,1) and (1,2) elements of U'*A and V'*B,
441: // * and (1,1) element of |U|'*|A| and |V|'*|B|.
442: // *
443: UA11 = CSR * A1 + SNR * A2;
444: UA12 = SNR * A3;
445: // *
446: VB11 = CSL * B1 + SNL * B2;
447: VB12 = SNL * B3;
448: // *
449: AUA11 = Math.Abs(CSR) * Math.Abs(A1) + Math.Abs(SNR) * Math.Abs(A2);
450: AVB11 = Math.Abs(CSL) * Math.Abs(B1) + Math.Abs(SNL) * Math.Abs(B2);
451: // *
452: // * zero (1,1) elements of U'*A and V'*B, and then swap.
453: // *
454: if ((Math.Abs(UA11) + Math.Abs(UA12)) != ZERO)
455: {
456: if (AUA11 / (Math.Abs(UA11) + Math.Abs(UA12)) <= AVB11 / (Math.Abs(VB11) + Math.Abs(VB12)))
457: {
458: this._dlartg.Run(UA12, UA11, ref CSQ, ref SNQ, ref R);
459: }
460: else
461: {
462: this._dlartg.Run(VB12, VB11, ref CSQ, ref SNQ, ref R);
463: }
464: }
465: else
466: {
467: this._dlartg.Run(VB12, VB11, ref CSQ, ref SNQ, ref R);
468: }
469: // *
470: CSU = SNR;
471: SNU = CSR;
472: CSV = SNL;
473: SNV = CSL;
474: // *
475: }
476: // *
477: }
478: // *
479: return;
480: // *
481: // * End of DLAGS2
482: // *
483:
484: #endregion
485:
486: }
487: }
488: }