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: /// DLASD2 merges the two sets of singular values together into a single
27: /// sorted set. Then it tries to deflate the size of the problem.
28: /// There are two ways in which deflation can occur: when two or more
29: /// singular values are close together or if there is a tiny entry in the
30: /// Z vector. For each such occurrence the order of the related secular
31: /// equation problem is reduced by one.
32: ///
33: /// DLASD2 is called from DLASD1.
34: ///
35: ///</summary>
36: public class DLASD2
37: {
38:
39:
40: #region Dependencies
41:
42: DLAMCH _dlamch; DLAPY2 _dlapy2; DCOPY _dcopy; DLACPY _dlacpy; DLAMRG _dlamrg; DLASET _dlaset; DROT _drot; XERBLA _xerbla;
43:
44: #endregion
45:
46:
47: #region Fields
48:
49: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; const double TWO = 2.0E+0; const double EIGHT = 8.0E+0;
50: int[] CTOT = new int[4]; int o_ctot = -1;int[] PSM = new int[4]; int o_psm = -1;
51: int CT = 0;int I = 0; int IDXI = 0; int IDXJ = 0; int IDXJP = 0; int J = 0; int JP = 0; int JPREV = 0; int K2 = 0;
52: int M = 0;int N = 0; int NLP1 = 0; int NLP2 = 0; double C = 0; double EPS = 0; double HLFTOL = 0; double S = 0;
53: double TAU = 0;double TOL = 0; double Z1 = 0;
54:
55: #endregion
56:
57: public DLASD2(DLAMCH dlamch, DLAPY2 dlapy2, DCOPY dcopy, DLACPY dlacpy, DLAMRG dlamrg, DLASET dlaset, DROT drot, XERBLA xerbla)
58: {
59:
60:
61: #region Set Dependencies
62:
63: this._dlamch = dlamch; this._dlapy2 = dlapy2; this._dcopy = dcopy; this._dlacpy = dlacpy; this._dlamrg = dlamrg;
64: this._dlaset = dlaset;this._drot = drot; this._xerbla = xerbla;
65:
66: #endregion
67:
68: }
69:
70: public DLASD2()
71: {
72:
73:
74: #region Dependencies (Initialization)
75:
76: LSAME lsame = new LSAME();
77: DLAMC3 dlamc3 = new DLAMC3();
78: DLAPY2 dlapy2 = new DLAPY2();
79: DCOPY dcopy = new DCOPY();
80: DLAMRG dlamrg = new DLAMRG();
81: DROT drot = new DROT();
82: XERBLA xerbla = new XERBLA();
83: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
84: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
85: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
86: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
87: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
88: DLACPY dlacpy = new DLACPY(lsame);
89: DLASET dlaset = new DLASET(lsame);
90:
91: #endregion
92:
93:
94: #region Set Dependencies
95:
96: this._dlamch = dlamch; this._dlapy2 = dlapy2; this._dcopy = dcopy; this._dlacpy = dlacpy; this._dlamrg = dlamrg;
97: this._dlaset = dlaset;this._drot = drot; this._xerbla = xerbla;
98:
99: #endregion
100:
101: }
102: /// <summary>
103: /// Purpose
104: /// =======
105: ///
106: /// DLASD2 merges the two sets of singular values together into a single
107: /// sorted set. Then it tries to deflate the size of the problem.
108: /// There are two ways in which deflation can occur: when two or more
109: /// singular values are close together or if there is a tiny entry in the
110: /// Z vector. For each such occurrence the order of the related secular
111: /// equation problem is reduced by one.
112: ///
113: /// DLASD2 is called from DLASD1.
114: ///
115: ///</summary>
116: /// <param name="NL">
117: /// (input) INTEGER
118: /// The row dimension of the upper block. NL .GE. 1.
119: ///</param>
120: /// <param name="NR">
121: /// (input) INTEGER
122: /// The row dimension of the lower block. NR .GE. 1.
123: ///</param>
124: /// <param name="SQRE">
125: /// (input) INTEGER
126: /// = 0: the lower block is an NR-by-NR square matrix.
127: /// = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
128: ///
129: /// The bidiagonal matrix has N = NL + NR + 1 rows and
130: /// M = N + SQRE .GE. N columns.
131: ///</param>
132: /// <param name="K">
133: /// (output) INTEGER
134: /// Contains the dimension of the non-deflated matrix,
135: /// This is the order of the related secular equation. 1 .LE. K .LE.N.
136: ///</param>
137: /// <param name="D">
138: /// (input/output) DOUBLE PRECISION array, dimension(N)
139: /// On entry D contains the singular values of the two submatrices
140: /// to be combined. On exit D contains the trailing (N-K) updated
141: /// singular values (those which were deflated) sorted into
142: /// increasing order.
143: ///</param>
144: /// <param name="Z">
145: /// (output) DOUBLE PRECISION array, dimension(N)
146: /// On exit Z contains the updating row vector in the secular
147: /// equation.
148: ///</param>
149: /// <param name="ALPHA">
150: /// (input) DOUBLE PRECISION
151: /// Contains the diagonal element associated with the added row.
152: ///</param>
153: /// <param name="BETA">
154: /// (input) DOUBLE PRECISION
155: /// Contains the off-diagonal element associated with the added
156: /// row.
157: ///</param>
158: /// <param name="U">
159: /// (input/output) DOUBLE PRECISION array, dimension(LDU,N)
160: /// On entry U contains the left singular vectors of two
161: /// submatrices in the two square blocks with corners at (1,1),
162: /// (NL, NL), and (NL+2, NL+2), (N,N).
163: /// On exit U contains the trailing (N-K) updated left singular
164: /// vectors (those which were deflated) in its last N-K columns.
165: ///</param>
166: /// <param name="LDU">
167: /// (input) INTEGER
168: /// The leading dimension of the array U. LDU .GE. N.
169: ///</param>
170: /// <param name="VT">
171: /// (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
172: /// On entry VT' contains the right singular vectors of two
173: /// submatrices in the two square blocks with corners at (1,1),
174: /// (NL+1, NL+1), and (NL+2, NL+2), (M,M).
175: /// On exit VT' contains the trailing (N-K) updated right singular
176: /// vectors (those which were deflated) in its last N-K columns.
177: /// In case SQRE =1, the last row of VT spans the right null
178: /// space.
179: ///</param>
180: /// <param name="LDVT">
181: /// (input) INTEGER
182: /// The leading dimension of the array VT. LDVT .GE. M.
183: ///</param>
184: /// <param name="DSIGMA">
185: /// (output) DOUBLE PRECISION array, dimension (N)
186: /// Contains a copy of the diagonal elements (K-1 singular values
187: /// and one zero) in the secular equation.
188: ///</param>
189: /// <param name="U2">
190: /// (output) DOUBLE PRECISION array, dimension(LDU2,N)
191: /// Contains a copy of the first K-1 left singular vectors which
192: /// will be used by DLASD3 in a matrix multiply (DGEMM) to solve
193: /// for the new left singular vectors. U2 is arranged into four
194: /// blocks. The first block contains a column with 1 at NL+1 and
195: /// zero everywhere else; the second block contains non-zero
196: /// entries only at and above NL; the third contains non-zero
197: /// entries only below NL+1; and the fourth is dense.
198: ///</param>
199: /// <param name="LDU2">
200: /// (input) INTEGER
201: /// The leading dimension of the array U2. LDU2 .GE. N.
202: ///</param>
203: /// <param name="VT2">
204: /// (output) DOUBLE PRECISION array, dimension(LDVT2,N)
205: /// VT2' contains a copy of the first K right singular vectors
206: /// which will be used by DLASD3 in a matrix multiply (DGEMM) to
207: /// solve for the new right singular vectors. VT2 is arranged into
208: /// three blocks. The first block contains a row that corresponds
209: /// to the special 0 diagonal element in SIGMA; the second block
210: /// contains non-zeros only at and before NL +1; the third block
211: /// contains non-zeros only at and after NL +2.
212: ///</param>
213: /// <param name="LDVT2">
214: /// (input) INTEGER
215: /// The leading dimension of the array VT2. LDVT2 .GE. M.
216: ///</param>
217: /// <param name="IDXP">
218: /// (workspace) INTEGER array dimension(N)
219: /// This will contain the permutation used to place deflated
220: /// values of D at the end of the array. On output IDXP(2:K)
221: /// points to the nondeflated D-values and IDXP(K+1:N)
222: /// points to the deflated singular values.
223: ///</param>
224: /// <param name="IDX">
225: /// (workspace) INTEGER array dimension(N)
226: /// This will contain the permutation used to sort the contents of
227: /// D into ascending order.
228: ///</param>
229: /// <param name="IDXC">
230: /// (output) INTEGER array dimension(N)
231: /// This will contain the permutation used to arrange the columns
232: /// of the deflated U matrix into three groups: the first group
233: /// contains non-zero entries only at and above NL, the second
234: /// contains non-zero entries only below NL+2, and the third is
235: /// dense.
236: ///</param>
237: /// <param name="IDXQ">
238: /// (input/output) INTEGER array dimension(N)
239: /// This contains the permutation which separately sorts the two
240: /// sub-problems in D into ascending order. Note that entries in
241: /// the first hlaf of this permutation must first be moved one
242: /// position backward; and entries in the second half
243: /// must first have NL+1 added to their values.
244: ///</param>
245: /// <param name="COLTYP">
246: /// (workspace/output) INTEGER array dimension(N)
247: /// As workspace, this will contain a label which will indicate
248: /// which of the following types a column in the U2 matrix or a
249: /// row in the VT2 matrix is:
250: /// 1 : non-zero in the upper half only
251: /// 2 : non-zero in the lower half only
252: /// 3 : dense
253: /// 4 : deflated
254: ///
255: /// On exit, it is an array of dimension 4, with COLTYP(I) being
256: /// the dimension of the I-th type columns.
257: ///</param>
258: /// <param name="INFO">
259: /// (output) INTEGER
260: /// = 0: successful exit.
261: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
262: ///</param>
263: public void Run(int NL, int NR, int SQRE, ref int K, ref double[] D, int offset_d, ref double[] Z, int offset_z
264: , double ALPHA, double BETA, ref double[] U, int offset_u, int LDU, ref double[] VT, int offset_vt, int LDVT
265: , ref double[] DSIGMA, int offset_dsigma, ref double[] U2, int offset_u2, int LDU2, ref double[] VT2, int offset_vt2, int LDVT2, ref int[] IDXP, int offset_idxp
266: , ref int[] IDX, int offset_idx, ref int[] IDXC, int offset_idxc, ref int[] IDXQ, int offset_idxq, ref int[] COLTYP, int offset_coltyp, ref int INFO)
267: {
268:
269: #region Array Index Correction
270:
271: int o_d = -1 + offset_d; int o_z = -1 + offset_z; int o_u = -1 - LDU + offset_u; int o_vt = -1 - LDVT + offset_vt;
272: int o_dsigma = -1 + offset_dsigma; int o_u2 = -1 - LDU2 + offset_u2; int o_vt2 = -1 - LDVT2 + offset_vt2;
273: int o_idxp = -1 + offset_idxp; int o_idx = -1 + offset_idx; int o_idxc = -1 + offset_idxc;
274: int o_idxq = -1 + offset_idxq; int o_coltyp = -1 + offset_coltyp;
275:
276: #endregion
277:
278:
279: #region Prolog
280:
281: // *
282: // * -- LAPACK auxiliary routine (version 3.1) --
283: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
284: // * November 2006
285: // *
286: // * .. Scalar Arguments ..
287: // * ..
288: // * .. Array Arguments ..
289: // * ..
290: // *
291: // * Purpose
292: // * =======
293: // *
294: // * DLASD2 merges the two sets of singular values together into a single
295: // * sorted set. Then it tries to deflate the size of the problem.
296: // * There are two ways in which deflation can occur: when two or more
297: // * singular values are close together or if there is a tiny entry in the
298: // * Z vector. For each such occurrence the order of the related secular
299: // * equation problem is reduced by one.
300: // *
301: // * DLASD2 is called from DLASD1.
302: // *
303: // * Arguments
304: // * =========
305: // *
306: // * NL (input) INTEGER
307: // * The row dimension of the upper block. NL >= 1.
308: // *
309: // * NR (input) INTEGER
310: // * The row dimension of the lower block. NR >= 1.
311: // *
312: // * SQRE (input) INTEGER
313: // * = 0: the lower block is an NR-by-NR square matrix.
314: // * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
315: // *
316: // * The bidiagonal matrix has N = NL + NR + 1 rows and
317: // * M = N + SQRE >= N columns.
318: // *
319: // * K (output) INTEGER
320: // * Contains the dimension of the non-deflated matrix,
321: // * This is the order of the related secular equation. 1 <= K <=N.
322: // *
323: // * D (input/output) DOUBLE PRECISION array, dimension(N)
324: // * On entry D contains the singular values of the two submatrices
325: // * to be combined. On exit D contains the trailing (N-K) updated
326: // * singular values (those which were deflated) sorted into
327: // * increasing order.
328: // *
329: // * Z (output) DOUBLE PRECISION array, dimension(N)
330: // * On exit Z contains the updating row vector in the secular
331: // * equation.
332: // *
333: // * ALPHA (input) DOUBLE PRECISION
334: // * Contains the diagonal element associated with the added row.
335: // *
336: // * BETA (input) DOUBLE PRECISION
337: // * Contains the off-diagonal element associated with the added
338: // * row.
339: // *
340: // * U (input/output) DOUBLE PRECISION array, dimension(LDU,N)
341: // * On entry U contains the left singular vectors of two
342: // * submatrices in the two square blocks with corners at (1,1),
343: // * (NL, NL), and (NL+2, NL+2), (N,N).
344: // * On exit U contains the trailing (N-K) updated left singular
345: // * vectors (those which were deflated) in its last N-K columns.
346: // *
347: // * LDU (input) INTEGER
348: // * The leading dimension of the array U. LDU >= N.
349: // *
350: // * VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
351: // * On entry VT' contains the right singular vectors of two
352: // * submatrices in the two square blocks with corners at (1,1),
353: // * (NL+1, NL+1), and (NL+2, NL+2), (M,M).
354: // * On exit VT' contains the trailing (N-K) updated right singular
355: // * vectors (those which were deflated) in its last N-K columns.
356: // * In case SQRE =1, the last row of VT spans the right null
357: // * space.
358: // *
359: // * LDVT (input) INTEGER
360: // * The leading dimension of the array VT. LDVT >= M.
361: // *
362: // * DSIGMA (output) DOUBLE PRECISION array, dimension (N)
363: // * Contains a copy of the diagonal elements (K-1 singular values
364: // * and one zero) in the secular equation.
365: // *
366: // * U2 (output) DOUBLE PRECISION array, dimension(LDU2,N)
367: // * Contains a copy of the first K-1 left singular vectors which
368: // * will be used by DLASD3 in a matrix multiply (DGEMM) to solve
369: // * for the new left singular vectors. U2 is arranged into four
370: // * blocks. The first block contains a column with 1 at NL+1 and
371: // * zero everywhere else; the second block contains non-zero
372: // * entries only at and above NL; the third contains non-zero
373: // * entries only below NL+1; and the fourth is dense.
374: // *
375: // * LDU2 (input) INTEGER
376: // * The leading dimension of the array U2. LDU2 >= N.
377: // *
378: // * VT2 (output) DOUBLE PRECISION array, dimension(LDVT2,N)
379: // * VT2' contains a copy of the first K right singular vectors
380: // * which will be used by DLASD3 in a matrix multiply (DGEMM) to
381: // * solve for the new right singular vectors. VT2 is arranged into
382: // * three blocks. The first block contains a row that corresponds
383: // * to the special 0 diagonal element in SIGMA; the second block
384: // * contains non-zeros only at and before NL +1; the third block
385: // * contains non-zeros only at and after NL +2.
386: // *
387: // * LDVT2 (input) INTEGER
388: // * The leading dimension of the array VT2. LDVT2 >= M.
389: // *
390: // * IDXP (workspace) INTEGER array dimension(N)
391: // * This will contain the permutation used to place deflated
392: // * values of D at the end of the array. On output IDXP(2:K)
393: // * points to the nondeflated D-values and IDXP(K+1:N)
394: // * points to the deflated singular values.
395: // *
396: // * IDX (workspace) INTEGER array dimension(N)
397: // * This will contain the permutation used to sort the contents of
398: // * D into ascending order.
399: // *
400: // * IDXC (output) INTEGER array dimension(N)
401: // * This will contain the permutation used to arrange the columns
402: // * of the deflated U matrix into three groups: the first group
403: // * contains non-zero entries only at and above NL, the second
404: // * contains non-zero entries only below NL+2, and the third is
405: // * dense.
406: // *
407: // * IDXQ (input/output) INTEGER array dimension(N)
408: // * This contains the permutation which separately sorts the two
409: // * sub-problems in D into ascending order. Note that entries in
410: // * the first hlaf of this permutation must first be moved one
411: // * position backward; and entries in the second half
412: // * must first have NL+1 added to their values.
413: // *
414: // * COLTYP (workspace/output) INTEGER array dimension(N)
415: // * As workspace, this will contain a label which will indicate
416: // * which of the following types a column in the U2 matrix or a
417: // * row in the VT2 matrix is:
418: // * 1 : non-zero in the upper half only
419: // * 2 : non-zero in the lower half only
420: // * 3 : dense
421: // * 4 : deflated
422: // *
423: // * On exit, it is an array of dimension 4, with COLTYP(I) being
424: // * the dimension of the I-th type columns.
425: // *
426: // * INFO (output) INTEGER
427: // * = 0: successful exit.
428: // * < 0: if INFO = -i, the i-th argument had an illegal value.
429: // *
430: // * Further Details
431: // * ===============
432: // *
433: // * Based on contributions by
434: // * Ming Gu and Huan Ren, Computer Science Division, University of
435: // * California at Berkeley, USA
436: // *
437: // * =====================================================================
438: // *
439: // * .. Parameters ..
440: // * ..
441: // * .. Local Arrays ..
442: // * ..
443: // * .. Local Scalars ..
444: // * ..
445: // * .. External Functions ..
446: // * ..
447: // * .. External Subroutines ..
448: // * ..
449: // * .. Intrinsic Functions ..
450: // INTRINSIC ABS, MAX;
451: // * ..
452: // * .. Executable Statements ..
453: // *
454: // * Test the input parameters.
455: // *
456:
457: #endregion
458:
459:
460: #region Body
461:
462: INFO = 0;
463: // *
464: if (NL < 1)
465: {
466: INFO = - 1;
467: }
468: else
469: {
470: if (NR < 1)
471: {
472: INFO = - 2;
473: }
474: else
475: {
476: if ((SQRE != 1) && (SQRE != 0))
477: {
478: INFO = - 3;
479: }
480: }
481: }
482: // *
483: N = NL + NR + 1;
484: M = N + SQRE;
485: // *
486: if (LDU < N)
487: {
488: INFO = - 10;
489: }
490: else
491: {
492: if (LDVT < M)
493: {
494: INFO = - 12;
495: }
496: else
497: {
498: if (LDU2 < N)
499: {
500: INFO = - 15;
501: }
502: else
503: {
504: if (LDVT2 < M)
505: {
506: INFO = - 17;
507: }
508: }
509: }
510: }
511: if (INFO != 0)
512: {
513: this._xerbla.Run("DLASD2", - INFO);
514: return;
515: }
516: // *
517: NLP1 = NL + 1;
518: NLP2 = NL + 2;
519: // *
520: // * Generate the first part of the vector Z; and move the singular
521: // * values in the first part of D one position backward.
522: // *
523: Z1 = ALPHA * VT[NLP1+NLP1 * LDVT + o_vt];
524: Z[1 + o_z] = Z1;
525: for (I = NL; I >= 1; I += - 1)
526: {
527: Z[I + 1 + o_z] = ALPHA * VT[I+NLP1 * LDVT + o_vt];
528: D[I + 1 + o_d] = D[I + o_d];
529: IDXQ[I + 1 + o_idxq] = IDXQ[I + o_idxq] + 1;
530: }
531: // *
532: // * Generate the second part of the vector Z.
533: // *
534: for (I = NLP2; I <= M; I++)
535: {
536: Z[I + o_z] = BETA * VT[I+NLP2 * LDVT + o_vt];
537: }
538: // *
539: // * Initialize some reference arrays.
540: // *
541: for (I = 2; I <= NLP1; I++)
542: {
543: COLTYP[I + o_coltyp] = 1;
544: }
545: for (I = NLP2; I <= N; I++)
546: {
547: COLTYP[I + o_coltyp] = 2;
548: }
549: // *
550: // * Sort the singular values into increasing order
551: // *
552: for (I = NLP2; I <= N; I++)
553: {
554: IDXQ[I + o_idxq] = IDXQ[I + o_idxq] + NLP1;
555: }
556: // *
557: // * DSIGMA, IDXC, IDXC, and the first column of U2
558: // * are used as storage space.
559: // *
560: for (I = 2; I <= N; I++)
561: {
562: DSIGMA[I + o_dsigma] = D[IDXQ[I + o_idxq] + o_d];
563: U2[I+1 * LDU2 + o_u2] = Z[IDXQ[I + o_idxq] + o_z];
564: IDXC[I + o_idxc] = COLTYP[IDXQ[I + o_idxq] + o_coltyp];
565: }
566: // *
567: this._dlamrg.Run(NL, NR, DSIGMA, 2 + o_dsigma, 1, 1, ref IDX, 2 + o_idx);
568: // *
569: for (I = 2; I <= N; I++)
570: {
571: IDXI = 1 + IDX[I + o_idx];
572: D[I + o_d] = DSIGMA[IDXI + o_dsigma];
573: Z[I + o_z] = U2[IDXI+1 * LDU2 + o_u2];
574: COLTYP[I + o_coltyp] = IDXC[IDXI + o_idxc];
575: }
576: // *
577: // * Calculate the allowable deflation tolerance
578: // *
579: EPS = this._dlamch.Run("Epsilon");
580: TOL = Math.Max(Math.Abs(ALPHA), Math.Abs(BETA));
581: TOL = EIGHT * EPS * Math.Max(Math.Abs(D[N + o_d]), TOL);
582: // *
583: // * There are 2 kinds of deflation -- first a value in the z-vector
584: // * is small, second two (or more) singular values are very close
585: // * together (their difference is small).
586: // *
587: // * If the value in the z-vector is small, we simply permute the
588: // * array so that the corresponding singular value is moved to the
589: // * end.
590: // *
591: // * If two values in the D-vector are close, we perform a two-sided
592: // * rotation designed to make one of the corresponding z-vector
593: // * entries zero, and then permute the array so that the deflated
594: // * singular value is moved to the end.
595: // *
596: // * If there are multiple singular values then the problem deflates.
597: // * Here the number of equal singular values are found. As each equal
598: // * singular value is found, an elementary reflector is computed to
599: // * rotate the corresponding singular subspace so that the
600: // * corresponding components of Z are zero in this new basis.
601: // *
602: K = 1;
603: K2 = N + 1;
604: for (J = 2; J <= N; J++)
605: {
606: if (Math.Abs(Z[J + o_z]) <= TOL)
607: {
608: // *
609: // * Deflate due to small z component.
610: // *
611: K2 = K2 - 1;
612: IDXP[K2 + o_idxp] = J;
613: COLTYP[J + o_coltyp] = 4;
614: if (J == N) goto LABEL120;
615: }
616: else
617: {
618: JPREV = J;
619: goto LABEL90;
620: }
621: }
622: LABEL90:;
623: J = JPREV;
624: LABEL100:;
625: J = J + 1;
626: if (J > N) goto LABEL110;
627: if (Math.Abs(Z[J + o_z]) <= TOL)
628: {
629: // *
630: // * Deflate due to small z component.
631: // *
632: K2 = K2 - 1;
633: IDXP[K2 + o_idxp] = J;
634: COLTYP[J + o_coltyp] = 4;
635: }
636: else
637: {
638: // *
639: // * Check if singular values are close enough to allow deflation.
640: // *
641: if (Math.Abs(D[J + o_d] - D[JPREV + o_d]) <= TOL)
642: {
643: // *
644: // * Deflation is possible.
645: // *
646: S = Z[JPREV + o_z];
647: C = Z[J + o_z];
648: // *
649: // * Find sqrt(a**2+b**2) without overflow or
650: // * destructive underflow.
651: // *
652: TAU = this._dlapy2.Run(C, S);
653: C = C / TAU;
654: S = - S / TAU;
655: Z[J + o_z] = TAU;
656: Z[JPREV + o_z] = ZERO;
657: // *
658: // * Apply back the Givens rotation to the left and right
659: // * singular vector matrices.
660: // *
661: IDXJP = IDXQ[IDX[JPREV + o_idx] + 1 + o_idxq];
662: IDXJ = IDXQ[IDX[J + o_idx] + 1 + o_idxq];
663: if (IDXJP <= NLP1)
664: {
665: IDXJP = IDXJP - 1;
666: }
667: if (IDXJ <= NLP1)
668: {
669: IDXJ = IDXJ - 1;
670: }
671: this._drot.Run(N, ref U, 1+IDXJP * LDU + o_u, 1, ref U, 1+IDXJ * LDU + o_u, 1, C
672: , S);
673: this._drot.Run(M, ref VT, IDXJP+1 * LDVT + o_vt, LDVT, ref VT, IDXJ+1 * LDVT + o_vt, LDVT, C
674: , S);
675: if (COLTYP[J + o_coltyp] != COLTYP[JPREV + o_coltyp])
676: {
677: COLTYP[J + o_coltyp] = 3;
678: }
679: COLTYP[JPREV + o_coltyp] = 4;
680: K2 = K2 - 1;
681: IDXP[K2 + o_idxp] = JPREV;
682: JPREV = J;
683: }
684: else
685: {
686: K = K + 1;
687: U2[K+1 * LDU2 + o_u2] = Z[JPREV + o_z];
688: DSIGMA[K + o_dsigma] = D[JPREV + o_d];
689: IDXP[K + o_idxp] = JPREV;
690: JPREV = J;
691: }
692: }
693: goto LABEL100;
694: LABEL110:;
695: // *
696: // * Record the last singular value.
697: // *
698: K = K + 1;
699: U2[K+1 * LDU2 + o_u2] = Z[JPREV + o_z];
700: DSIGMA[K + o_dsigma] = D[JPREV + o_d];
701: IDXP[K + o_idxp] = JPREV;
702: // *
703: LABEL120:;
704: // *
705: // * Count up the total number of the various types of columns, then
706: // * form a permutation which positions the four column types into
707: // * four groups of uniform structure (although one or more of these
708: // * groups may be empty).
709: // *
710: for (J = 1; J <= 4; J++)
711: {
712: CTOT[J + o_ctot] = 0;
713: }
714: for (J = 2; J <= N; J++)
715: {
716: CT = COLTYP[J + o_coltyp];
717: CTOT[CT + o_ctot] = CTOT[CT + o_ctot] + 1;
718: }
719: // *
720: // * PSM(*) = Position in SubMatrix (of types 1 through 4)
721: // *
722: PSM[1 + o_psm] = 2;
723: PSM[2 + o_psm] = 2 + CTOT[1 + o_ctot];
724: PSM[3 + o_psm] = PSM[2 + o_psm] + CTOT[2 + o_ctot];
725: PSM[4 + o_psm] = PSM[3 + o_psm] + CTOT[3 + o_ctot];
726: // *
727: // * Fill out the IDXC array so that the permutation which it induces
728: // * will place all type-1 columns first, all type-2 columns next,
729: // * then all type-3's, and finally all type-4's, starting from the
730: // * second column. This applies similarly to the rows of VT.
731: // *
732: for (J = 2; J <= N; J++)
733: {
734: JP = IDXP[J + o_idxp];
735: CT = COLTYP[JP + o_coltyp];
736: IDXC[PSM[CT + o_psm] + o_idxc] = J;
737: PSM[CT + o_psm] = PSM[CT + o_psm] + 1;
738: }
739: // *
740: // * Sort the singular values and corresponding singular vectors into
741: // * DSIGMA, U2, and VT2 respectively. The singular values/vectors
742: // * which were not deflated go into the first K slots of DSIGMA, U2,
743: // * and VT2 respectively, while those which were deflated go into the
744: // * last N - K slots, except that the first column/row will be treated
745: // * separately.
746: // *
747: for (J = 2; J <= N; J++)
748: {
749: JP = IDXP[J + o_idxp];
750: DSIGMA[J + o_dsigma] = D[JP + o_d];
751: IDXJ = IDXQ[IDX[IDXP[IDXC[J + o_idxc] + o_idxp] + o_idx] + 1 + o_idxq];
752: if (IDXJ <= NLP1)
753: {
754: IDXJ = IDXJ - 1;
755: }
756: this._dcopy.Run(N, U, 1+IDXJ * LDU + o_u, 1, ref U2, 1+J * LDU2 + o_u2, 1);
757: this._dcopy.Run(M, VT, IDXJ+1 * LDVT + o_vt, LDVT, ref VT2, J+1 * LDVT2 + o_vt2, LDVT2);
758: }
759: // *
760: // * Determine DSIGMA(1), DSIGMA(2) and Z(1)
761: // *
762: DSIGMA[1 + o_dsigma] = ZERO;
763: HLFTOL = TOL / TWO;
764: if (Math.Abs(DSIGMA[2 + o_dsigma]) <= HLFTOL) DSIGMA[2 + o_dsigma] = HLFTOL;
765: if (M > N)
766: {
767: Z[1 + o_z] = this._dlapy2.Run(Z1, Z[M + o_z]);
768: if (Z[1 + o_z] <= TOL)
769: {
770: C = ONE;
771: S = ZERO;
772: Z[1 + o_z] = TOL;
773: }
774: else
775: {
776: C = Z1 / Z[1 + o_z];
777: S = Z[M + o_z] / Z[1 + o_z];
778: }
779: }
780: else
781: {
782: if (Math.Abs(Z1) <= TOL)
783: {
784: Z[1 + o_z] = TOL;
785: }
786: else
787: {
788: Z[1 + o_z] = Z1;
789: }
790: }
791: // *
792: // * Move the rest of the updating row to Z.
793: // *
794: this._dcopy.Run(K - 1, U2, 2+1 * LDU2 + o_u2, 1, ref Z, 2 + o_z, 1);
795: // *
796: // * Determine the first column of U2, the first row of VT2 and the
797: // * last row of VT.
798: // *
799: this._dlaset.Run("A", N, 1, ZERO, ZERO, ref U2, offset_u2
800: , LDU2);
801: U2[NLP1+1 * LDU2 + o_u2] = ONE;
802: if (M > N)
803: {
804: for (I = 1; I <= NLP1; I++)
805: {
806: VT[M+I * LDVT + o_vt] = - S * VT[NLP1+I * LDVT + o_vt];
807: VT2[1+I * LDVT2 + o_vt2] = C * VT[NLP1+I * LDVT + o_vt];
808: }
809: for (I = NLP2; I <= M; I++)
810: {
811: VT2[1+I * LDVT2 + o_vt2] = S * VT[M+I * LDVT + o_vt];
812: VT[M+I * LDVT + o_vt] = C * VT[M+I * LDVT + o_vt];
813: }
814: }
815: else
816: {
817: this._dcopy.Run(M, VT, NLP1+1 * LDVT + o_vt, LDVT, ref VT2, 1+1 * LDVT2 + o_vt2, LDVT2);
818: }
819: if (M > N)
820: {
821: this._dcopy.Run(M, VT, M+1 * LDVT + o_vt, LDVT, ref VT2, M+1 * LDVT2 + o_vt2, LDVT2);
822: }
823: // *
824: // * The deflated singular values and their corresponding vectors go
825: // * into the back of D, U, and V respectively.
826: // *
827: if (N > K)
828: {
829: this._dcopy.Run(N - K, DSIGMA, K + 1 + o_dsigma, 1, ref D, K + 1 + o_d, 1);
830: this._dlacpy.Run("A", N, N - K, U2, 1+(K + 1) * LDU2 + o_u2, LDU2, ref U, 1+(K + 1) * LDU + o_u
831: , LDU);
832: this._dlacpy.Run("A", N - K, M, VT2, K + 1+1 * LDVT2 + o_vt2, LDVT2, ref VT, K + 1+1 * LDVT + o_vt
833: , LDVT);
834: }
835: // *
836: // * Copy CTOT into COLTYP for referencing in DLASD3.
837: // *
838: for (J = 1; J <= 4; J++)
839: {
840: COLTYP[J + o_coltyp] = CTOT[J + o_ctot];
841: }
842: // *
843: return;
844: // *
845: // * End of DLASD2
846: // *
847:
848: #endregion
849:
850: }
851: }
852: }