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 routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DLASQ1 computes the singular values of a real N-by-N bidiagonal
27: /// matrix with diagonal D and off-diagonal E. The singular values
28: /// are computed to high relative accuracy, in the absence of
29: /// denormalization, underflow and overflow. The algorithm was first
30: /// presented in
31: ///
32: /// "Accurate singular values and differential qd algorithms" by K. V.
33: /// Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
34: /// 1994,
35: ///
36: /// and the present implementation is described in "An implementation of
37: /// the dqds Algorithm (Positive Case)", LAPACK Working Note.
38: ///
39: ///</summary>
40: public class DLASQ1
41: {
42:
43:
44: #region Dependencies
45:
46: DCOPY _dcopy; DLAS2 _dlas2; DLASCL _dlascl; DLASQ2 _dlasq2; DLASRT _dlasrt; XERBLA _xerbla; DLAMCH _dlamch;
47:
48: #endregion
49:
50:
51: #region Fields
52:
53: const double ZERO = 0.0E0; int I = 0; int IINFO = 0; double EPS = 0; double SCALE = 0; double SAFMIN = 0;
54: double SIGMN = 0;double SIGMX = 0;
55:
56: #endregion
57:
58: public DLASQ1(DCOPY dcopy, DLAS2 dlas2, DLASCL dlascl, DLASQ2 dlasq2, DLASRT dlasrt, XERBLA xerbla, DLAMCH dlamch)
59: {
60:
61:
62: #region Set Dependencies
63:
64: this._dcopy = dcopy; this._dlas2 = dlas2; this._dlascl = dlascl; this._dlasq2 = dlasq2; this._dlasrt = dlasrt;
65: this._xerbla = xerbla;this._dlamch = dlamch;
66:
67: #endregion
68:
69: }
70:
71: public DLASQ1()
72: {
73:
74:
75: #region Dependencies (Initialization)
76:
77: DCOPY dcopy = new DCOPY();
78: DLAS2 dlas2 = new DLAS2();
79: LSAME lsame = new LSAME();
80: DLAMC3 dlamc3 = new DLAMC3();
81: XERBLA xerbla = new XERBLA();
82: DLASQ5 dlasq5 = new DLASQ5();
83: DLAZQ4 dlazq4 = new DLAZQ4();
84: IEEECK ieeeck = new IEEECK();
85: IPARMQ iparmq = new IPARMQ();
86: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
87: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
88: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
89: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
90: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
91: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
92: DLASQ6 dlasq6 = new DLASQ6(dlamch);
93: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
94: DLASRT dlasrt = new DLASRT(lsame, xerbla);
95: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
96: DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
97:
98: #endregion
99:
100:
101: #region Set Dependencies
102:
103: this._dcopy = dcopy; this._dlas2 = dlas2; this._dlascl = dlascl; this._dlasq2 = dlasq2; this._dlasrt = dlasrt;
104: this._xerbla = xerbla;this._dlamch = dlamch;
105:
106: #endregion
107:
108: }
109: /// <summary>
110: /// Purpose
111: /// =======
112: ///
113: /// DLASQ1 computes the singular values of a real N-by-N bidiagonal
114: /// matrix with diagonal D and off-diagonal E. The singular values
115: /// are computed to high relative accuracy, in the absence of
116: /// denormalization, underflow and overflow. The algorithm was first
117: /// presented in
118: ///
119: /// "Accurate singular values and differential qd algorithms" by K. V.
120: /// Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
121: /// 1994,
122: ///
123: /// and the present implementation is described in "An implementation of
124: /// the dqds Algorithm (Positive Case)", LAPACK Working Note.
125: ///
126: ///</summary>
127: /// <param name="N">
128: /// (input) INTEGER
129: /// The number of rows and columns in the matrix. N .GE. 0.
130: ///</param>
131: /// <param name="D">
132: /// (input/output) DOUBLE PRECISION array, dimension (N)
133: /// On entry, D contains the diagonal elements of the
134: /// bidiagonal matrix whose SVD is desired. On normal exit,
135: /// D contains the singular values in decreasing order.
136: ///</param>
137: /// <param name="E">
138: /// (input/output) DOUBLE PRECISION array, dimension (N)
139: /// On entry, elements E(1:N-1) contain the off-diagonal elements
140: /// of the bidiagonal matrix whose SVD is desired.
141: /// On exit, E is overwritten.
142: ///</param>
143: /// <param name="WORK">
144: /// (workspace) DOUBLE PRECISION array, dimension (4*N)
145: ///</param>
146: /// <param name="INFO">
147: /// (output) INTEGER
148: /// = 0: successful exit
149: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
150: /// .GT. 0: the algorithm failed
151: /// = 1, a split was marked by a positive value in E
152: /// = 2, current block of Z not diagonalized after 30*N
153: /// iterations (in inner while loop)
154: /// = 3, termination criterion of outer while loop not met
155: /// (program created more than N unreduced blocks)
156: ///</param>
157: public void Run(int N, ref double[] D, int offset_d, double[] E, int offset_e, ref double[] WORK, int offset_work, ref int INFO)
158: {
159:
160: #region Array Index Correction
161:
162: int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_work = -1 + offset_work;
163:
164: #endregion
165:
166:
167: #region Prolog
168:
169: // *
170: // * -- LAPACK routine (version 3.1) --
171: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
172: // * November 2006
173: // *
174: // * .. Scalar Arguments ..
175: // * ..
176: // * .. Array Arguments ..
177: // * ..
178: // *
179: // * Purpose
180: // * =======
181: // *
182: // * DLASQ1 computes the singular values of a real N-by-N bidiagonal
183: // * matrix with diagonal D and off-diagonal E. The singular values
184: // * are computed to high relative accuracy, in the absence of
185: // * denormalization, underflow and overflow. The algorithm was first
186: // * presented in
187: // *
188: // * "Accurate singular values and differential qd algorithms" by K. V.
189: // * Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
190: // * 1994,
191: // *
192: // * and the present implementation is described in "An implementation of
193: // * the dqds Algorithm (Positive Case)", LAPACK Working Note.
194: // *
195: // * Arguments
196: // * =========
197: // *
198: // * N (input) INTEGER
199: // * The number of rows and columns in the matrix. N >= 0.
200: // *
201: // * D (input/output) DOUBLE PRECISION array, dimension (N)
202: // * On entry, D contains the diagonal elements of the
203: // * bidiagonal matrix whose SVD is desired. On normal exit,
204: // * D contains the singular values in decreasing order.
205: // *
206: // * E (input/output) DOUBLE PRECISION array, dimension (N)
207: // * On entry, elements E(1:N-1) contain the off-diagonal elements
208: // * of the bidiagonal matrix whose SVD is desired.
209: // * On exit, E is overwritten.
210: // *
211: // * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
212: // *
213: // * INFO (output) INTEGER
214: // * = 0: successful exit
215: // * < 0: if INFO = -i, the i-th argument had an illegal value
216: // * > 0: the algorithm failed
217: // * = 1, a split was marked by a positive value in E
218: // * = 2, current block of Z not diagonalized after 30*N
219: // * iterations (in inner while loop)
220: // * = 3, termination criterion of outer while loop not met
221: // * (program created more than N unreduced blocks)
222: // *
223: // * =====================================================================
224: // *
225: // * .. Parameters ..
226: // * ..
227: // * .. Local Scalars ..
228: // * ..
229: // * .. External Subroutines ..
230: // * ..
231: // * .. External Functions ..
232: // * ..
233: // * .. Intrinsic Functions ..
234: // INTRINSIC ABS, MAX, SQRT;
235: // * ..
236: // * .. Executable Statements ..
237: // *
238:
239: #endregion
240:
241:
242: #region Body
243:
244: INFO = 0;
245: if (N < 0)
246: {
247: INFO = - 2;
248: this._xerbla.Run("DLASQ1", - INFO);
249: return;
250: }
251: else
252: {
253: if (N == 0)
254: {
255: return;
256: }
257: else
258: {
259: if (N == 1)
260: {
261: D[1 + o_d] = Math.Abs(D[1 + o_d]);
262: return;
263: }
264: else
265: {
266: if (N == 2)
267: {
268: this._dlas2.Run(D[1 + o_d], E[1 + o_e], D[2 + o_d], ref SIGMN, ref SIGMX);
269: D[1 + o_d] = SIGMX;
270: D[2 + o_d] = SIGMN;
271: return;
272: }
273: }
274: }
275: }
276: // *
277: // * Estimate the largest singular value.
278: // *
279: SIGMX = ZERO;
280: for (I = 1; I <= N - 1; I++)
281: {
282: D[I + o_d] = Math.Abs(D[I + o_d]);
283: SIGMX = Math.Max(SIGMX, Math.Abs(E[I + o_e]));
284: }
285: D[N + o_d] = Math.Abs(D[N + o_d]);
286: // *
287: // * Early return if SIGMX is zero (matrix is already diagonal).
288: // *
289: if (SIGMX == ZERO)
290: {
291: this._dlasrt.Run("D", N, ref D, offset_d, ref IINFO);
292: return;
293: }
294: // *
295: for (I = 1; I <= N; I++)
296: {
297: SIGMX = Math.Max(SIGMX, D[I + o_d]);
298: }
299: // *
300: // * Copy D and E into WORK (in the Z format) and scale (squaring the
301: // * input data makes scaling by a power of the radix pointless).
302: // *
303: EPS = this._dlamch.Run("Precision");
304: SAFMIN = this._dlamch.Run("Safe minimum");
305: SCALE = Math.Sqrt(EPS / SAFMIN);
306: this._dcopy.Run(N, D, offset_d, 1, ref WORK, 1 + o_work, 2);
307: this._dcopy.Run(N - 1, E, offset_e, 1, ref WORK, 2 + o_work, 2);
308: this._dlascl.Run("G", 0, 0, SIGMX, SCALE, 2 * N - 1
309: , 1, ref WORK, offset_work, 2 * N - 1, ref IINFO);
310: // *
311: // * Compute the q's and e's.
312: // *
313: for (I = 1; I <= 2 * N - 1; I++)
314: {
315: WORK[I + o_work] = Math.Pow(WORK[I + o_work],2);
316: }
317: WORK[2 * N + o_work] = ZERO;
318: // *
319: this._dlasq2.Run(N, ref WORK, offset_work, ref INFO);
320: // *
321: if (INFO == 0)
322: {
323: for (I = 1; I <= N; I++)
324: {
325: D[I + o_d] = Math.Sqrt(WORK[I + o_work]);
326: }
327: this._dlascl.Run("G", 0, 0, SCALE, SIGMX, N
328: , 1, ref D, offset_d, N, ref IINFO);
329: }
330: // *
331: return;
332: // *
333: // * End of DLASQ1
334: // *
335:
336: #endregion
337:
338: }
339: }
340: }