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.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// January 2007
23: /// Purpose
24: /// =======
25: ///
26: /// DBDSQR computes the singular values and, optionally, the right and/or
27: /// left singular vectors from the singular value decomposition (SVD) of
28: /// a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
29: /// zero-shift QR algorithm. The SVD of B has the form
30: ///
31: /// B = Q * S * P**T
32: ///
33: /// where S is the diagonal matrix of singular values, Q is an orthogonal
34: /// matrix of left singular vectors, and P is an orthogonal matrix of
35: /// right singular vectors. If left singular vectors are requested, this
36: /// subroutine actually returns U*Q instead of Q, and, if right singular
37: /// vectors are requested, this subroutine returns P**T*VT instead of
38: /// P**T, for given real input matrices U and VT. When U and VT are the
39: /// orthogonal matrices that reduce a general matrix A to bidiagonal
40: /// form: A = U*B*VT, as computed by DGEBRD, then
41: ///
42: /// A = (U*Q) * S * (P**T*VT)
43: ///
44: /// is the SVD of A. Optionally, the subroutine may also compute Q**T*C
45: /// for a given real input matrix C.
46: ///
47: /// See "Computing Small Singular Values of Bidiagonal Matrices With
48: /// Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
49: /// LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
50: /// no. 5, pp. 873-912, Sept 1990) and
51: /// "Accurate singular values and differential qd algorithms," by
52: /// B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
53: /// Department, University of California at Berkeley, July 1992
54: /// for a detailed description of the algorithm.
55: ///
56: ///</summary>
57: public class DBDSQR
58: {
59:
60:
61: #region Dependencies
62:
63: LSAME _lsame; DLAMCH _dlamch; DLARTG _dlartg; DLAS2 _dlas2; DLASQ1 _dlasq1; DLASR _dlasr; DLASV2 _dlasv2; DROT _drot;
64: DSCAL _dscal;DSWAP _dswap; XERBLA _xerbla;
65:
66: #endregion
67:
68:
69: #region Fields
70:
71: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double NEGONE = - 1.0E0; const double HNDRTH = 0.01E0;
72: const double TEN = 10.0E0;const double HNDRD = 100.0E0; const double MEIGTH = - 0.125E0; const int MAXITR = 6;
73: bool LOWER = false;bool ROTATE = false; int I = 0; int IDIR = 0; int ISUB = 0; int ITER = 0; int J = 0; int LL = 0;
74: int LLL = 0;int M = 0; int MAXIT = 0; int NM1 = 0; int NM12 = 0; int NM13 = 0; int OLDLL = 0; int OLDM = 0;
75: double ABSE = 0;double ABSS = 0; double COSL = 0; double COSR = 0; double CS = 0; double EPS = 0; double F = 0;
76: double G = 0;double H = 0; double MU = 0; double OLDCS = 0; double OLDSN = 0; double R = 0; double SHIFT = 0;
77: double SIGMN = 0;double SIGMX = 0; double SINL = 0; double SINR = 0; double SLL = 0; double SMAX = 0; double SMIN = 0;
78: double SMINL = 0;double SMINOA = 0; double SN = 0; double THRESH = 0; double TOL = 0; double TOLMUL = 0; double UNFL = 0;
79:
80: #endregion
81:
82: public DBDSQR(LSAME lsame, DLAMCH dlamch, DLARTG dlartg, DLAS2 dlas2, DLASQ1 dlasq1, DLASR dlasr, DLASV2 dlasv2, DROT drot, DSCAL dscal, DSWAP dswap
83: , XERBLA xerbla)
84: {
85:
86:
87: #region Set Dependencies
88:
89: this._lsame = lsame; this._dlamch = dlamch; this._dlartg = dlartg; this._dlas2 = dlas2; this._dlasq1 = dlasq1;
90: this._dlasr = dlasr;this._dlasv2 = dlasv2; this._drot = drot; this._dscal = dscal; this._dswap = dswap;
91: this._xerbla = xerbla;
92:
93: #endregion
94:
95: }
96:
97: public DBDSQR()
98: {
99:
100:
101: #region Dependencies (Initialization)
102:
103: LSAME lsame = new LSAME();
104: DLAMC3 dlamc3 = new DLAMC3();
105: DLAS2 dlas2 = new DLAS2();
106: DCOPY dcopy = new DCOPY();
107: XERBLA xerbla = new XERBLA();
108: DLASQ5 dlasq5 = new DLASQ5();
109: DLAZQ4 dlazq4 = new DLAZQ4();
110: IEEECK ieeeck = new IEEECK();
111: IPARMQ iparmq = new IPARMQ();
112: DROT drot = new DROT();
113: DSCAL dscal = new DSCAL();
114: DSWAP dswap = new DSWAP();
115: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
116: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
117: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
118: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
119: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
120: DLARTG dlartg = new DLARTG(dlamch);
121: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
122: DLASQ6 dlasq6 = new DLASQ6(dlamch);
123: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
124: DLASRT dlasrt = new DLASRT(lsame, xerbla);
125: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
126: DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
127: DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
128: DLASR dlasr = new DLASR(lsame, xerbla);
129: DLASV2 dlasv2 = new DLASV2(dlamch);
130:
131: #endregion
132:
133:
134: #region Set Dependencies
135:
136: this._lsame = lsame; this._dlamch = dlamch; this._dlartg = dlartg; this._dlas2 = dlas2; this._dlasq1 = dlasq1;
137: this._dlasr = dlasr;this._dlasv2 = dlasv2; this._drot = drot; this._dscal = dscal; this._dswap = dswap;
138: this._xerbla = xerbla;
139:
140: #endregion
141:
142: }
143: /// <summary>
144: /// Purpose
145: /// =======
146: ///
147: /// DBDSQR computes the singular values and, optionally, the right and/or
148: /// left singular vectors from the singular value decomposition (SVD) of
149: /// a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
150: /// zero-shift QR algorithm. The SVD of B has the form
151: ///
152: /// B = Q * S * P**T
153: ///
154: /// where S is the diagonal matrix of singular values, Q is an orthogonal
155: /// matrix of left singular vectors, and P is an orthogonal matrix of
156: /// right singular vectors. If left singular vectors are requested, this
157: /// subroutine actually returns U*Q instead of Q, and, if right singular
158: /// vectors are requested, this subroutine returns P**T*VT instead of
159: /// P**T, for given real input matrices U and VT. When U and VT are the
160: /// orthogonal matrices that reduce a general matrix A to bidiagonal
161: /// form: A = U*B*VT, as computed by DGEBRD, then
162: ///
163: /// A = (U*Q) * S * (P**T*VT)
164: ///
165: /// is the SVD of A. Optionally, the subroutine may also compute Q**T*C
166: /// for a given real input matrix C.
167: ///
168: /// See "Computing Small Singular Values of Bidiagonal Matrices With
169: /// Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
170: /// LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
171: /// no. 5, pp. 873-912, Sept 1990) and
172: /// "Accurate singular values and differential qd algorithms," by
173: /// B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
174: /// Department, University of California at Berkeley, July 1992
175: /// for a detailed description of the algorithm.
176: ///
177: ///</summary>
178: /// <param name="UPLO">
179: /// (input) CHARACTER*1
180: /// = 'U': B is upper bidiagonal;
181: /// = 'L': B is lower bidiagonal.
182: ///</param>
183: /// <param name="N">
184: /// (input) INTEGER
185: /// The order of the matrix B. N .GE. 0.
186: ///</param>
187: /// <param name="NCVT">
188: /// (input) INTEGER
189: /// The number of columns of the matrix VT. NCVT .GE. 0.
190: ///</param>
191: /// <param name="NRU">
192: /// (input) INTEGER
193: /// The number of rows of the matrix U. NRU .GE. 0.
194: ///</param>
195: /// <param name="NCC">
196: /// (input) INTEGER
197: /// The number of columns of the matrix C. NCC .GE. 0.
198: ///</param>
199: /// <param name="D">
200: /// (input/output) DOUBLE PRECISION array, dimension (N)
201: /// On entry, the n diagonal elements of the bidiagonal matrix B.
202: /// On exit, if INFO=0, the singular values of B in decreasing
203: /// order.
204: ///</param>
205: /// <param name="E">
206: /// (input/output) DOUBLE PRECISION array, dimension (N-1)
207: /// On entry, the N-1 offdiagonal elements of the bidiagonal
208: /// matrix B.
209: /// On exit, if INFO = 0, E is destroyed; if INFO .GT. 0, D and E
210: /// will contain the diagonal and superdiagonal elements of a
211: /// bidiagonal matrix orthogonally equivalent to the one given
212: /// as input.
213: ///</param>
214: /// <param name="VT">
215: /// (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
216: /// On entry, an N-by-NCVT matrix VT.
217: /// On exit, VT is overwritten by P**T * VT.
218: /// Not referenced if NCVT = 0.
219: ///</param>
220: /// <param name="LDVT">
221: /// (input) INTEGER
222: /// The leading dimension of the array VT.
223: /// LDVT .GE. max(1,N) if NCVT .GT. 0; LDVT .GE. 1 if NCVT = 0.
224: ///</param>
225: /// <param name="U">
226: /// (input/output) DOUBLE PRECISION array, dimension (LDU, N)
227: /// On entry, an NRU-by-N matrix U.
228: /// On exit, U is overwritten by U * Q.
229: /// Not referenced if NRU = 0.
230: ///</param>
231: /// <param name="LDU">
232: /// (input) INTEGER
233: /// The leading dimension of the array U. LDU .GE. max(1,NRU).
234: ///</param>
235: /// <param name="C">
236: /// (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
237: /// On entry, an N-by-NCC matrix C.
238: /// On exit, C is overwritten by Q**T * C.
239: /// Not referenced if NCC = 0.
240: ///</param>
241: /// <param name="LDC">
242: /// (input) INTEGER
243: /// The leading dimension of the array C.
244: /// LDC .GE. max(1,N) if NCC .GT. 0; LDC .GE.1 if NCC = 0.
245: ///</param>
246: /// <param name="WORK">
247: /// (workspace) DOUBLE PRECISION array, dimension (2*N)
248: /// if NCVT = NRU = NCC = 0, (max(1, 4*N)) otherwise
249: ///</param>
250: /// <param name="INFO">
251: /// (output) INTEGER
252: /// = 0: successful exit
253: /// .LT. 0: If INFO = -i, the i-th argument had an illegal value
254: /// .GT. 0: the algorithm did not converge; D and E contain the
255: /// elements of a bidiagonal matrix which is orthogonally
256: /// similar to the input matrix B; if INFO = i, i
257: /// elements of E have not converged to zero.
258: ///</param>
259: public void Run(string UPLO, int N, int NCVT, int NRU, int NCC, ref double[] D, int offset_d
260: , ref double[] E, int offset_e, ref double[] VT, int offset_vt, int LDVT, ref double[] U, int offset_u, int LDU, ref double[] C, int offset_c
261: , int LDC, ref double[] WORK, int offset_work, ref int INFO)
262: {
263:
264: #region Array Index Correction
265:
266: int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_vt = -1 - LDVT + offset_vt; int o_u = -1 - LDU + offset_u;
267: int o_c = -1 - LDC + offset_c; int o_work = -1 + offset_work;
268:
269: #endregion
270:
271:
272: #region Strings
273:
274: UPLO = UPLO.Substring(0, 1);
275:
276: #endregion
277:
278:
279: #region Prolog
280:
281: // *
282: // * -- LAPACK routine (version 3.1.1) --
283: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
284: // * January 2007
285: // *
286: // * .. Scalar Arguments ..
287: // * ..
288: // * .. Array Arguments ..
289: // * ..
290: // *
291: // * Purpose
292: // * =======
293: // *
294: // * DBDSQR computes the singular values and, optionally, the right and/or
295: // * left singular vectors from the singular value decomposition (SVD) of
296: // * a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
297: // * zero-shift QR algorithm. The SVD of B has the form
298: // *
299: // * B = Q * S * P**T
300: // *
301: // * where S is the diagonal matrix of singular values, Q is an orthogonal
302: // * matrix of left singular vectors, and P is an orthogonal matrix of
303: // * right singular vectors. If left singular vectors are requested, this
304: // * subroutine actually returns U*Q instead of Q, and, if right singular
305: // * vectors are requested, this subroutine returns P**T*VT instead of
306: // * P**T, for given real input matrices U and VT. When U and VT are the
307: // * orthogonal matrices that reduce a general matrix A to bidiagonal
308: // * form: A = U*B*VT, as computed by DGEBRD, then
309: // *
310: // * A = (U*Q) * S * (P**T*VT)
311: // *
312: // * is the SVD of A. Optionally, the subroutine may also compute Q**T*C
313: // * for a given real input matrix C.
314: // *
315: // * See "Computing Small Singular Values of Bidiagonal Matrices With
316: // * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
317: // * LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
318: // * no. 5, pp. 873-912, Sept 1990) and
319: // * "Accurate singular values and differential qd algorithms," by
320: // * B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
321: // * Department, University of California at Berkeley, July 1992
322: // * for a detailed description of the algorithm.
323: // *
324: // * Arguments
325: // * =========
326: // *
327: // * UPLO (input) CHARACTER*1
328: // * = 'U': B is upper bidiagonal;
329: // * = 'L': B is lower bidiagonal.
330: // *
331: // * N (input) INTEGER
332: // * The order of the matrix B. N >= 0.
333: // *
334: // * NCVT (input) INTEGER
335: // * The number of columns of the matrix VT. NCVT >= 0.
336: // *
337: // * NRU (input) INTEGER
338: // * The number of rows of the matrix U. NRU >= 0.
339: // *
340: // * NCC (input) INTEGER
341: // * The number of columns of the matrix C. NCC >= 0.
342: // *
343: // * D (input/output) DOUBLE PRECISION array, dimension (N)
344: // * On entry, the n diagonal elements of the bidiagonal matrix B.
345: // * On exit, if INFO=0, the singular values of B in decreasing
346: // * order.
347: // *
348: // * E (input/output) DOUBLE PRECISION array, dimension (N-1)
349: // * On entry, the N-1 offdiagonal elements of the bidiagonal
350: // * matrix B.
351: // * On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
352: // * will contain the diagonal and superdiagonal elements of a
353: // * bidiagonal matrix orthogonally equivalent to the one given
354: // * as input.
355: // *
356: // * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
357: // * On entry, an N-by-NCVT matrix VT.
358: // * On exit, VT is overwritten by P**T * VT.
359: // * Not referenced if NCVT = 0.
360: // *
361: // * LDVT (input) INTEGER
362: // * The leading dimension of the array VT.
363: // * LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
364: // *
365: // * U (input/output) DOUBLE PRECISION array, dimension (LDU, N)
366: // * On entry, an NRU-by-N matrix U.
367: // * On exit, U is overwritten by U * Q.
368: // * Not referenced if NRU = 0.
369: // *
370: // * LDU (input) INTEGER
371: // * The leading dimension of the array U. LDU >= max(1,NRU).
372: // *
373: // * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
374: // * On entry, an N-by-NCC matrix C.
375: // * On exit, C is overwritten by Q**T * C.
376: // * Not referenced if NCC = 0.
377: // *
378: // * LDC (input) INTEGER
379: // * The leading dimension of the array C.
380: // * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
381: // *
382: // * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
383: // * if NCVT = NRU = NCC = 0, (max(1, 4*N)) otherwise
384: // *
385: // * INFO (output) INTEGER
386: // * = 0: successful exit
387: // * < 0: If INFO = -i, the i-th argument had an illegal value
388: // * > 0: the algorithm did not converge; D and E contain the
389: // * elements of a bidiagonal matrix which is orthogonally
390: // * similar to the input matrix B; if INFO = i, i
391: // * elements of E have not converged to zero.
392: // *
393: // * Internal Parameters
394: // * ===================
395: // *
396: // * TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
397: // * TOLMUL controls the convergence criterion of the QR loop.
398: // * If it is positive, TOLMUL*EPS is the desired relative
399: // * precision in the computed singular values.
400: // * If it is negative, abs(TOLMUL*EPS*sigma_max) is the
401: // * desired absolute accuracy in the computed singular
402: // * values (corresponds to relative accuracy
403: // * abs(TOLMUL*EPS) in the largest singular value.
404: // * abs(TOLMUL) should be between 1 and 1/EPS, and preferably
405: // * between 10 (for fast convergence) and .1/EPS
406: // * (for there to be some accuracy in the results).
407: // * Default is to lose at either one eighth or 2 of the
408: // * available decimal digits in each computed singular value
409: // * (whichever is smaller).
410: // *
411: // * MAXITR INTEGER, default = 6
412: // * MAXITR controls the maximum number of passes of the
413: // * algorithm through its inner loop. The algorithms stops
414: // * (and so fails to converge) if the number of passes
415: // * through the inner loop exceeds MAXITR*N**2.
416: // *
417: // * =====================================================================
418: // *
419: // * .. Parameters ..
420: // * ..
421: // * .. Local Scalars ..
422: // * ..
423: // * .. External Functions ..
424: // * ..
425: // * .. External Subroutines ..
426: // * ..
427: // * .. Intrinsic Functions ..
428: // INTRINSIC ABS, DBLE, MAX, MIN, SIGN, SQRT;
429: // * ..
430: // * .. Executable Statements ..
431: // *
432: // * Test the input parameters.
433: // *
434:
435: #endregion
436:
437:
438: #region Body
439:
440: INFO = 0;
441: LOWER = this._lsame.Run(UPLO, "L");
442: if (!this._lsame.Run(UPLO, "U") && !LOWER)
443: {
444: INFO = - 1;
445: }
446: else
447: {
448: if (N < 0)
449: {
450: INFO = - 2;
451: }
452: else
453: {
454: if (NCVT < 0)
455: {
456: INFO = - 3;
457: }
458: else
459: {
460: if (NRU < 0)
461: {
462: INFO = - 4;
463: }
464: else
465: {
466: if (NCC < 0)
467: {
468: INFO = - 5;
469: }
470: else
471: {
472: if ((NCVT == 0 && LDVT < 1) || (NCVT > 0 && LDVT < Math.Max(1, N)))
473: {
474: INFO = - 9;
475: }
476: else
477: {
478: if (LDU < Math.Max(1, NRU))
479: {
480: INFO = - 11;
481: }
482: else
483: {
484: if ((NCC == 0 && LDC < 1) || (NCC > 0 && LDC < Math.Max(1, N)))
485: {
486: INFO = - 13;
487: }
488: }
489: }
490: }
491: }
492: }
493: }
494: }
495: if (INFO != 0)
496: {
497: this._xerbla.Run("DBDSQR", - INFO);
498: return;
499: }
500: if (N == 0) return;
501: if (N == 1) goto LABEL160;
502: // *
503: // * ROTATE is true if any singular vectors desired, false otherwise
504: // *
505: ROTATE = (NCVT > 0) || (NRU > 0) || (NCC > 0);
506: // *
507: // * If no singular vectors desired, use qd algorithm
508: // *
509: if (!ROTATE)
510: {
511: this._dlasq1.Run(N, ref D, offset_d, E, offset_e, ref WORK, offset_work, ref INFO);
512: return;
513: }
514: // *
515: NM1 = N - 1;
516: NM12 = NM1 + NM1;
517: NM13 = NM12 + NM1;
518: IDIR = 0;
519: // *
520: // * Get machine constants
521: // *
522: EPS = this._dlamch.Run("Epsilon");
523: UNFL = this._dlamch.Run("Safe minimum");
524: // *
525: // * If matrix lower bidiagonal, rotate to be upper bidiagonal
526: // * by applying Givens rotations on the left
527: // *
528: if (LOWER)
529: {
530: for (I = 1; I <= N - 1; I++)
531: {
532: this._dlartg.Run(D[I + o_d], E[I + o_e], ref CS, ref SN, ref R);
533: D[I + o_d] = R;
534: E[I + o_e] = SN * D[I + 1 + o_d];
535: D[I + 1 + o_d] = CS * D[I + 1 + o_d];
536: WORK[I + o_work] = CS;
537: WORK[NM1 + I + o_work] = SN;
538: }
539: // *
540: // * Update singular vectors if desired
541: // *
542: if (NRU > 0)
543: {
544: this._dlasr.Run("R", "V", "F", NRU, N, WORK, 1 + o_work
545: , WORK, N + o_work, ref U, offset_u, LDU);
546: }
547: if (NCC > 0)
548: {
549: this._dlasr.Run("L", "V", "F", N, NCC, WORK, 1 + o_work
550: , WORK, N + o_work, ref C, offset_c, LDC);
551: }
552: }
553: // *
554: // * Compute singular values to relative accuracy TOL
555: // * (By setting TOL to be negative, algorithm will compute
556: // * singular values to absolute accuracy ABS(TOL)*norm(input matrix))
557: // *
558: TOLMUL = Math.Max(TEN, Math.Min(HNDRD, Math.Pow(EPS,MEIGTH)));
559: TOL = TOLMUL * EPS;
560: // *
561: // * Compute approximate maximum, minimum singular values
562: // *
563: SMAX = ZERO;
564: for (I = 1; I <= N; I++)
565: {
566: SMAX = Math.Max(SMAX, Math.Abs(D[I + o_d]));
567: }
568: for (I = 1; I <= N - 1; I++)
569: {
570: SMAX = Math.Max(SMAX, Math.Abs(E[I + o_e]));
571: }
572: SMINL = ZERO;
573: if (TOL >= ZERO)
574: {
575: // *
576: // * Relative accuracy desired
577: // *
578: SMINOA = Math.Abs(D[1 + o_d]);
579: if (SMINOA == ZERO) goto LABEL50;
580: MU = SMINOA;
581: for (I = 2; I <= N; I++)
582: {
583: MU = Math.Abs(D[I + o_d]) * (MU / (MU + Math.Abs(E[I - 1 + o_e])));
584: SMINOA = Math.Min(SMINOA, MU);
585: if (SMINOA == ZERO) goto LABEL50;
586: }
587: LABEL50:;
588: SMINOA = SMINOA / Math.Sqrt(Convert.ToDouble(N));
589: THRESH = Math.Max(TOL * SMINOA, MAXITR * N * N * UNFL);
590: }
591: else
592: {
593: // *
594: // * Absolute accuracy desired
595: // *
596: THRESH = Math.Max(Math.Abs(TOL) * SMAX, MAXITR * N * N * UNFL);
597: }
598: // *
599: // * Prepare for main iteration loop for the singular values
600: // * (MAXIT is the maximum number of passes through the inner
601: // * loop permitted before nonconvergence signalled.)
602: // *
603: MAXIT = MAXITR * N * N;
604: ITER = 0;
605: OLDLL = - 1;
606: OLDM = - 1;
607: // *
608: // * M points to last element of unconverged part of matrix
609: // *
610: M = N;
611: // *
612: // * Begin main iteration loop
613: // *
614: LABEL60:;
615: // *
616: // * Check for convergence or exceeding iteration count
617: // *
618: if (M <= 1) goto LABEL160;
619: if (ITER > MAXIT) goto LABEL200;
620: // *
621: // * Find diagonal block of matrix to work on
622: // *
623: if (TOL < ZERO && Math.Abs(D[M + o_d]) <= THRESH) D[M + o_d] = ZERO;
624: SMAX = Math.Abs(D[M + o_d]);
625: SMIN = SMAX;
626: for (LLL = 1; LLL <= M - 1; LLL++)
627: {
628: LL = M - LLL;
629: ABSS = Math.Abs(D[LL + o_d]);
630: ABSE = Math.Abs(E[LL + o_e]);
631: if (TOL < ZERO && ABSS <= THRESH) D[LL + o_d] = ZERO;
632: if (ABSE <= THRESH) goto LABEL80;
633: SMIN = Math.Min(SMIN, ABSS);
634: SMAX = Math.Max(SMAX, Math.Max(ABSS, ABSE));
635: }
636: LL = 0;
637: goto LABEL90;
638: LABEL80:;
639: E[LL + o_e] = ZERO;
640: // *
641: // * Matrix splits since E(LL) = 0
642: // *
643: if (LL == M - 1)
644: {
645: // *
646: // * Convergence of bottom singular value, return to top of loop
647: // *
648: M = M - 1;
649: goto LABEL60;
650: }
651: LABEL90:;
652: LL = LL + 1;
653: // *
654: // * E(LL) through E(M-1) are nonzero, E(LL-1) is zero
655: // *
656: if (LL == M - 1)
657: {
658: // *
659: // * 2 by 2 block, handle separately
660: // *
661: this._dlasv2.Run(D[M - 1 + o_d], E[M - 1 + o_e], D[M + o_d], ref SIGMN, ref SIGMX, ref SINR
662: , ref COSR, ref SINL, ref COSL);
663: D[M - 1 + o_d] = SIGMX;
664: E[M - 1 + o_e] = ZERO;
665: D[M + o_d] = SIGMN;
666: // *
667: // * Compute singular vectors, if desired
668: // *
669: if (NCVT > 0)
670: {
671: this._drot.Run(NCVT, ref VT, M - 1+1 * LDVT + o_vt, LDVT, ref VT, M+1 * LDVT + o_vt, LDVT, COSR
672: , SINR);
673: }
674: if (NRU > 0)
675: {
676: this._drot.Run(NRU, ref U, 1+(M - 1) * LDU + o_u, 1, ref U, 1+M * LDU + o_u, 1, COSL
677: , SINL);
678: }
679: if (NCC > 0)
680: {
681: this._drot.Run(NCC, ref C, M - 1+1 * LDC + o_c, LDC, ref C, M+1 * LDC + o_c, LDC, COSL
682: , SINL);
683: }
684: M = M - 2;
685: goto LABEL60;
686: }
687: // *
688: // * If working on new submatrix, choose shift direction
689: // * (from larger end diagonal element towards smaller)
690: // *
691: if (LL > OLDM || M < OLDLL)
692: {
693: if (Math.Abs(D[LL + o_d]) >= Math.Abs(D[M + o_d]))
694: {
695: // *
696: // * Chase bulge from top (big end) to bottom (small end)
697: // *
698: IDIR = 1;
699: }
700: else
701: {
702: // *
703: // * Chase bulge from bottom (big end) to top (small end)
704: // *
705: IDIR = 2;
706: }
707: }
708: // *
709: // * Apply convergence tests
710: // *
711: if (IDIR == 1)
712: {
713: // *
714: // * Run convergence test in forward direction
715: // * First apply standard test to bottom of matrix
716: // *
717: if (Math.Abs(E[M - 1 + o_e]) <= Math.Abs(TOL) * Math.Abs(D[M + o_d]) || (TOL < ZERO && Math.Abs(E[M - 1 + o_e]) <= THRESH))
718: {
719: E[M - 1 + o_e] = ZERO;
720: goto LABEL60;
721: }
722: // *
723: if (TOL >= ZERO)
724: {
725: // *
726: // * If relative accuracy desired,
727: // * apply convergence criterion forward
728: // *
729: MU = Math.Abs(D[LL + o_d]);
730: SMINL = MU;
731: for (LLL = LL; LLL <= M - 1; LLL++)
732: {
733: if (Math.Abs(E[LLL + o_e]) <= TOL * MU)
734: {
735: E[LLL + o_e] = ZERO;
736: goto LABEL60;
737: }
738: MU = Math.Abs(D[LLL + 1 + o_d]) * (MU / (MU + Math.Abs(E[LLL + o_e])));
739: SMINL = Math.Min(SMINL, MU);
740: }
741: }
742: // *
743: }
744: else
745: {
746: // *
747: // * Run convergence test in backward direction
748: // * First apply standard test to top of matrix
749: // *
750: if (Math.Abs(E[LL + o_e]) <= Math.Abs(TOL) * Math.Abs(D[LL + o_d]) || (TOL < ZERO && Math.Abs(E[LL + o_e]) <= THRESH))
751: {
752: E[LL + o_e] = ZERO;
753: goto LABEL60;
754: }
755: // *
756: if (TOL >= ZERO)
757: {
758: // *
759: // * If relative accuracy desired,
760: // * apply convergence criterion backward
761: // *
762: MU = Math.Abs(D[M + o_d]);
763: SMINL = MU;
764: for (LLL = M - 1; LLL >= LL; LLL += - 1)
765: {
766: if (Math.Abs(E[LLL + o_e]) <= TOL * MU)
767: {
768: E[LLL + o_e] = ZERO;
769: goto LABEL60;
770: }
771: MU = Math.Abs(D[LLL + o_d]) * (MU / (MU + Math.Abs(E[LLL + o_e])));
772: SMINL = Math.Min(SMINL, MU);
773: }
774: }
775: }
776: OLDLL = LL;
777: OLDM = M;
778: // *
779: // * Compute shift. First, test if shifting would ruin relative
780: // * accuracy, and if so set the shift to zero.
781: // *
782: if (TOL >= ZERO && N * TOL * (SMINL / SMAX) <= Math.Max(EPS, HNDRTH * TOL))
783: {
784: // *
785: // * Use a zero shift to avoid loss of relative accuracy
786: // *
787: SHIFT = ZERO;
788: }
789: else
790: {
791: // *
792: // * Compute the shift from 2-by-2 block at end of matrix
793: // *
794: if (IDIR == 1)
795: {
796: SLL = Math.Abs(D[LL + o_d]);
797: this._dlas2.Run(D[M - 1 + o_d], E[M - 1 + o_e], D[M + o_d], ref SHIFT, ref R);
798: }
799: else
800: {
801: SLL = Math.Abs(D[M + o_d]);
802: this._dlas2.Run(D[LL + o_d], E[LL + o_e], D[LL + 1 + o_d], ref SHIFT, ref R);
803: }
804: // *
805: // * Test if shift negligible, and if so set to zero
806: // *
807: if (SLL > ZERO)
808: {
809: if (Math.Pow(SHIFT / SLL,2) < EPS) SHIFT = ZERO;
810: }
811: }
812: // *
813: // * Increment iteration count
814: // *
815: ITER = ITER + M - LL;
816: // *
817: // * If SHIFT = 0, do simplified QR iteration
818: // *
819: if (SHIFT == ZERO)
820: {
821: if (IDIR == 1)
822: {
823: // *
824: // * Chase bulge from top to bottom
825: // * Save cosines and sines for later singular vector updates
826: // *
827: CS = ONE;
828: OLDCS = ONE;
829: for (I = LL; I <= M - 1; I++)
830: {
831: this._dlartg.Run(D[I + o_d] * CS, E[I + o_e], ref CS, ref SN, ref R);
832: if (I > LL) E[I - 1 + o_e] = OLDSN * R;
833: this._dlartg.Run(OLDCS * R, D[I + 1 + o_d] * SN, ref OLDCS, ref OLDSN, ref D[I + o_d]);
834: WORK[I - LL + 1 + o_work] = CS;
835: WORK[I - LL + 1 + NM1 + o_work] = SN;
836: WORK[I - LL + 1 + NM12 + o_work] = OLDCS;
837: WORK[I - LL + 1 + NM13 + o_work] = OLDSN;
838: }
839: H = D[M + o_d] * CS;
840: D[M + o_d] = H * OLDCS;
841: E[M - 1 + o_e] = H * OLDSN;
842: // *
843: // * Update singular vectors
844: // *
845: if (NCVT > 0)
846: {
847: this._dlasr.Run("L", "V", "F", M - LL + 1, NCVT, WORK, 1 + o_work
848: , WORK, N + o_work, ref VT, LL+1 * LDVT + o_vt, LDVT);
849: }
850: if (NRU > 0)
851: {
852: this._dlasr.Run("R", "V", "F", NRU, M - LL + 1, WORK, NM12 + 1 + o_work
853: , WORK, NM13 + 1 + o_work, ref U, 1+LL * LDU + o_u, LDU);
854: }
855: if (NCC > 0)
856: {
857: this._dlasr.Run("L", "V", "F", M - LL + 1, NCC, WORK, NM12 + 1 + o_work
858: , WORK, NM13 + 1 + o_work, ref C, LL+1 * LDC + o_c, LDC);
859: }
860: // *
861: // * Test convergence
862: // *
863: if (Math.Abs(E[M - 1 + o_e]) <= THRESH) E[M - 1 + o_e] = ZERO;
864: // *
865: }
866: else
867: {
868: // *
869: // * Chase bulge from bottom to top
870: // * Save cosines and sines for later singular vector updates
871: // *
872: CS = ONE;
873: OLDCS = ONE;
874: for (I = M; I >= LL + 1; I += - 1)
875: {
876: this._dlartg.Run(D[I + o_d] * CS, E[I - 1 + o_e], ref CS, ref SN, ref R);
877: if (I < M) E[I + o_e] = OLDSN * R;
878: this._dlartg.Run(OLDCS * R, D[I - 1 + o_d] * SN, ref OLDCS, ref OLDSN, ref D[I + o_d]);
879: WORK[I - LL + o_work] = CS;
880: WORK[I - LL + NM1 + o_work] = - SN;
881: WORK[I - LL + NM12 + o_work] = OLDCS;
882: WORK[I - LL + NM13 + o_work] = - OLDSN;
883: }
884: H = D[LL + o_d] * CS;
885: D[LL + o_d] = H * OLDCS;
886: E[LL + o_e] = H * OLDSN;
887: // *
888: // * Update singular vectors
889: // *
890: if (NCVT > 0)
891: {
892: this._dlasr.Run("L", "V", "B", M - LL + 1, NCVT, WORK, NM12 + 1 + o_work
893: , WORK, NM13 + 1 + o_work, ref VT, LL+1 * LDVT + o_vt, LDVT);
894: }
895: if (NRU > 0)
896: {
897: this._dlasr.Run("R", "V", "B", NRU, M - LL + 1, WORK, 1 + o_work
898: , WORK, N + o_work, ref U, 1+LL * LDU + o_u, LDU);
899: }
900: if (NCC > 0)
901: {
902: this._dlasr.Run("L", "V", "B", M - LL + 1, NCC, WORK, 1 + o_work
903: , WORK, N + o_work, ref C, LL+1 * LDC + o_c, LDC);
904: }
905: // *
906: // * Test convergence
907: // *
908: if (Math.Abs(E[LL + o_e]) <= THRESH) E[LL + o_e] = ZERO;
909: }
910: }
911: else
912: {
913: // *
914: // * Use nonzero shift
915: // *
916: if (IDIR == 1)
917: {
918: // *
919: // * Chase bulge from top to bottom
920: // * Save cosines and sines for later singular vector updates
921: // *
922: F = (Math.Abs(D[LL + o_d]) - SHIFT) * (FortranLib.Sign(ONE,D[LL + o_d]) + SHIFT / D[LL + o_d]);
923: G = E[LL + o_e];
924: for (I = LL; I <= M - 1; I++)
925: {
926: this._dlartg.Run(F, G, ref COSR, ref SINR, ref R);
927: if (I > LL) E[I - 1 + o_e] = R;
928: F = COSR * D[I + o_d] + SINR * E[I + o_e];
929: E[I + o_e] = COSR * E[I + o_e] - SINR * D[I + o_d];
930: G = SINR * D[I + 1 + o_d];
931: D[I + 1 + o_d] = COSR * D[I + 1 + o_d];
932: this._dlartg.Run(F, G, ref COSL, ref SINL, ref R);
933: D[I + o_d] = R;
934: F = COSL * E[I + o_e] + SINL * D[I + 1 + o_d];
935: D[I + 1 + o_d] = COSL * D[I + 1 + o_d] - SINL * E[I + o_e];
936: if (I < M - 1)
937: {
938: G = SINL * E[I + 1 + o_e];
939: E[I + 1 + o_e] = COSL * E[I + 1 + o_e];
940: }
941: WORK[I - LL + 1 + o_work] = COSR;
942: WORK[I - LL + 1 + NM1 + o_work] = SINR;
943: WORK[I - LL + 1 + NM12 + o_work] = COSL;
944: WORK[I - LL + 1 + NM13 + o_work] = SINL;
945: }
946: E[M - 1 + o_e] = F;
947: // *
948: // * Update singular vectors
949: // *
950: if (NCVT > 0)
951: {
952: this._dlasr.Run("L", "V", "F", M - LL + 1, NCVT, WORK, 1 + o_work
953: , WORK, N + o_work, ref VT, LL+1 * LDVT + o_vt, LDVT);
954: }
955: if (NRU > 0)
956: {
957: this._dlasr.Run("R", "V", "F", NRU, M - LL + 1, WORK, NM12 + 1 + o_work
958: , WORK, NM13 + 1 + o_work, ref U, 1+LL * LDU + o_u, LDU);
959: }
960: if (NCC > 0)
961: {
962: this._dlasr.Run("L", "V", "F", M - LL + 1, NCC, WORK, NM12 + 1 + o_work
963: , WORK, NM13 + 1 + o_work, ref C, LL+1 * LDC + o_c, LDC);
964: }
965: // *
966: // * Test convergence
967: // *
968: if (Math.Abs(E[M - 1 + o_e]) <= THRESH) E[M - 1 + o_e] = ZERO;
969: // *
970: }
971: else
972: {
973: // *
974: // * Chase bulge from bottom to top
975: // * Save cosines and sines for later singular vector updates
976: // *
977: F = (Math.Abs(D[M + o_d]) - SHIFT) * (FortranLib.Sign(ONE,D[M + o_d]) + SHIFT / D[M + o_d]);
978: G = E[M - 1 + o_e];
979: for (I = M; I >= LL + 1; I += - 1)
980: {
981: this._dlartg.Run(F, G, ref COSR, ref SINR, ref R);
982: if (I < M) E[I + o_e] = R;
983: F = COSR * D[I + o_d] + SINR * E[I - 1 + o_e];
984: E[I - 1 + o_e] = COSR * E[I - 1 + o_e] - SINR * D[I + o_d];
985: G = SINR * D[I - 1 + o_d];
986: D[I - 1 + o_d] = COSR * D[I - 1 + o_d];
987: this._dlartg.Run(F, G, ref COSL, ref SINL, ref R);
988: D[I + o_d] = R;
989: F = COSL * E[I - 1 + o_e] + SINL * D[I - 1 + o_d];
990: D[I - 1 + o_d] = COSL * D[I - 1 + o_d] - SINL * E[I - 1 + o_e];
991: if (I > LL + 1)
992: {
993: G = SINL * E[I - 2 + o_e];
994: E[I - 2 + o_e] = COSL * E[I - 2 + o_e];
995: }
996: WORK[I - LL + o_work] = COSR;
997: WORK[I - LL + NM1 + o_work] = - SINR;
998: WORK[I - LL + NM12 + o_work] = COSL;
999: WORK[I - LL + NM13 + o_work] = - SINL;
1000: }
1001: E[LL + o_e] = F;
1002: // *
1003: // * Test convergence
1004: // *
1005: if (Math.Abs(E[LL + o_e]) <= THRESH) E[LL + o_e] = ZERO;
1006: // *
1007: // * Update singular vectors if desired
1008: // *
1009: if (NCVT > 0)
1010: {
1011: this._dlasr.Run("L", "V", "B", M - LL + 1, NCVT, WORK, NM12 + 1 + o_work
1012: , WORK, NM13 + 1 + o_work, ref VT, LL+1 * LDVT + o_vt, LDVT);
1013: }
1014: if (NRU > 0)
1015: {
1016: this._dlasr.Run("R", "V", "B", NRU, M - LL + 1, WORK, 1 + o_work
1017: , WORK, N + o_work, ref U, 1+LL * LDU + o_u, LDU);
1018: }
1019: if (NCC > 0)
1020: {
1021: this._dlasr.Run("L", "V", "B", M - LL + 1, NCC, WORK, 1 + o_work
1022: , WORK, N + o_work, ref C, LL+1 * LDC + o_c, LDC);
1023: }
1024: }
1025: }
1026: // *
1027: // * QR iteration finished, go back and check convergence
1028: // *
1029: goto LABEL60;
1030: // *
1031: // * All singular values converged, so make them positive
1032: // *
1033: LABEL160:;
1034: for (I = 1; I <= N; I++)
1035: {
1036: if (D[I + o_d] < ZERO)
1037: {
1038: D[I + o_d] = - D[I + o_d];
1039: // *
1040: // * Change sign of singular vectors, if desired
1041: // *
1042: if (NCVT > 0) this._dscal.Run(NCVT, NEGONE, ref VT, I+1 * LDVT + o_vt, LDVT);
1043: }
1044: }
1045: // *
1046: // * Sort the singular values into decreasing order (insertion sort on
1047: // * singular values, but only one transposition per singular vector)
1048: // *
1049: for (I = 1; I <= N - 1; I++)
1050: {
1051: // *
1052: // * Scan for smallest D(I)
1053: // *
1054: ISUB = 1;
1055: SMIN = D[1 + o_d];
1056: for (J = 2; J <= N + 1 - I; J++)
1057: {
1058: if (D[J + o_d] <= SMIN)
1059: {
1060: ISUB = J;
1061: SMIN = D[J + o_d];
1062: }
1063: }
1064: if (ISUB != N + 1 - I)
1065: {
1066: // *
1067: // * Swap singular values and vectors
1068: // *
1069: D[ISUB + o_d] = D[N + 1 - I + o_d];
1070: D[N + 1 - I + o_d] = SMIN;
1071: if (NCVT > 0) this._dswap.Run(NCVT, ref VT, ISUB+1 * LDVT + o_vt, LDVT, ref VT, N + 1 - I+1 * LDVT + o_vt, LDVT);
1072: if (NRU > 0) this._dswap.Run(NRU, ref U, 1+ISUB * LDU + o_u, 1, ref U, 1+(N + 1 - I) * LDU + o_u, 1);
1073: if (NCC > 0) this._dswap.Run(NCC, ref C, ISUB+1 * LDC + o_c, LDC, ref C, N + 1 - I+1 * LDC + o_c, LDC);
1074: }
1075: }
1076: goto LABEL220;
1077: // *
1078: // * Maximum number of iterations exceeded, failure to converge
1079: // *
1080: LABEL200:;
1081: INFO = 0;
1082: for (I = 1; I <= N - 1; I++)
1083: {
1084: if (E[I + o_e] != ZERO) INFO = INFO + 1;
1085: }
1086: LABEL220:;
1087: return;
1088: // *
1089: // * End of DBDSQR
1090: // *
1091:
1092: #endregion
1093:
1094: }
1095: }
1096: }