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