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: /// DBDSDC computes the singular value decomposition (SVD) of a real
27: /// N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
28: /// using a divide and conquer method, where S is a diagonal matrix
29: /// with non-negative diagonal elements (the singular values of B), and
30: /// U and VT are orthogonal matrices of left and right singular vectors,
31: /// respectively. DBDSDC can be used to compute all singular values,
32: /// and optionally, singular vectors or singular vectors in compact form.
33: ///
34: /// This code makes very mild assumptions about floating point
35: /// arithmetic. It will work on machines with a guard digit in
36: /// add/subtract, or on those binary machines without guard digits
37: /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
38: /// It could conceivably fail on hexadecimal or decimal machines
39: /// without guard digits, but we know of none. See DLASD3 for details.
40: ///
41: /// The code currently calls DLASDQ if singular values only are desired.
42: /// However, it can be slightly modified to compute singular values
43: /// using the divide and conquer method.
44: ///
45: ///</summary>
46: public class DBDSDC
47: {
48:
49:
50: #region Dependencies
51:
52: LSAME _lsame; ILAENV _ilaenv; DLAMCH _dlamch; DLANST _dlanst; DCOPY _dcopy; DLARTG _dlartg; DLASCL _dlascl;
53: DLASD0 _dlasd0;DLASDA _dlasda; DLASDQ _dlasdq; DLASET _dlaset; DLASR _dlasr; DSWAP _dswap; XERBLA _xerbla;
54:
55: #endregion
56:
57:
58: #region Fields
59:
60: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; const double TWO = 2.0E+0; int DIFL = 0; int DIFR = 0;
61: int GIVCOL = 0;int GIVNUM = 0; int GIVPTR = 0; int I = 0; int IC = 0; int ICOMPQ = 0; int IERR = 0; int II = 0;
62: int IS = 0;int IU = 0; int IUPLO = 0; int IVT = 0; int J = 0; int K = 0; int KK = 0; int MLVL = 0; int NM1 = 0;
63: int NSIZE = 0;int PERM = 0; int POLES = 0; int QSTART = 0; int SMLSIZ = 0; int SMLSZP = 0; int SQRE = 0; int START = 0;
64: int WSTART = 0;int Z = 0; double CS = 0; double EPS = 0; double ORGNRM = 0; double P = 0; double R = 0; double SN = 0;
65:
66: #endregion
67:
68: public DBDSDC(LSAME lsame, ILAENV ilaenv, DLAMCH dlamch, DLANST dlanst, DCOPY dcopy, DLARTG dlartg, DLASCL dlascl, DLASD0 dlasd0, DLASDA dlasda, DLASDQ dlasdq
69: , DLASET dlaset, DLASR dlasr, DSWAP dswap, XERBLA xerbla)
70: {
71:
72:
73: #region Set Dependencies
74:
75: this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlanst = dlanst; this._dcopy = dcopy;
76: this._dlartg = dlartg;this._dlascl = dlascl; this._dlasd0 = dlasd0; this._dlasda = dlasda; this._dlasdq = dlasdq;
77: this._dlaset = dlaset;this._dlasr = dlasr; this._dswap = dswap; this._xerbla = xerbla;
78:
79: #endregion
80:
81: }
82:
83: public DBDSDC()
84: {
85:
86:
87: #region Dependencies (Initialization)
88:
89: LSAME lsame = new LSAME();
90: IEEECK ieeeck = new IEEECK();
91: IPARMQ iparmq = new IPARMQ();
92: DLAMC3 dlamc3 = new DLAMC3();
93: DLASSQ dlassq = new DLASSQ();
94: DCOPY dcopy = new DCOPY();
95: XERBLA xerbla = new XERBLA();
96: DLAMRG dlamrg = new DLAMRG();
97: DLAPY2 dlapy2 = new DLAPY2();
98: DROT drot = new DROT();
99: DNRM2 dnrm2 = new DNRM2();
100: DLASD5 dlasd5 = new DLASD5();
101: DLAS2 dlas2 = new DLAS2();
102: DLASQ5 dlasq5 = new DLASQ5();
103: DLAZQ4 dlazq4 = new DLAZQ4();
104: DSCAL dscal = new DSCAL();
105: DSWAP dswap = new DSWAP();
106: DLASDT dlasdt = new DLASDT();
107: DDOT ddot = new DDOT();
108: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
109: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
110: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
111: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
112: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
113: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
114: DLANST dlanst = new DLANST(lsame, dlassq);
115: DLARTG dlartg = new DLARTG(dlamch);
116: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
117: DLACPY dlacpy = new DLACPY(lsame);
118: DLASET dlaset = new DLASET(lsame);
119: DLASD2 dlasd2 = new DLASD2(dlamch, dlapy2, dcopy, dlacpy, dlamrg, dlaset, drot, xerbla);
120: DGEMM dgemm = new DGEMM(lsame, xerbla);
121: DLAED6 dlaed6 = new DLAED6(dlamch);
122: DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
123: DLASD3 dlasd3 = new DLASD3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlascl, dlasd4, xerbla);
124: DLASD1 dlasd1 = new DLASD1(dlamrg, dlascl, dlasd2, dlasd3, xerbla);
125: DLASQ6 dlasq6 = new DLASQ6(dlamch);
126: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
127: DLASRT dlasrt = new DLASRT(lsame, xerbla);
128: DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
129: DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
130: DLASR dlasr = new DLASR(lsame, xerbla);
131: DLASV2 dlasv2 = new DLASV2(dlamch);
132: DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap
133: , xerbla);
134: DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);
135: DLASD0 dlasd0 = new DLASD0(dlasd1, dlasdq, dlasdt, xerbla);
136: DLASD7 dlasd7 = new DLASD7(dcopy, dlamrg, drot, xerbla, dlamch, dlapy2);
137: DLASD8 dlasd8 = new DLASD8(dcopy, dlascl, dlasd4, dlaset, xerbla, ddot, dlamc3, dnrm2);
138: DLASD6 dlasd6 = new DLASD6(dcopy, dlamrg, dlascl, dlasd7, dlasd8, xerbla);
139: DLASDA dlasda = new DLASDA(dcopy, dlasd6, dlasdq, dlasdt, dlaset, xerbla);
140:
141: #endregion
142:
143:
144: #region Set Dependencies
145:
146: this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlanst = dlanst; this._dcopy = dcopy;
147: this._dlartg = dlartg;this._dlascl = dlascl; this._dlasd0 = dlasd0; this._dlasda = dlasda; this._dlasdq = dlasdq;
148: this._dlaset = dlaset;this._dlasr = dlasr; this._dswap = dswap; this._xerbla = xerbla;
149:
150: #endregion
151:
152: }
153: /// <summary>
154: /// Purpose
155: /// =======
156: ///
157: /// DBDSDC computes the singular value decomposition (SVD) of a real
158: /// N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
159: /// using a divide and conquer method, where S is a diagonal matrix
160: /// with non-negative diagonal elements (the singular values of B), and
161: /// U and VT are orthogonal matrices of left and right singular vectors,
162: /// respectively. DBDSDC can be used to compute all singular values,
163: /// and optionally, singular vectors or singular vectors in compact form.
164: ///
165: /// This code makes very mild assumptions about floating point
166: /// arithmetic. It will work on machines with a guard digit in
167: /// add/subtract, or on those binary machines without guard digits
168: /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
169: /// It could conceivably fail on hexadecimal or decimal machines
170: /// without guard digits, but we know of none. See DLASD3 for details.
171: ///
172: /// The code currently calls DLASDQ if singular values only are desired.
173: /// However, it can be slightly modified to compute singular values
174: /// using the divide and conquer method.
175: ///
176: ///</summary>
177: /// <param name="UPLO">
178: /// (input) CHARACTER*1
179: /// = 'U': B is upper bidiagonal.
180: /// = 'L': B is lower bidiagonal.
181: ///</param>
182: /// <param name="COMPQ">
183: /// (input) CHARACTER*1
184: /// Specifies whether singular vectors are to be computed
185: /// as follows:
186: /// = 'N': Compute singular values only;
187: /// = 'P': Compute singular values and compute singular
188: /// vectors in compact form;
189: /// = 'I': Compute singular values and singular vectors.
190: ///</param>
191: /// <param name="N">
192: /// (input) INTEGER
193: /// The order of the matrix B. N .GE. 0.
194: ///</param>
195: /// <param name="D">
196: /// (input/output) DOUBLE PRECISION array, dimension (N)
197: /// On entry, the n diagonal elements of the bidiagonal matrix B.
198: /// On exit, if INFO=0, the singular values of B.
199: ///</param>
200: /// <param name="E">
201: /// (input/output) DOUBLE PRECISION array, dimension (N-1)
202: /// On entry, the elements of E contain the offdiagonal
203: /// elements of the bidiagonal matrix whose SVD is desired.
204: /// On exit, E has been destroyed.
205: ///</param>
206: /// <param name="U">
207: /// (output) DOUBLE PRECISION array, dimension (LDU,N)
208: /// If COMPQ = 'I', then:
209: /// On exit, if INFO = 0, U contains the left singular vectors
210: /// of the bidiagonal matrix.
211: /// For other values of COMPQ, U is not referenced.
212: ///</param>
213: /// <param name="LDU">
214: /// (input) INTEGER
215: /// The leading dimension of the array U. LDU .GE. 1.
216: /// If singular vectors are desired, then LDU .GE. max( 1, N ).
217: ///</param>
218: /// <param name="VT">
219: /// (output) DOUBLE PRECISION array, dimension (LDVT,N)
220: /// If COMPQ = 'I', then:
221: /// On exit, if INFO = 0, VT' contains the right singular
222: /// vectors of the bidiagonal matrix.
223: /// For other values of COMPQ, VT is not referenced.
224: ///</param>
225: /// <param name="LDVT">
226: /// (input) INTEGER
227: /// The leading dimension of the array VT. LDVT .GE. 1.
228: /// If singular vectors are desired, then LDVT .GE. max( 1, N ).
229: ///</param>
230: /// <param name="Q">
231: /// (output) DOUBLE PRECISION array, dimension (LDQ)
232: /// If COMPQ = 'P', then:
233: /// On exit, if INFO = 0, Q and IQ contain the left
234: /// and right singular vectors in a compact form,
235: /// requiring O(N log N) space instead of 2*N**2.
236: /// In particular, Q contains all the DOUBLE PRECISION data in
237: /// LDQ .GE. N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
238: /// words of memory, where SMLSIZ is returned by ILAENV and
239: /// is equal to the maximum size of the subproblems at the
240: /// bottom of the computation tree (usually about 25).
241: /// For other values of COMPQ, Q is not referenced.
242: ///</param>
243: /// <param name="IQ">
244: /// (output) INTEGER array, dimension (LDIQ)
245: /// If COMPQ = 'P', then:
246: /// On exit, if INFO = 0, Q and IQ contain the left
247: /// and right singular vectors in a compact form,
248: /// requiring O(N log N) space instead of 2*N**2.
249: /// In particular, IQ contains all INTEGER data in
250: /// LDIQ .GE. N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
251: /// words of memory, where SMLSIZ is returned by ILAENV and
252: /// is equal to the maximum size of the subproblems at the
253: /// bottom of the computation tree (usually about 25).
254: /// For other values of COMPQ, IQ is not referenced.
255: ///</param>
256: /// <param name="WORK">
257: /// (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
258: /// If COMPQ = 'N' then LWORK .GE. (4 * N).
259: /// If COMPQ = 'P' then LWORK .GE. (6 * N).
260: /// If COMPQ = 'I' then LWORK .GE. (3 * N**2 + 4 * N).
261: ///</param>
262: /// <param name="IWORK">
263: /// (workspace) INTEGER array, dimension (8*N)
264: ///</param>
265: /// <param name="INFO">
266: /// (output) INTEGER
267: /// = 0: successful exit.
268: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
269: /// .GT. 0: The algorithm failed to compute an singular value.
270: /// The update process of divide and conquer failed.
271: ///</param>
272: public void Run(string UPLO, string COMPQ, int N, ref double[] D, int offset_d, ref double[] E, int offset_e, ref double[] U, int offset_u
273: , int LDU, ref double[] VT, int offset_vt, int LDVT, ref double[] Q, int offset_q, ref int[] IQ, int offset_iq, ref double[] WORK, int offset_work
274: , ref int[] IWORK, int offset_iwork, ref int INFO)
275: {
276:
277: #region Array Index Correction
278:
279: 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;
280: int o_q = -1 + offset_q; int o_iq = -1 + offset_iq; int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
281:
282: #endregion
283:
284:
285: #region Strings
286:
287: UPLO = UPLO.Substring(0, 1); COMPQ = COMPQ.Substring(0, 1);
288:
289: #endregion
290:
291:
292: #region Prolog
293:
294: // *
295: // * -- LAPACK routine (version 3.1) --
296: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
297: // * November 2006
298: // *
299: // * .. Scalar Arguments ..
300: // * ..
301: // * .. Array Arguments ..
302: // * ..
303: // *
304: // * Purpose
305: // * =======
306: // *
307: // * DBDSDC computes the singular value decomposition (SVD) of a real
308: // * N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
309: // * using a divide and conquer method, where S is a diagonal matrix
310: // * with non-negative diagonal elements (the singular values of B), and
311: // * U and VT are orthogonal matrices of left and right singular vectors,
312: // * respectively. DBDSDC can be used to compute all singular values,
313: // * and optionally, singular vectors or singular vectors in compact form.
314: // *
315: // * This code makes very mild assumptions about floating point
316: // * arithmetic. It will work on machines with a guard digit in
317: // * add/subtract, or on those binary machines without guard digits
318: // * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
319: // * It could conceivably fail on hexadecimal or decimal machines
320: // * without guard digits, but we know of none. See DLASD3 for details.
321: // *
322: // * The code currently calls DLASDQ if singular values only are desired.
323: // * However, it can be slightly modified to compute singular values
324: // * using the divide and conquer method.
325: // *
326: // * Arguments
327: // * =========
328: // *
329: // * UPLO (input) CHARACTER*1
330: // * = 'U': B is upper bidiagonal.
331: // * = 'L': B is lower bidiagonal.
332: // *
333: // * COMPQ (input) CHARACTER*1
334: // * Specifies whether singular vectors are to be computed
335: // * as follows:
336: // * = 'N': Compute singular values only;
337: // * = 'P': Compute singular values and compute singular
338: // * vectors in compact form;
339: // * = 'I': Compute singular values and singular vectors.
340: // *
341: // * N (input) INTEGER
342: // * The order of the matrix B. N >= 0.
343: // *
344: // * D (input/output) DOUBLE PRECISION array, dimension (N)
345: // * On entry, the n diagonal elements of the bidiagonal matrix B.
346: // * On exit, if INFO=0, the singular values of B.
347: // *
348: // * E (input/output) DOUBLE PRECISION array, dimension (N-1)
349: // * On entry, the elements of E contain the offdiagonal
350: // * elements of the bidiagonal matrix whose SVD is desired.
351: // * On exit, E has been destroyed.
352: // *
353: // * U (output) DOUBLE PRECISION array, dimension (LDU,N)
354: // * If COMPQ = 'I', then:
355: // * On exit, if INFO = 0, U contains the left singular vectors
356: // * of the bidiagonal matrix.
357: // * For other values of COMPQ, U is not referenced.
358: // *
359: // * LDU (input) INTEGER
360: // * The leading dimension of the array U. LDU >= 1.
361: // * If singular vectors are desired, then LDU >= max( 1, N ).
362: // *
363: // * VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
364: // * If COMPQ = 'I', then:
365: // * On exit, if INFO = 0, VT' contains the right singular
366: // * vectors of the bidiagonal matrix.
367: // * For other values of COMPQ, VT is not referenced.
368: // *
369: // * LDVT (input) INTEGER
370: // * The leading dimension of the array VT. LDVT >= 1.
371: // * If singular vectors are desired, then LDVT >= max( 1, N ).
372: // *
373: // * Q (output) DOUBLE PRECISION array, dimension (LDQ)
374: // * If COMPQ = 'P', then:
375: // * On exit, if INFO = 0, Q and IQ contain the left
376: // * and right singular vectors in a compact form,
377: // * requiring O(N log N) space instead of 2*N**2.
378: // * In particular, Q contains all the DOUBLE PRECISION data in
379: // * LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
380: // * words of memory, where SMLSIZ is returned by ILAENV and
381: // * is equal to the maximum size of the subproblems at the
382: // * bottom of the computation tree (usually about 25).
383: // * For other values of COMPQ, Q is not referenced.
384: // *
385: // * IQ (output) INTEGER array, dimension (LDIQ)
386: // * If COMPQ = 'P', then:
387: // * On exit, if INFO = 0, Q and IQ contain the left
388: // * and right singular vectors in a compact form,
389: // * requiring O(N log N) space instead of 2*N**2.
390: // * In particular, IQ contains all INTEGER data in
391: // * LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
392: // * words of memory, where SMLSIZ is returned by ILAENV and
393: // * is equal to the maximum size of the subproblems at the
394: // * bottom of the computation tree (usually about 25).
395: // * For other values of COMPQ, IQ is not referenced.
396: // *
397: // * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
398: // * If COMPQ = 'N' then LWORK >= (4 * N).
399: // * If COMPQ = 'P' then LWORK >= (6 * N).
400: // * If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
401: // *
402: // * IWORK (workspace) INTEGER array, dimension (8*N)
403: // *
404: // * INFO (output) INTEGER
405: // * = 0: successful exit.
406: // * < 0: if INFO = -i, the i-th argument had an illegal value.
407: // * > 0: The algorithm failed to compute an singular value.
408: // * The update process of divide and conquer failed.
409: // *
410: // * Further Details
411: // * ===============
412: // *
413: // * Based on contributions by
414: // * Ming Gu and Huan Ren, Computer Science Division, University of
415: // * California at Berkeley, USA
416: // *
417: // * =====================================================================
418: // * Changed dimension statement in comment describing E from (N) to
419: // * (N-1). Sven, 17 Feb 05.
420: // * =====================================================================
421: // *
422: // * .. Parameters ..
423: // * ..
424: // * .. Local Scalars ..
425: // * ..
426: // * .. External Functions ..
427: // * ..
428: // * .. External Subroutines ..
429: // * ..
430: // * .. Intrinsic Functions ..
431: // INTRINSIC ABS, DBLE, INT, LOG, SIGN;
432: // * ..
433: // * .. Executable Statements ..
434: // *
435: // * Test the input parameters.
436: // *
437:
438: #endregion
439:
440:
441: #region Body
442:
443: INFO = 0;
444: // *
445: IUPLO = 0;
446: if (this._lsame.Run(UPLO, "U")) IUPLO = 1;
447: if (this._lsame.Run(UPLO, "L")) IUPLO = 2;
448: if (this._lsame.Run(COMPQ, "N"))
449: {
450: ICOMPQ = 0;
451: }
452: else
453: {
454: if (this._lsame.Run(COMPQ, "P"))
455: {
456: ICOMPQ = 1;
457: }
458: else
459: {
460: if (this._lsame.Run(COMPQ, "I"))
461: {
462: ICOMPQ = 2;
463: }
464: else
465: {
466: ICOMPQ = - 1;
467: }
468: }
469: }
470: if (IUPLO == 0)
471: {
472: INFO = - 1;
473: }
474: else
475: {
476: if (ICOMPQ < 0)
477: {
478: INFO = - 2;
479: }
480: else
481: {
482: if (N < 0)
483: {
484: INFO = - 3;
485: }
486: else
487: {
488: if ((LDU < 1) || ((ICOMPQ == 2) && (LDU < N)))
489: {
490: INFO = - 7;
491: }
492: else
493: {
494: if ((LDVT < 1) || ((ICOMPQ == 2) && (LDVT < N)))
495: {
496: INFO = - 9;
497: }
498: }
499: }
500: }
501: }
502: if (INFO != 0)
503: {
504: this._xerbla.Run("DBDSDC", - INFO);
505: return;
506: }
507: // *
508: // * Quick return if possible
509: // *
510: if (N == 0) return;
511: SMLSIZ = this._ilaenv.Run(9, "DBDSDC", " ", 0, 0, 0, 0);
512: if (N == 1)
513: {
514: if (ICOMPQ == 1)
515: {
516: Q[1 + o_q] = FortranLib.Sign(ONE,D[1 + o_d]);
517: Q[1 + SMLSIZ * N + o_q] = ONE;
518: }
519: else
520: {
521: if (ICOMPQ == 2)
522: {
523: U[1+1 * LDU + o_u] = FortranLib.Sign(ONE,D[1 + o_d]);
524: VT[1+1 * LDVT + o_vt] = ONE;
525: }
526: }
527: D[1 + o_d] = Math.Abs(D[1 + o_d]);
528: return;
529: }
530: NM1 = N - 1;
531: // *
532: // * If matrix lower bidiagonal, rotate to be upper bidiagonal
533: // * by applying Givens rotations on the left
534: // *
535: WSTART = 1;
536: QSTART = 3;
537: if (ICOMPQ == 1)
538: {
539: this._dcopy.Run(N, D, offset_d, 1, ref Q, 1 + o_q, 1);
540: this._dcopy.Run(N - 1, E, offset_e, 1, ref Q, N + 1 + o_q, 1);
541: }
542: if (IUPLO == 2)
543: {
544: QSTART = 5;
545: WSTART = 2 * N - 1;
546: for (I = 1; I <= N - 1; I++)
547: {
548: this._dlartg.Run(D[I + o_d], E[I + o_e], ref CS, ref SN, ref R);
549: D[I + o_d] = R;
550: E[I + o_e] = SN * D[I + 1 + o_d];
551: D[I + 1 + o_d] = CS * D[I + 1 + o_d];
552: if (ICOMPQ == 1)
553: {
554: Q[I + 2 * N + o_q] = CS;
555: Q[I + 3 * N + o_q] = SN;
556: }
557: else
558: {
559: if (ICOMPQ == 2)
560: {
561: WORK[I + o_work] = CS;
562: WORK[NM1 + I + o_work] = - SN;
563: }
564: }
565: }
566: }
567: // *
568: // * If ICOMPQ = 0, use DLASDQ to compute the singular values.
569: // *
570: if (ICOMPQ == 0)
571: {
572: this._dlasdq.Run("U", 0, N, 0, 0, 0
573: , ref D, offset_d, ref E, offset_e, ref VT, offset_vt, LDVT, ref U, offset_u, LDU
574: , ref U, offset_u, LDU, ref WORK, WSTART + o_work, ref INFO);
575: goto LABEL40;
576: }
577: // *
578: // * If N is smaller than the minimum divide size SMLSIZ, then solve
579: // * the problem with another solver.
580: // *
581: if (N <= SMLSIZ)
582: {
583: if (ICOMPQ == 2)
584: {
585: this._dlaset.Run("A", N, N, ZERO, ONE, ref U, offset_u
586: , LDU);
587: this._dlaset.Run("A", N, N, ZERO, ONE, ref VT, offset_vt
588: , LDVT);
589: this._dlasdq.Run("U", 0, N, N, N, 0
590: , ref D, offset_d, ref E, offset_e, ref VT, offset_vt, LDVT, ref U, offset_u, LDU
591: , ref U, offset_u, LDU, ref WORK, WSTART + o_work, ref INFO);
592: }
593: else
594: {
595: if (ICOMPQ == 1)
596: {
597: IU = 1;
598: IVT = IU + N;
599: this._dlaset.Run("A", N, N, ZERO, ONE, ref Q, IU + (QSTART - 1) * N + o_q
600: , N);
601: this._dlaset.Run("A", N, N, ZERO, ONE, ref Q, IVT + (QSTART - 1) * N + o_q
602: , N);
603: this._dlasdq.Run("U", 0, N, N, N, 0
604: , ref D, offset_d, ref E, offset_e, ref Q, IVT + (QSTART - 1) * N + o_q, N, ref Q, IU + (QSTART - 1) * N + o_q, N
605: , ref Q, IU + (QSTART - 1) * N + o_q, N, ref WORK, WSTART + o_work, ref INFO);
606: }
607: }
608: goto LABEL40;
609: }
610: // *
611: if (ICOMPQ == 2)
612: {
613: this._dlaset.Run("A", N, N, ZERO, ONE, ref U, offset_u
614: , LDU);
615: this._dlaset.Run("A", N, N, ZERO, ONE, ref VT, offset_vt
616: , LDVT);
617: }
618: // *
619: // * Scale.
620: // *
621: ORGNRM = this._dlanst.Run("M", N, D, offset_d, E, offset_e);
622: if (ORGNRM == ZERO) return;
623: this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N
624: , 1, ref D, offset_d, N, ref IERR);
625: this._dlascl.Run("G", 0, 0, ORGNRM, ONE, NM1
626: , 1, ref E, offset_e, NM1, ref IERR);
627: // *
628: EPS = this._dlamch.Run("Epsilon");
629: // *
630: MLVL = Convert.ToInt32(Math.Truncate(Math.Log(Convert.ToDouble(N) / Convert.ToDouble(SMLSIZ + 1)) / Math.Log(TWO))) + 1;
631: SMLSZP = SMLSIZ + 1;
632: // *
633: if (ICOMPQ == 1)
634: {
635: IU = 1;
636: IVT = 1 + SMLSIZ;
637: DIFL = IVT + SMLSZP;
638: DIFR = DIFL + MLVL;
639: Z = DIFR + MLVL * 2;
640: IC = Z + MLVL;
641: IS = IC + 1;
642: POLES = IS + 1;
643: GIVNUM = POLES + 2 * MLVL;
644: // *
645: K = 1;
646: GIVPTR = 2;
647: PERM = 3;
648: GIVCOL = PERM + MLVL;
649: }
650: // *
651: for (I = 1; I <= N; I++)
652: {
653: if (Math.Abs(D[I + o_d]) < EPS)
654: {
655: D[I + o_d] = FortranLib.Sign(EPS,D[I + o_d]);
656: }
657: }
658: // *
659: START = 1;
660: SQRE = 0;
661: // *
662: for (I = 1; I <= NM1; I++)
663: {
664: if ((Math.Abs(E[I + o_e]) < EPS) || (I == NM1))
665: {
666: // *
667: // * Subproblem found. First determine its size and then
668: // * apply divide and conquer on it.
669: // *
670: if (I < NM1)
671: {
672: // *
673: // * A subproblem with E(I) small for I < NM1.
674: // *
675: NSIZE = I - START + 1;
676: }
677: else
678: {
679: if (Math.Abs(E[I + o_e]) >= EPS)
680: {
681: // *
682: // * A subproblem with E(NM1) not too small but I = NM1.
683: // *
684: NSIZE = N - START + 1;
685: }
686: else
687: {
688: // *
689: // * A subproblem with E(NM1) small. This implies an
690: // * 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
691: // * first.
692: // *
693: NSIZE = I - START + 1;
694: if (ICOMPQ == 2)
695: {
696: U[N+N * LDU + o_u] = FortranLib.Sign(ONE,D[N + o_d]);
697: VT[N+N * LDVT + o_vt] = ONE;
698: }
699: else
700: {
701: if (ICOMPQ == 1)
702: {
703: Q[N + (QSTART - 1) * N + o_q] = FortranLib.Sign(ONE,D[N + o_d]);
704: Q[N + (SMLSIZ + QSTART - 1) * N + o_q] = ONE;
705: }
706: }
707: D[N + o_d] = Math.Abs(D[N + o_d]);
708: }
709: }
710: if (ICOMPQ == 2)
711: {
712: this._dlasd0.Run(NSIZE, SQRE, ref D, START + o_d, ref E, START + o_e, ref U, START+START * LDU + o_u, LDU
713: , ref VT, START+START * LDVT + o_vt, LDVT, SMLSIZ, ref IWORK, offset_iwork, ref WORK, WSTART + o_work, ref INFO);
714: }
715: else
716: {
717: this._dlasda.Run(ICOMPQ, SMLSIZ, NSIZE, SQRE, ref D, START + o_d, ref E, START + o_e
718: , ref Q, START + (IU + QSTART - 2) * N + o_q, N, ref Q, START + (IVT + QSTART - 2) * N + o_q, ref IQ, START + K * N + o_iq, ref Q, START + (DIFL + QSTART - 2) * N + o_q, ref Q, START + (DIFR + QSTART - 2) * N + o_q
719: , ref Q, START + (Z + QSTART - 2) * N + o_q, ref Q, START + (POLES + QSTART - 2) * N + o_q, ref IQ, START + GIVPTR * N + o_iq, ref IQ, START + GIVCOL * N + o_iq, N, ref IQ, START + PERM * N + o_iq
720: , ref Q, START + (GIVNUM + QSTART - 2) * N + o_q, ref Q, START + (IC + QSTART - 2) * N + o_q, ref Q, START + (IS + QSTART - 2) * N + o_q, ref WORK, WSTART + o_work, ref IWORK, offset_iwork, ref INFO);
721: if (INFO != 0)
722: {
723: return;
724: }
725: }
726: START = I + 1;
727: }
728: }
729: // *
730: // * Unscale
731: // *
732: this._dlascl.Run("G", 0, 0, ONE, ORGNRM, N
733: , 1, ref D, offset_d, N, ref IERR);
734: LABEL40:;
735: // *
736: // * Use Selection Sort to minimize swaps of singular vectors
737: // *
738: for (II = 2; II <= N; II++)
739: {
740: I = II - 1;
741: KK = I;
742: P = D[I + o_d];
743: for (J = II; J <= N; J++)
744: {
745: if (D[J + o_d] > P)
746: {
747: KK = J;
748: P = D[J + o_d];
749: }
750: }
751: if (KK != I)
752: {
753: D[KK + o_d] = D[I + o_d];
754: D[I + o_d] = P;
755: if (ICOMPQ == 1)
756: {
757: IQ[I + o_iq] = KK;
758: }
759: else
760: {
761: if (ICOMPQ == 2)
762: {
763: this._dswap.Run(N, ref U, 1+I * LDU + o_u, 1, ref U, 1+KK * LDU + o_u, 1);
764: this._dswap.Run(N, ref VT, I+1 * LDVT + o_vt, LDVT, ref VT, KK+1 * LDVT + o_vt, LDVT);
765: }
766: }
767: }
768: else
769: {
770: if (ICOMPQ == 1)
771: {
772: IQ[I + o_iq] = I;
773: }
774: }
775: }
776: // *
777: // * If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
778: // *
779: if (ICOMPQ == 1)
780: {
781: if (IUPLO == 1)
782: {
783: IQ[N + o_iq] = 1;
784: }
785: else
786: {
787: IQ[N + o_iq] = 0;
788: }
789: }
790: // *
791: // * If B is lower bidiagonal, update U by those Givens rotations
792: // * which rotated B to be upper bidiagonal
793: // *
794: if ((IUPLO == 2) && (ICOMPQ == 2))
795: {
796: this._dlasr.Run("L", "V", "B", N, N, WORK, 1 + o_work
797: , WORK, N + o_work, ref U, offset_u, LDU);
798: }
799: // *
800: return;
801: // *
802: // * End of DBDSDC
803: // *
804:
805: #endregion
806:
807: }
808: }
809: }