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: ///</summary>
24: public class DLAQR1
25: {
26:
27:
28: #region Fields
29:
30: const double ZERO = 0.0E0; double H21S = 0; double H31S = 0; double S = 0;
31:
32: #endregion
33:
34: public DLAQR1()
35: {
36:
37: }
38:
39: /// <param name="N">
40: /// (input) integer
41: /// Order of the matrix H. N must be either 2 or 3.
42: ///</param>
43: /// <param name="H">
44: /// (input) DOUBLE PRECISION array of dimension (LDH,N)
45: /// The 2-by-2 or 3-by-3 matrix H in (*).
46: ///</param>
47: /// <param name="LDH">
48: /// (input) integer
49: /// The leading dimension of H as declared in
50: /// the calling procedure. LDH.GE.N
51: ///</param>
52: /// <param name="SR1">
53: /// (input) DOUBLE PRECISION
54: ///</param>
55: /// <param name="SI1">
56: /// The shifts in (*).
57: ///</param>
58: /// <param name="V">
59: /// (output) DOUBLE PRECISION array of dimension N
60: /// A scalar multiple of the first column of the
61: /// matrix K in (*).
62: ///</param>
63: public void Run(int N, double[] H, int offset_h, int LDH, double SR1, double SI1, double SR2
64: , double SI2, ref double[] V, int offset_v)
65: {
66:
67: #region Array Index Correction
68:
69: int o_h = -1 - LDH + offset_h; int o_v = -1 + offset_v;
70:
71: #endregion
72:
73:
74: #region Prolog
75:
76: // *
77: // * -- LAPACK auxiliary routine (version 3.1) --
78: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
79: // * November 2006
80: // *
81: // * .. Scalar Arguments ..
82: // * ..
83: // * .. Array Arguments ..
84: // * ..
85: // *
86: // * Given a 2-by-2 or 3-by-3 matrix H, DLAQR1 sets v to a
87: // * scalar multiple of the first column of the product
88: // *
89: // * (*) K = (H - (sr1 + i*si1)*I)*(H - (sr2 + i*si2)*I)
90: // *
91: // * scaling to avoid overflows and most underflows. It
92: // * is assumed that either
93: // *
94: // * 1) sr1 = sr2 and si1 = -si2
95: // * or
96: // * 2) si1 = si2 = 0.
97: // *
98: // * This is useful for starting double implicit shift bulges
99: // * in the QR algorithm.
100: // *
101: // *
102: // * N (input) integer
103: // * Order of the matrix H. N must be either 2 or 3.
104: // *
105: // * H (input) DOUBLE PRECISION array of dimension (LDH,N)
106: // * The 2-by-2 or 3-by-3 matrix H in (*).
107: // *
108: // * LDH (input) integer
109: // * The leading dimension of H as declared in
110: // * the calling procedure. LDH.GE.N
111: // *
112: // * SR1 (input) DOUBLE PRECISION
113: // * SI1 The shifts in (*).
114: // * SR2
115: // * SI2
116: // *
117: // * V (output) DOUBLE PRECISION array of dimension N
118: // * A scalar multiple of the first column of the
119: // * matrix K in (*).
120: // *
121: // * ================================================================
122: // * Based on contributions by
123: // * Karen Braman and Ralph Byers, Department of Mathematics,
124: // * University of Kansas, USA
125: // *
126: // * ================================================================
127: // *
128: // * .. Parameters ..
129: // * ..
130: // * .. Local Scalars ..
131: // * ..
132: // * .. Intrinsic Functions ..
133: // INTRINSIC ABS;
134: // * ..
135: // * .. Executable Statements ..
136:
137: #endregion
138:
139:
140: #region Body
141:
142: if (N == 2)
143: {
144: S = Math.Abs(H[1+1 * LDH + o_h] - SR2) + Math.Abs(SI2) + Math.Abs(H[2+1 * LDH + o_h]);
145: if (S == ZERO)
146: {
147: V[1 + o_v] = ZERO;
148: V[2 + o_v] = ZERO;
149: }
150: else
151: {
152: H21S = H[2+1 * LDH + o_h] / S;
153: V[1 + o_v] = H21S * H[1+2 * LDH + o_h] + (H[1+1 * LDH + o_h] - SR1) * ((H[1+1 * LDH + o_h] - SR2) / S) - SI1 * (SI2 / S);
154: V[2 + o_v] = H21S * (H[1+1 * LDH + o_h] + H[2+2 * LDH + o_h] - SR1 - SR2);
155: }
156: }
157: else
158: {
159: S = Math.Abs(H[1+1 * LDH + o_h] - SR2) + Math.Abs(SI2) + Math.Abs(H[2+1 * LDH + o_h]) + Math.Abs(H[3+1 * LDH + o_h]);
160: if (S == ZERO)
161: {
162: V[1 + o_v] = ZERO;
163: V[2 + o_v] = ZERO;
164: V[3 + o_v] = ZERO;
165: }
166: else
167: {
168: H21S = H[2+1 * LDH + o_h] / S;
169: H31S = H[3+1 * LDH + o_h] / S;
170: V[1 + o_v] = (H[1+1 * LDH + o_h] - SR1) * ((H[1+1 * LDH + o_h] - SR2) / S) - SI1 * (SI2 / S) + H[1+2 * LDH + o_h] * H21S + H[1+3 * LDH + o_h] * H31S;
171: V[2 + o_v] = H21S * (H[1+1 * LDH + o_h] + H[2+2 * LDH + o_h] - SR1 - SR2) + H[2+3 * LDH + o_h] * H31S;
172: V[3 + o_v] = H31S * (H[1+1 * LDH + o_h] + H[3+3 * LDH + o_h] - SR1 - SR2) + H21S * H[3+2 * LDH + o_h];
173: }
174: }
175:
176: #endregion
177:
178: }
179: }
180: }