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, DLASDA computes the singular
27: /// value decomposition (SVD) of a real upper bidiagonal N-by-M matrix
28: /// B with diagonal D and offdiagonal E, where M = N + SQRE. The
29: /// algorithm computes the singular values in the SVD B = U * S * VT.
30: /// The orthogonal matrices U and VT are optionally computed in
31: /// compact form.
32: ///
33: /// A related subroutine, DLASD0, computes the singular values and
34: /// the singular vectors in explicit form.
35: ///
36: ///</summary>
37: public class DLASDA
38: {
39:
40:
41: #region Dependencies
42:
43: DCOPY _dcopy; DLASD6 _dlasd6; DLASDQ _dlasdq; DLASDT _dlasdt; DLASET _dlaset; XERBLA _xerbla;
44:
45: #endregion
46:
47:
48: #region Fields
49:
50: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int I = 0; int I1 = 0; int IC = 0; int IDXQ = 0; int IDXQI = 0;
51: int IM1 = 0;int INODE = 0; int ITEMP = 0; int IWK = 0; int J = 0; int LF = 0; int LL = 0; int LVL = 0; int LVL2 = 0;
52: int M = 0;int NCC = 0; int ND = 0; int NDB1 = 0; int NDIML = 0; int NDIMR = 0; int NL = 0; int NLF = 0; int NLP1 = 0;
53: int NLVL = 0;int NR = 0; int NRF = 0; int NRP1 = 0; int NRU = 0; int NWORK1 = 0; int NWORK2 = 0; int SMLSZP = 0;
54: int SQREI = 0;int VF = 0; int VFI = 0; int VL = 0; int VLI = 0; double ALPHA = 0; double BETA = 0;
55:
56: #endregion
57:
58: public DLASDA(DCOPY dcopy, DLASD6 dlasd6, DLASDQ dlasdq, DLASDT dlasdt, DLASET dlaset, XERBLA xerbla)
59: {
60:
61:
62: #region Set Dependencies
63:
64: this._dcopy = dcopy; this._dlasd6 = dlasd6; this._dlasdq = dlasdq; this._dlasdt = dlasdt; this._dlaset = dlaset;
65: this._xerbla = xerbla;
66:
67: #endregion
68:
69: }
70:
71: public DLASDA()
72: {
73:
74:
75: #region Dependencies (Initialization)
76:
77: DCOPY dcopy = new DCOPY();
78: DLAMRG dlamrg = new DLAMRG();
79: LSAME lsame = new LSAME();
80: DLAMC3 dlamc3 = new DLAMC3();
81: XERBLA xerbla = new XERBLA();
82: DROT drot = new DROT();
83: DLAPY2 dlapy2 = new DLAPY2();
84: DLASD5 dlasd5 = new DLASD5();
85: DDOT ddot = new DDOT();
86: DNRM2 dnrm2 = new DNRM2();
87: DLAS2 dlas2 = new DLAS2();
88: DLASQ5 dlasq5 = new DLASQ5();
89: DLAZQ4 dlazq4 = new DLAZQ4();
90: IEEECK ieeeck = new IEEECK();
91: IPARMQ iparmq = new IPARMQ();
92: DSCAL dscal = new DSCAL();
93: DSWAP dswap = new DSWAP();
94: DLASDT dlasdt = new DLASDT();
95: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
96: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
97: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
98: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
99: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
100: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
101: DLASD7 dlasd7 = new DLASD7(dcopy, dlamrg, drot, xerbla, dlamch, dlapy2);
102: DLAED6 dlaed6 = new DLAED6(dlamch);
103: DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
104: DLASET dlaset = new DLASET(lsame);
105: DLASD8 dlasd8 = new DLASD8(dcopy, dlascl, dlasd4, dlaset, xerbla, ddot, dlamc3, dnrm2);
106: DLASD6 dlasd6 = new DLASD6(dcopy, dlamrg, dlascl, dlasd7, dlasd8, xerbla);
107: DLARTG dlartg = new DLARTG(dlamch);
108: DLASQ6 dlasq6 = new DLASQ6(dlamch);
109: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
110: DLASRT dlasrt = new DLASRT(lsame, xerbla);
111: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
112: DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
113: DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
114: DLASR dlasr = new DLASR(lsame, xerbla);
115: DLASV2 dlasv2 = new DLASV2(dlamch);
116: DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap
117: , xerbla);
118: DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);
119:
120: #endregion
121:
122:
123: #region Set Dependencies
124:
125: this._dcopy = dcopy; this._dlasd6 = dlasd6; this._dlasdq = dlasdq; this._dlasdt = dlasdt; this._dlaset = dlaset;
126: this._xerbla = xerbla;
127:
128: #endregion
129:
130: }
131: /// <summary>
132: /// Purpose
133: /// =======
134: ///
135: /// Using a divide and conquer approach, DLASDA computes the singular
136: /// value decomposition (SVD) of a real upper bidiagonal N-by-M matrix
137: /// B with diagonal D and offdiagonal E, where M = N + SQRE. The
138: /// algorithm computes the singular values in the SVD B = U * S * VT.
139: /// The orthogonal matrices U and VT are optionally computed in
140: /// compact form.
141: ///
142: /// A related subroutine, DLASD0, computes the singular values and
143: /// the singular vectors in explicit form.
144: ///
145: ///</summary>
146: /// <param name="ICOMPQ">
147: /// (input) INTEGER
148: /// Specifies whether singular vectors are to be computed
149: /// in compact form, as follows
150: /// = 0: Compute singular values only.
151: /// = 1: Compute singular vectors of upper bidiagonal
152: /// matrix in compact form.
153: ///</param>
154: /// <param name="SMLSIZ">
155: /// (input) INTEGER
156: /// The maximum size of the subproblems at the bottom of the
157: /// computation tree.
158: ///</param>
159: /// <param name="N">
160: /// (input) INTEGER
161: /// The row dimension of the upper bidiagonal matrix. This is
162: /// also the dimension of the main diagonal array D.
163: ///</param>
164: /// <param name="SQRE">
165: /// (input) INTEGER
166: /// Specifies the column dimension of the bidiagonal matrix.
167: /// = 0: The bidiagonal matrix has column dimension M = N;
168: /// = 1: The bidiagonal matrix has column dimension M = N + 1.
169: ///</param>
170: /// <param name="D">
171: /// (input/output) DOUBLE PRECISION array, dimension ( N )
172: /// On entry D contains the main diagonal of the bidiagonal
173: /// matrix. On exit D, if INFO = 0, contains its singular values.
174: ///</param>
175: /// <param name="E">
176: /// (input) DOUBLE PRECISION array, dimension ( M-1 )
177: /// Contains the subdiagonal entries of the bidiagonal matrix.
178: /// On exit, E has been destroyed.
179: ///</param>
180: /// <param name="U">
181: /// (output) DOUBLE PRECISION array,
182: /// dimension ( LDU, SMLSIZ ) if ICOMPQ = 1, and not referenced
183: /// if ICOMPQ = 0. If ICOMPQ = 1, on exit, U contains the left
184: /// singular vector matrices of all subproblems at the bottom
185: /// level.
186: ///</param>
187: /// <param name="LDU">
188: /// (input) INTEGER, LDU = .GT. N.
189: /// The leading dimension of arrays U, VT, DIFL, DIFR, POLES,
190: /// GIVNUM, and Z.
191: ///</param>
192: /// <param name="VT">
193: /// (output) DOUBLE PRECISION array,
194: /// dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced
195: /// if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT' contains the right
196: /// singular vector matrices of all subproblems at the bottom
197: /// level.
198: ///</param>
199: /// <param name="K">
200: /// (output) INTEGER array,
201: /// dimension ( N ) if ICOMPQ = 1 and dimension 1 if ICOMPQ = 0.
202: /// If ICOMPQ = 1, on exit, K(I) is the dimension of the I-th
203: /// secular equation on the computation tree.
204: ///</param>
205: /// <param name="DIFL">
206: /// (output) DOUBLE PRECISION array, dimension ( LDU, NLVL ),
207: /// where NLVL = floor(log_2 (N/SMLSIZ))).
208: ///</param>
209: /// <param name="DIFR">
210: /// (output) DOUBLE PRECISION array,
211: /// dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1 and
212: /// dimension ( N ) if ICOMPQ = 0.
213: /// If ICOMPQ = 1, on exit, DIFL(1:N, I) and DIFR(1:N, 2 * I - 1)
214: /// record distances between singular values on the I-th
215: /// level and singular values on the (I -1)-th level, and
216: /// DIFR(1:N, 2 * I ) contains the normalizing factors for
217: /// the right singular vector matrix. See DLASD8 for details.
218: ///</param>
219: /// <param name="Z">
220: /// (output) DOUBLE PRECISION array,
221: /// dimension ( LDU, NLVL ) if ICOMPQ = 1 and
222: /// dimension ( N ) if ICOMPQ = 0.
223: /// The first K elements of Z(1, I) contain the components of
224: /// the deflation-adjusted updating row vector for subproblems
225: /// on the I-th level.
226: ///</param>
227: /// <param name="POLES">
228: /// (output) DOUBLE PRECISION array,
229: /// dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced
230: /// if ICOMPQ = 0. If ICOMPQ = 1, on exit, POLES(1, 2*I - 1) and
231: /// POLES(1, 2*I) contain the new and old singular values
232: /// involved in the secular equations on the I-th level.
233: ///</param>
234: /// <param name="GIVPTR">
235: /// (output) INTEGER array,
236: /// dimension ( N ) if ICOMPQ = 1, and not referenced if
237: /// ICOMPQ = 0. If ICOMPQ = 1, on exit, GIVPTR( I ) records
238: /// the number of Givens rotations performed on the I-th
239: /// problem on the computation tree.
240: ///</param>
241: /// <param name="GIVCOL">
242: /// (output) INTEGER array,
243: /// dimension ( LDGCOL, 2 * NLVL ) if ICOMPQ = 1, and not
244: /// referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
245: /// GIVCOL(1, 2 *I - 1) and GIVCOL(1, 2 *I) record the locations
246: /// of Givens rotations performed on the I-th level on the
247: /// computation tree.
248: ///</param>
249: /// <param name="LDGCOL">
250: /// (input) INTEGER, LDGCOL = .GT. N.
251: /// The leading dimension of arrays GIVCOL and PERM.
252: ///</param>
253: /// <param name="PERM">
254: /// (output) INTEGER array,
255: /// dimension ( LDGCOL, NLVL ) if ICOMPQ = 1, and not referenced
256: /// if ICOMPQ = 0. If ICOMPQ = 1, on exit, PERM(1, I) records
257: /// permutations done on the I-th level of the computation tree.
258: ///</param>
259: /// <param name="GIVNUM">
260: /// (output) DOUBLE PRECISION array,
261: /// dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not
262: /// referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
263: /// GIVNUM(1, 2 *I - 1) and GIVNUM(1, 2 *I) record the C- and S-
264: /// values of Givens rotations performed on the I-th level on
265: /// the computation tree.
266: ///</param>
267: /// <param name="C">
268: /// (output) DOUBLE PRECISION array,
269: /// dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0.
270: /// If ICOMPQ = 1 and the I-th subproblem is not square, on exit,
271: /// C( I ) contains the C-value of a Givens rotation related to
272: /// the right null space of the I-th subproblem.
273: ///</param>
274: /// <param name="S">
275: /// (output) DOUBLE PRECISION array, dimension ( N ) if
276: /// ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1
277: /// and the I-th subproblem is not square, on exit, S( I )
278: /// contains the S-value of a Givens rotation related to
279: /// the right null space of the I-th subproblem.
280: ///</param>
281: /// <param name="WORK">
282: /// (workspace) DOUBLE PRECISION array, dimension
283: /// (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)).
284: ///</param>
285: /// <param name="IWORK">
286: /// (workspace) INTEGER array.
287: /// Dimension must be at least (7 * N).
288: ///</param>
289: /// <param name="INFO">
290: /// (output) INTEGER
291: /// = 0: successful exit.
292: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
293: /// .GT. 0: if INFO = 1, an singular value did not converge
294: ///</param>
295: public void Run(int ICOMPQ, int SMLSIZ, int N, int SQRE, ref double[] D, int offset_d, ref double[] E, int offset_e
296: , ref double[] U, int offset_u, int LDU, ref double[] VT, int offset_vt, ref int[] K, int offset_k, ref double[] DIFL, int offset_difl, ref double[] DIFR, int offset_difr
297: , ref double[] Z, int offset_z, ref double[] POLES, int offset_poles, ref int[] GIVPTR, int offset_givptr, ref int[] GIVCOL, int offset_givcol, int LDGCOL, ref int[] PERM, int offset_perm
298: , ref double[] GIVNUM, int offset_givnum, ref double[] C, int offset_c, ref double[] S, int offset_s, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork, ref int INFO)
299: {
300:
301: #region Array Index Correction
302:
303: int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_u = -1 - LDU + offset_u; int o_vt = -1 - LDU + offset_vt;
304: int o_k = -1 + offset_k; int o_difl = -1 - LDU + offset_difl; int o_difr = -1 - LDU + offset_difr;
305: int o_z = -1 - LDU + offset_z; int o_poles = -1 - LDU + offset_poles; int o_givptr = -1 + offset_givptr;
306: int o_givcol = -1 - LDGCOL + offset_givcol; int o_perm = -1 - LDGCOL + offset_perm;
307: int o_givnum = -1 - LDU + offset_givnum; int o_c = -1 + offset_c; int o_s = -1 + offset_s;
308: int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
309:
310: #endregion
311:
312:
313: #region Prolog
314:
315: // *
316: // * -- LAPACK auxiliary routine (version 3.1) --
317: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
318: // * November 2006
319: // *
320: // * .. Scalar Arguments ..
321: // * ..
322: // * .. Array Arguments ..
323: // * ..
324: // *
325: // * Purpose
326: // * =======
327: // *
328: // * Using a divide and conquer approach, DLASDA computes the singular
329: // * value decomposition (SVD) of a real upper bidiagonal N-by-M matrix
330: // * B with diagonal D and offdiagonal E, where M = N + SQRE. The
331: // * algorithm computes the singular values in the SVD B = U * S * VT.
332: // * The orthogonal matrices U and VT are optionally computed in
333: // * compact form.
334: // *
335: // * A related subroutine, DLASD0, computes the singular values and
336: // * the singular vectors in explicit form.
337: // *
338: // * Arguments
339: // * =========
340: // *
341: // * ICOMPQ (input) INTEGER
342: // * Specifies whether singular vectors are to be computed
343: // * in compact form, as follows
344: // * = 0: Compute singular values only.
345: // * = 1: Compute singular vectors of upper bidiagonal
346: // * matrix in compact form.
347: // *
348: // * SMLSIZ (input) INTEGER
349: // * The maximum size of the subproblems at the bottom of the
350: // * computation tree.
351: // *
352: // * N (input) INTEGER
353: // * The row dimension of the upper bidiagonal matrix. This is
354: // * also the dimension of the main diagonal array D.
355: // *
356: // * SQRE (input) INTEGER
357: // * Specifies the column dimension of the bidiagonal matrix.
358: // * = 0: The bidiagonal matrix has column dimension M = N;
359: // * = 1: The bidiagonal matrix has column dimension M = N + 1.
360: // *
361: // * D (input/output) DOUBLE PRECISION array, dimension ( N )
362: // * On entry D contains the main diagonal of the bidiagonal
363: // * matrix. On exit D, if INFO = 0, contains its singular values.
364: // *
365: // * E (input) DOUBLE PRECISION array, dimension ( M-1 )
366: // * Contains the subdiagonal entries of the bidiagonal matrix.
367: // * On exit, E has been destroyed.
368: // *
369: // * U (output) DOUBLE PRECISION array,
370: // * dimension ( LDU, SMLSIZ ) if ICOMPQ = 1, and not referenced
371: // * if ICOMPQ = 0. If ICOMPQ = 1, on exit, U contains the left
372: // * singular vector matrices of all subproblems at the bottom
373: // * level.
374: // *
375: // * LDU (input) INTEGER, LDU = > N.
376: // * The leading dimension of arrays U, VT, DIFL, DIFR, POLES,
377: // * GIVNUM, and Z.
378: // *
379: // * VT (output) DOUBLE PRECISION array,
380: // * dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced
381: // * if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT' contains the right
382: // * singular vector matrices of all subproblems at the bottom
383: // * level.
384: // *
385: // * K (output) INTEGER array,
386: // * dimension ( N ) if ICOMPQ = 1 and dimension 1 if ICOMPQ = 0.
387: // * If ICOMPQ = 1, on exit, K(I) is the dimension of the I-th
388: // * secular equation on the computation tree.
389: // *
390: // * DIFL (output) DOUBLE PRECISION array, dimension ( LDU, NLVL ),
391: // * where NLVL = floor(log_2 (N/SMLSIZ))).
392: // *
393: // * DIFR (output) DOUBLE PRECISION array,
394: // * dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1 and
395: // * dimension ( N ) if ICOMPQ = 0.
396: // * If ICOMPQ = 1, on exit, DIFL(1:N, I) and DIFR(1:N, 2 * I - 1)
397: // * record distances between singular values on the I-th
398: // * level and singular values on the (I -1)-th level, and
399: // * DIFR(1:N, 2 * I ) contains the normalizing factors for
400: // * the right singular vector matrix. See DLASD8 for details.
401: // *
402: // * Z (output) DOUBLE PRECISION array,
403: // * dimension ( LDU, NLVL ) if ICOMPQ = 1 and
404: // * dimension ( N ) if ICOMPQ = 0.
405: // * The first K elements of Z(1, I) contain the components of
406: // * the deflation-adjusted updating row vector for subproblems
407: // * on the I-th level.
408: // *
409: // * POLES (output) DOUBLE PRECISION array,
410: // * dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced
411: // * if ICOMPQ = 0. If ICOMPQ = 1, on exit, POLES(1, 2*I - 1) and
412: // * POLES(1, 2*I) contain the new and old singular values
413: // * involved in the secular equations on the I-th level.
414: // *
415: // * GIVPTR (output) INTEGER array,
416: // * dimension ( N ) if ICOMPQ = 1, and not referenced if
417: // * ICOMPQ = 0. If ICOMPQ = 1, on exit, GIVPTR( I ) records
418: // * the number of Givens rotations performed on the I-th
419: // * problem on the computation tree.
420: // *
421: // * GIVCOL (output) INTEGER array,
422: // * dimension ( LDGCOL, 2 * NLVL ) if ICOMPQ = 1, and not
423: // * referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
424: // * GIVCOL(1, 2 *I - 1) and GIVCOL(1, 2 *I) record the locations
425: // * of Givens rotations performed on the I-th level on the
426: // * computation tree.
427: // *
428: // * LDGCOL (input) INTEGER, LDGCOL = > N.
429: // * The leading dimension of arrays GIVCOL and PERM.
430: // *
431: // * PERM (output) INTEGER array,
432: // * dimension ( LDGCOL, NLVL ) if ICOMPQ = 1, and not referenced
433: // * if ICOMPQ = 0. If ICOMPQ = 1, on exit, PERM(1, I) records
434: // * permutations done on the I-th level of the computation tree.
435: // *
436: // * GIVNUM (output) DOUBLE PRECISION array,
437: // * dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not
438: // * referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
439: // * GIVNUM(1, 2 *I - 1) and GIVNUM(1, 2 *I) record the C- and S-
440: // * values of Givens rotations performed on the I-th level on
441: // * the computation tree.
442: // *
443: // * C (output) DOUBLE PRECISION array,
444: // * dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0.
445: // * If ICOMPQ = 1 and the I-th subproblem is not square, on exit,
446: // * C( I ) contains the C-value of a Givens rotation related to
447: // * the right null space of the I-th subproblem.
448: // *
449: // * S (output) DOUBLE PRECISION array, dimension ( N ) if
450: // * ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1
451: // * and the I-th subproblem is not square, on exit, S( I )
452: // * contains the S-value of a Givens rotation related to
453: // * the right null space of the I-th subproblem.
454: // *
455: // * WORK (workspace) DOUBLE PRECISION array, dimension
456: // * (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)).
457: // *
458: // * IWORK (workspace) INTEGER array.
459: // * Dimension must be at least (7 * N).
460: // *
461: // * INFO (output) INTEGER
462: // * = 0: successful exit.
463: // * < 0: if INFO = -i, the i-th argument had an illegal value.
464: // * > 0: if INFO = 1, an singular value did not converge
465: // *
466: // * Further Details
467: // * ===============
468: // *
469: // * Based on contributions by
470: // * Ming Gu and Huan Ren, Computer Science Division, University of
471: // * California at Berkeley, USA
472: // *
473: // * =====================================================================
474: // *
475: // * .. Parameters ..
476: // * ..
477: // * .. Local Scalars ..
478: // * ..
479: // * .. External Subroutines ..
480: // * ..
481: // * .. Executable Statements ..
482: // *
483: // * Test the input parameters.
484: // *
485:
486: #endregion
487:
488:
489: #region Body
490:
491: INFO = 0;
492: // *
493: if ((ICOMPQ < 0) || (ICOMPQ > 1))
494: {
495: INFO = - 1;
496: }
497: else
498: {
499: if (SMLSIZ < 3)
500: {
501: INFO = - 2;
502: }
503: else
504: {
505: if (N < 0)
506: {
507: INFO = - 3;
508: }
509: else
510: {
511: if ((SQRE < 0) || (SQRE > 1))
512: {
513: INFO = - 4;
514: }
515: else
516: {
517: if (LDU < (N + SQRE))
518: {
519: INFO = - 8;
520: }
521: else
522: {
523: if (LDGCOL < N)
524: {
525: INFO = - 17;
526: }
527: }
528: }
529: }
530: }
531: }
532: if (INFO != 0)
533: {
534: this._xerbla.Run("DLASDA", - INFO);
535: return;
536: }
537: // *
538: M = N + SQRE;
539: // *
540: // * If the input matrix is too small, call DLASDQ to find the SVD.
541: // *
542: if (N <= SMLSIZ)
543: {
544: if (ICOMPQ == 0)
545: {
546: this._dlasdq.Run("U", SQRE, N, 0, 0, 0
547: , ref D, offset_d, ref E, offset_e, ref VT, offset_vt, LDU, ref U, offset_u, LDU
548: , ref U, offset_u, LDU, ref WORK, offset_work, ref INFO);
549: }
550: else
551: {
552: this._dlasdq.Run("U", SQRE, N, M, N, 0
553: , ref D, offset_d, ref E, offset_e, ref VT, offset_vt, LDU, ref U, offset_u, LDU
554: , ref U, offset_u, LDU, ref WORK, offset_work, ref INFO);
555: }
556: return;
557: }
558: // *
559: // * Book-keeping and set up the computation tree.
560: // *
561: INODE = 1;
562: NDIML = INODE + N;
563: NDIMR = NDIML + N;
564: IDXQ = NDIMR + N;
565: IWK = IDXQ + N;
566: // *
567: NCC = 0;
568: NRU = 0;
569: // *
570: SMLSZP = SMLSIZ + 1;
571: VF = 1;
572: VL = VF + M;
573: NWORK1 = VL + M;
574: NWORK2 = NWORK1 + SMLSZP * SMLSZP;
575: // *
576: this._dlasdt.Run(N, ref NLVL, ref ND, ref IWORK, INODE + o_iwork, ref IWORK, NDIML + o_iwork, ref IWORK, NDIMR + o_iwork
577: , SMLSIZ);
578: // *
579: // * for the nodes on bottom level of the tree, solve
580: // * their subproblems by DLASDQ.
581: // *
582: NDB1 = (ND + 1) / 2;
583: for (I = NDB1; I <= ND; I++)
584: {
585: // *
586: // * IC : center row of each node
587: // * NL : number of rows of left subproblem
588: // * NR : number of rows of right subproblem
589: // * NLF: starting row of the left subproblem
590: // * NRF: starting row of the right subproblem
591: // *
592: I1 = I - 1;
593: IC = IWORK[INODE + I1 + o_iwork];
594: NL = IWORK[NDIML + I1 + o_iwork];
595: NLP1 = NL + 1;
596: NR = IWORK[NDIMR + I1 + o_iwork];
597: NLF = IC - NL;
598: NRF = IC + 1;
599: IDXQI = IDXQ + NLF - 2;
600: VFI = VF + NLF - 1;
601: VLI = VL + NLF - 1;
602: SQREI = 1;
603: if (ICOMPQ == 0)
604: {
605: this._dlaset.Run("A", NLP1, NLP1, ZERO, ONE, ref WORK, NWORK1 + o_work
606: , SMLSZP);
607: this._dlasdq.Run("U", SQREI, NL, NLP1, NRU, NCC
608: , ref D, NLF + o_d, ref E, NLF + o_e, ref WORK, NWORK1 + o_work, SMLSZP, ref WORK, NWORK2 + o_work, NL
609: , ref WORK, NWORK2 + o_work, NL, ref WORK, NWORK2 + o_work, ref INFO);
610: ITEMP = NWORK1 + NL * SMLSZP;
611: this._dcopy.Run(NLP1, WORK, NWORK1 + o_work, 1, ref WORK, VFI + o_work, 1);
612: this._dcopy.Run(NLP1, WORK, ITEMP + o_work, 1, ref WORK, VLI + o_work, 1);
613: }
614: else
615: {
616: this._dlaset.Run("A", NL, NL, ZERO, ONE, ref U, NLF+1 * LDU + o_u
617: , LDU);
618: this._dlaset.Run("A", NLP1, NLP1, ZERO, ONE, ref VT, NLF+1 * LDU + o_vt
619: , LDU);
620: this._dlasdq.Run("U", SQREI, NL, NLP1, NL, NCC
621: , ref D, NLF + o_d, ref E, NLF + o_e, ref VT, NLF+1 * LDU + o_vt, LDU, ref U, NLF+1 * LDU + o_u, LDU
622: , ref U, NLF+1 * LDU + o_u, LDU, ref WORK, NWORK1 + o_work, ref INFO);
623: this._dcopy.Run(NLP1, VT, NLF+1 * LDU + o_vt, 1, ref WORK, VFI + o_work, 1);
624: this._dcopy.Run(NLP1, VT, NLF+NLP1 * LDU + o_vt, 1, ref WORK, VLI + o_work, 1);
625: }
626: if (INFO != 0)
627: {
628: return;
629: }
630: for (J = 1; J <= NL; J++)
631: {
632: IWORK[IDXQI + J + o_iwork] = J;
633: }
634: if ((I == ND) && (SQRE == 0))
635: {
636: SQREI = 0;
637: }
638: else
639: {
640: SQREI = 1;
641: }
642: IDXQI = IDXQI + NLP1;
643: VFI = VFI + NLP1;
644: VLI = VLI + NLP1;
645: NRP1 = NR + SQREI;
646: if (ICOMPQ == 0)
647: {
648: this._dlaset.Run("A", NRP1, NRP1, ZERO, ONE, ref WORK, NWORK1 + o_work
649: , SMLSZP);
650: this._dlasdq.Run("U", SQREI, NR, NRP1, NRU, NCC
651: , ref D, NRF + o_d, ref E, NRF + o_e, ref WORK, NWORK1 + o_work, SMLSZP, ref WORK, NWORK2 + o_work, NR
652: , ref WORK, NWORK2 + o_work, NR, ref WORK, NWORK2 + o_work, ref INFO);
653: ITEMP = NWORK1 + (NRP1 - 1) * SMLSZP;
654: this._dcopy.Run(NRP1, WORK, NWORK1 + o_work, 1, ref WORK, VFI + o_work, 1);
655: this._dcopy.Run(NRP1, WORK, ITEMP + o_work, 1, ref WORK, VLI + o_work, 1);
656: }
657: else
658: {
659: this._dlaset.Run("A", NR, NR, ZERO, ONE, ref U, NRF+1 * LDU + o_u
660: , LDU);
661: this._dlaset.Run("A", NRP1, NRP1, ZERO, ONE, ref VT, NRF+1 * LDU + o_vt
662: , LDU);
663: this._dlasdq.Run("U", SQREI, NR, NRP1, NR, NCC
664: , ref D, NRF + o_d, ref E, NRF + o_e, ref VT, NRF+1 * LDU + o_vt, LDU, ref U, NRF+1 * LDU + o_u, LDU
665: , ref U, NRF+1 * LDU + o_u, LDU, ref WORK, NWORK1 + o_work, ref INFO);
666: this._dcopy.Run(NRP1, VT, NRF+1 * LDU + o_vt, 1, ref WORK, VFI + o_work, 1);
667: this._dcopy.Run(NRP1, VT, NRF+NRP1 * LDU + o_vt, 1, ref WORK, VLI + o_work, 1);
668: }
669: if (INFO != 0)
670: {
671: return;
672: }
673: for (J = 1; J <= NR; J++)
674: {
675: IWORK[IDXQI + J + o_iwork] = J;
676: }
677: }
678: // *
679: // * Now conquer each subproblem bottom-up.
680: // *
681: J = (int)Math.Pow(2, NLVL);
682: for (LVL = NLVL; LVL >= 1; LVL += - 1)
683: {
684: LVL2 = LVL * 2 - 1;
685: // *
686: // * Find the first node LF and last node LL on
687: // * the current level LVL.
688: // *
689: if (LVL == 1)
690: {
691: LF = 1;
692: LL = 1;
693: }
694: else
695: {
696: LF = (int)Math.Pow(2, LVL - 1);
697: LL = 2 * LF - 1;
698: }
699: for (I = LF; I <= LL; I++)
700: {
701: IM1 = I - 1;
702: IC = IWORK[INODE + IM1 + o_iwork];
703: NL = IWORK[NDIML + IM1 + o_iwork];
704: NR = IWORK[NDIMR + IM1 + o_iwork];
705: NLF = IC - NL;
706: NRF = IC + 1;
707: if (I == LL)
708: {
709: SQREI = SQRE;
710: }
711: else
712: {
713: SQREI = 1;
714: }
715: VFI = VF + NLF - 1;
716: VLI = VL + NLF - 1;
717: IDXQI = IDXQ + NLF - 1;
718: ALPHA = D[IC + o_d];
719: BETA = E[IC + o_e];
720: if (ICOMPQ == 0)
721: {
722: this._dlasd6.Run(ICOMPQ, NL, NR, SQREI, ref D, NLF + o_d, ref WORK, VFI + o_work
723: , ref WORK, VLI + o_work, ref ALPHA, ref BETA, ref IWORK, IDXQI + o_iwork, ref PERM, offset_perm, ref GIVPTR[1 + o_givptr]
724: , ref GIVCOL, offset_givcol, LDGCOL, ref GIVNUM, offset_givnum, LDU, ref POLES, offset_poles, ref DIFL, offset_difl
725: , ref DIFR, offset_difr, ref Z, offset_z, ref K[1 + o_k], ref C[1 + o_c], ref S[1 + o_s], ref WORK, NWORK1 + o_work
726: , ref IWORK, IWK + o_iwork, ref INFO);
727: }
728: else
729: {
730: J = J - 1;
731: this._dlasd6.Run(ICOMPQ, NL, NR, SQREI, ref D, NLF + o_d, ref WORK, VFI + o_work
732: , ref WORK, VLI + o_work, ref ALPHA, ref BETA, ref IWORK, IDXQI + o_iwork, ref PERM, NLF+LVL * LDGCOL + o_perm, ref GIVPTR[J + o_givptr]
733: , ref GIVCOL, NLF+LVL2 * LDGCOL + o_givcol, LDGCOL, ref GIVNUM, NLF+LVL2 * LDU + o_givnum, LDU, ref POLES, NLF+LVL2 * LDU + o_poles, ref DIFL, NLF+LVL * LDU + o_difl
734: , ref DIFR, NLF+LVL2 * LDU + o_difr, ref Z, NLF+LVL * LDU + o_z, ref K[J + o_k], ref C[J + o_c], ref S[J + o_s], ref WORK, NWORK1 + o_work
735: , ref IWORK, IWK + o_iwork, ref INFO);
736: }
737: if (INFO != 0)
738: {
739: return;
740: }
741: }
742: }
743: // *
744: return;
745: // *
746: // * End of DLASDA
747: // *
748:
749: #endregion
750:
751: }
752: }
753: }