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: /// Using a divide and conquer approach, DLASD0 computes the singular
27: /// value decomposition (SVD) of a real upper bidiagonal N-by-M
28: /// matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
29: /// The algorithm computes orthogonal matrices U and VT such that
30: /// B = U * S * VT. The singular values S are overwritten on D.
31: ///
32: /// A related subroutine, DLASDA, computes only the singular values,
33: /// and optionally, the singular vectors in compact form.
34: ///
35: ///</summary>
36: public class DLASD0
37: {
38:
39:
40: #region Dependencies
41:
42: DLASD1 _dlasd1; DLASDQ _dlasdq; DLASDT _dlasdt; XERBLA _xerbla;
43:
44: #endregion
45:
46:
47: #region Fields
48:
49: int I = 0; int I1 = 0; int IC = 0; int IDXQ = 0; int IDXQC = 0; int IM1 = 0; int INODE = 0; int ITEMP = 0; int IWK = 0;
50: int J = 0;int LF = 0; int LL = 0; int LVL = 0; int M = 0; int NCC = 0; int ND = 0; int NDB1 = 0; int NDIML = 0;
51: int NDIMR = 0;int NL = 0; int NLF = 0; int NLP1 = 0; int NLVL = 0; int NR = 0; int NRF = 0; int NRP1 = 0; int SQREI = 0;
52: double ALPHA = 0;double BETA = 0;
53:
54: #endregion
55:
56: public DLASD0(DLASD1 dlasd1, DLASDQ dlasdq, DLASDT dlasdt, XERBLA xerbla)
57: {
58:
59:
60: #region Set Dependencies
61:
62: this._dlasd1 = dlasd1; this._dlasdq = dlasdq; this._dlasdt = dlasdt; this._xerbla = xerbla;
63:
64: #endregion
65:
66: }
67:
68: public DLASD0()
69: {
70:
71:
72: #region Dependencies (Initialization)
73:
74: DLAMRG dlamrg = new DLAMRG();
75: LSAME lsame = new LSAME();
76: DLAMC3 dlamc3 = new DLAMC3();
77: XERBLA xerbla = new XERBLA();
78: DLAPY2 dlapy2 = new DLAPY2();
79: DCOPY dcopy = new DCOPY();
80: DROT drot = new DROT();
81: DNRM2 dnrm2 = new DNRM2();
82: DLASD5 dlasd5 = new DLASD5();
83: DLAS2 dlas2 = new DLAS2();
84: DLASQ5 dlasq5 = new DLASQ5();
85: DLAZQ4 dlazq4 = new DLAZQ4();
86: IEEECK ieeeck = new IEEECK();
87: IPARMQ iparmq = new IPARMQ();
88: DSCAL dscal = new DSCAL();
89: DSWAP dswap = new DSWAP();
90: DLASDT dlasdt = new DLASDT();
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: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
97: DLACPY dlacpy = new DLACPY(lsame);
98: DLASET dlaset = new DLASET(lsame);
99: DLASD2 dlasd2 = new DLASD2(dlamch, dlapy2, dcopy, dlacpy, dlamrg, dlaset, drot, xerbla);
100: DGEMM dgemm = new DGEMM(lsame, xerbla);
101: DLAED6 dlaed6 = new DLAED6(dlamch);
102: DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
103: DLASD3 dlasd3 = new DLASD3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlascl, dlasd4, xerbla);
104: DLASD1 dlasd1 = new DLASD1(dlamrg, dlascl, dlasd2, dlasd3, xerbla);
105: DLARTG dlartg = new DLARTG(dlamch);
106: DLASQ6 dlasq6 = new DLASQ6(dlamch);
107: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
108: DLASRT dlasrt = new DLASRT(lsame, xerbla);
109: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
110: DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
111: DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
112: DLASR dlasr = new DLASR(lsame, xerbla);
113: DLASV2 dlasv2 = new DLASV2(dlamch);
114: DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap
115: , xerbla);
116: DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);
117:
118: #endregion
119:
120:
121: #region Set Dependencies
122:
123: this._dlasd1 = dlasd1; this._dlasdq = dlasdq; this._dlasdt = dlasdt; this._xerbla = xerbla;
124:
125: #endregion
126:
127: }
128: /// <summary>
129: /// Purpose
130: /// =======
131: ///
132: /// Using a divide and conquer approach, DLASD0 computes the singular
133: /// value decomposition (SVD) of a real upper bidiagonal N-by-M
134: /// matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
135: /// The algorithm computes orthogonal matrices U and VT such that
136: /// B = U * S * VT. The singular values S are overwritten on D.
137: ///
138: /// A related subroutine, DLASDA, computes only the singular values,
139: /// and optionally, the singular vectors in compact form.
140: ///
141: ///</summary>
142: /// <param name="N">
143: /// (input) INTEGER
144: /// On entry, the row dimension of the upper bidiagonal matrix.
145: /// This is also the dimension of the main diagonal array D.
146: ///</param>
147: /// <param name="SQRE">
148: /// (input) INTEGER
149: /// Specifies the column dimension of the bidiagonal matrix.
150: /// = 0: The bidiagonal matrix has column dimension M = N;
151: /// = 1: The bidiagonal matrix has column dimension M = N+1;
152: ///</param>
153: /// <param name="D">
154: /// (input/output) DOUBLE PRECISION array, dimension (N)
155: /// On entry D contains the main diagonal of the bidiagonal
156: /// matrix.
157: /// On exit D, if INFO = 0, contains its singular values.
158: ///</param>
159: /// <param name="E">
160: /// (input) DOUBLE PRECISION array, dimension (M-1)
161: /// Contains the subdiagonal entries of the bidiagonal matrix.
162: /// On exit, E has been destroyed.
163: ///</param>
164: /// <param name="U">
165: /// (output) DOUBLE PRECISION array, dimension at least (LDQ, N)
166: /// On exit, U contains the left singular vectors.
167: ///</param>
168: /// <param name="LDU">
169: /// (input) INTEGER
170: /// On entry, leading dimension of U.
171: ///</param>
172: /// <param name="VT">
173: /// (output) DOUBLE PRECISION array, dimension at least (LDVT, M)
174: /// On exit, VT' contains the right singular vectors.
175: ///</param>
176: /// <param name="LDVT">
177: /// (input) INTEGER
178: /// On entry, leading dimension of VT.
179: ///</param>
180: /// <param name="SMLSIZ">
181: /// (input) INTEGER
182: /// On entry, maximum size of the subproblems at the
183: /// bottom of the computation tree.
184: ///</param>
185: /// <param name="IWORK">
186: /// (workspace) INTEGER work array.
187: /// Dimension must be at least (8 * N)
188: ///</param>
189: /// <param name="WORK">
190: /// (workspace) DOUBLE PRECISION work array.
191: /// Dimension must be at least (3 * M**2 + 2 * M)
192: ///</param>
193: /// <param name="INFO">
194: /// (output) INTEGER
195: /// = 0: successful exit.
196: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
197: /// .GT. 0: if INFO = 1, an singular value did not converge
198: ///</param>
199: public void Run(int N, int SQRE, ref double[] D, int offset_d, ref double[] E, int offset_e, ref double[] U, int offset_u, int LDU
200: , ref double[] VT, int offset_vt, int LDVT, int SMLSIZ, ref int[] IWORK, int offset_iwork, ref double[] WORK, int offset_work, ref int INFO)
201: {
202:
203: #region Array Index Correction
204:
205: int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_u = -1 - LDU + offset_u; int o_vt = -1 - LDVT + offset_vt;
206: int o_iwork = -1 + offset_iwork; int o_work = -1 + offset_work;
207:
208: #endregion
209:
210:
211: #region Prolog
212:
213: // *
214: // * -- LAPACK auxiliary routine (version 3.1) --
215: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
216: // * November 2006
217: // *
218: // * .. Scalar Arguments ..
219: // * ..
220: // * .. Array Arguments ..
221: // * ..
222: // *
223: // * Purpose
224: // * =======
225: // *
226: // * Using a divide and conquer approach, DLASD0 computes the singular
227: // * value decomposition (SVD) of a real upper bidiagonal N-by-M
228: // * matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
229: // * The algorithm computes orthogonal matrices U and VT such that
230: // * B = U * S * VT. The singular values S are overwritten on D.
231: // *
232: // * A related subroutine, DLASDA, computes only the singular values,
233: // * and optionally, the singular vectors in compact form.
234: // *
235: // * Arguments
236: // * =========
237: // *
238: // * N (input) INTEGER
239: // * On entry, the row dimension of the upper bidiagonal matrix.
240: // * This is also the dimension of the main diagonal array D.
241: // *
242: // * SQRE (input) INTEGER
243: // * Specifies the column dimension of the bidiagonal matrix.
244: // * = 0: The bidiagonal matrix has column dimension M = N;
245: // * = 1: The bidiagonal matrix has column dimension M = N+1;
246: // *
247: // * D (input/output) DOUBLE PRECISION array, dimension (N)
248: // * On entry D contains the main diagonal of the bidiagonal
249: // * matrix.
250: // * On exit D, if INFO = 0, contains its singular values.
251: // *
252: // * E (input) DOUBLE PRECISION array, dimension (M-1)
253: // * Contains the subdiagonal entries of the bidiagonal matrix.
254: // * On exit, E has been destroyed.
255: // *
256: // * U (output) DOUBLE PRECISION array, dimension at least (LDQ, N)
257: // * On exit, U contains the left singular vectors.
258: // *
259: // * LDU (input) INTEGER
260: // * On entry, leading dimension of U.
261: // *
262: // * VT (output) DOUBLE PRECISION array, dimension at least (LDVT, M)
263: // * On exit, VT' contains the right singular vectors.
264: // *
265: // * LDVT (input) INTEGER
266: // * On entry, leading dimension of VT.
267: // *
268: // * SMLSIZ (input) INTEGER
269: // * On entry, maximum size of the subproblems at the
270: // * bottom of the computation tree.
271: // *
272: // * IWORK (workspace) INTEGER work array.
273: // * Dimension must be at least (8 * N)
274: // *
275: // * WORK (workspace) DOUBLE PRECISION work array.
276: // * Dimension must be at least (3 * M**2 + 2 * M)
277: // *
278: // * INFO (output) INTEGER
279: // * = 0: successful exit.
280: // * < 0: if INFO = -i, the i-th argument had an illegal value.
281: // * > 0: if INFO = 1, an singular value did not converge
282: // *
283: // * Further Details
284: // * ===============
285: // *
286: // * Based on contributions by
287: // * Ming Gu and Huan Ren, Computer Science Division, University of
288: // * California at Berkeley, USA
289: // *
290: // * =====================================================================
291: // *
292: // * .. Local Scalars ..
293: // * ..
294: // * .. External Subroutines ..
295: // * ..
296: // * .. Executable Statements ..
297: // *
298: // * Test the input parameters.
299: // *
300:
301: #endregion
302:
303:
304: #region Body
305:
306: INFO = 0;
307: // *
308: if (N < 0)
309: {
310: INFO = - 1;
311: }
312: else
313: {
314: if ((SQRE < 0) || (SQRE > 1))
315: {
316: INFO = - 2;
317: }
318: }
319: // *
320: M = N + SQRE;
321: // *
322: if (LDU < N)
323: {
324: INFO = - 6;
325: }
326: else
327: {
328: if (LDVT < M)
329: {
330: INFO = - 8;
331: }
332: else
333: {
334: if (SMLSIZ < 3)
335: {
336: INFO = - 9;
337: }
338: }
339: }
340: if (INFO != 0)
341: {
342: this._xerbla.Run("DLASD0", - INFO);
343: return;
344: }
345: // *
346: // * If the input matrix is too small, call DLASDQ to find the SVD.
347: // *
348: if (N <= SMLSIZ)
349: {
350: this._dlasdq.Run("U", SQRE, N, M, N, 0
351: , ref D, offset_d, ref E, offset_e, ref VT, offset_vt, LDVT, ref U, offset_u, LDU
352: , ref U, offset_u, LDU, ref WORK, offset_work, ref INFO);
353: return;
354: }
355: // *
356: // * Set up the computation tree.
357: // *
358: INODE = 1;
359: NDIML = INODE + N;
360: NDIMR = NDIML + N;
361: IDXQ = NDIMR + N;
362: IWK = IDXQ + N;
363: this._dlasdt.Run(N, ref NLVL, ref ND, ref IWORK, INODE + o_iwork, ref IWORK, NDIML + o_iwork, ref IWORK, NDIMR + o_iwork
364: , SMLSIZ);
365: // *
366: // * For the nodes on bottom level of the tree, solve
367: // * their subproblems by DLASDQ.
368: // *
369: NDB1 = (ND + 1) / 2;
370: NCC = 0;
371: for (I = NDB1; I <= ND; I++)
372: {
373: // *
374: // * IC : center row of each node
375: // * NL : number of rows of left subproblem
376: // * NR : number of rows of right subproblem
377: // * NLF: starting row of the left subproblem
378: // * NRF: starting row of the right subproblem
379: // *
380: I1 = I - 1;
381: IC = IWORK[INODE + I1 + o_iwork];
382: NL = IWORK[NDIML + I1 + o_iwork];
383: NLP1 = NL + 1;
384: NR = IWORK[NDIMR + I1 + o_iwork];
385: NRP1 = NR + 1;
386: NLF = IC - NL;
387: NRF = IC + 1;
388: SQREI = 1;
389: this._dlasdq.Run("U", SQREI, NL, NLP1, NL, NCC
390: , ref D, NLF + o_d, ref E, NLF + o_e, ref VT, NLF+NLF * LDVT + o_vt, LDVT, ref U, NLF+NLF * LDU + o_u, LDU
391: , ref U, NLF+NLF * LDU + o_u, LDU, ref WORK, offset_work, ref INFO);
392: if (INFO != 0)
393: {
394: return;
395: }
396: ITEMP = IDXQ + NLF - 2;
397: for (J = 1; J <= NL; J++)
398: {
399: IWORK[ITEMP + J + o_iwork] = J;
400: }
401: if (I == ND)
402: {
403: SQREI = SQRE;
404: }
405: else
406: {
407: SQREI = 1;
408: }
409: NRP1 = NR + SQREI;
410: this._dlasdq.Run("U", SQREI, NR, NRP1, NR, NCC
411: , ref D, NRF + o_d, ref E, NRF + o_e, ref VT, NRF+NRF * LDVT + o_vt, LDVT, ref U, NRF+NRF * LDU + o_u, LDU
412: , ref U, NRF+NRF * LDU + o_u, LDU, ref WORK, offset_work, ref INFO);
413: if (INFO != 0)
414: {
415: return;
416: }
417: ITEMP = IDXQ + IC;
418: for (J = 1; J <= NR; J++)
419: {
420: IWORK[ITEMP + J - 1 + o_iwork] = J;
421: }
422: }
423: // *
424: // * Now conquer each subproblem bottom-up.
425: // *
426: for (LVL = NLVL; LVL >= 1; LVL += - 1)
427: {
428: // *
429: // * Find the first node LF and last node LL on the
430: // * current level LVL.
431: // *
432: if (LVL == 1)
433: {
434: LF = 1;
435: LL = 1;
436: }
437: else
438: {
439: LF = (int)Math.Pow(2, LVL - 1);
440: LL = 2 * LF - 1;
441: }
442: for (I = LF; I <= LL; I++)
443: {
444: IM1 = I - 1;
445: IC = IWORK[INODE + IM1 + o_iwork];
446: NL = IWORK[NDIML + IM1 + o_iwork];
447: NR = IWORK[NDIMR + IM1 + o_iwork];
448: NLF = IC - NL;
449: if ((SQRE == 0) && (I == LL))
450: {
451: SQREI = SQRE;
452: }
453: else
454: {
455: SQREI = 1;
456: }
457: IDXQC = IDXQ + NLF - 1;
458: ALPHA = D[IC + o_d];
459: BETA = E[IC + o_e];
460: this._dlasd1.Run(NL, NR, SQREI, ref D, NLF + o_d, ref ALPHA, ref BETA
461: , ref U, NLF+NLF * LDU + o_u, LDU, ref VT, NLF+NLF * LDVT + o_vt, LDVT, ref IWORK, IDXQC + o_iwork, ref IWORK, IWK + o_iwork
462: , ref WORK, offset_work, ref INFO);
463: if (INFO != 0)
464: {
465: return;
466: }
467: }
468: }
469: // *
470: return;
471: // *
472: // * End of DLASD0
473: // *
474:
475: #endregion
476:
477: }
478: }
479: }