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: /// DLASD3 finds all the square roots of the roots of the secular
27: /// equation, as defined by the values in D and Z. It makes the
28: /// appropriate calls to DLASD4 and then updates the singular
29: /// vectors by matrix multiplication.
30: ///
31: /// This code makes very mild assumptions about floating point
32: /// arithmetic. It will work on machines with a guard digit in
33: /// add/subtract, or on those binary machines without guard digits
34: /// which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
35: /// It could conceivably fail on hexadecimal or decimal machines
36: /// without guard digits, but we know of none.
37: ///
38: /// DLASD3 is called from DLASD1.
39: ///
40: ///</summary>
41: public class DLASD3
42: {
43:
44:
45: #region Dependencies
46:
47: DLAMC3 _dlamc3; DNRM2 _dnrm2; DCOPY _dcopy; DGEMM _dgemm; DLACPY _dlacpy; DLASCL _dlascl; DLASD4 _dlasd4; XERBLA _xerbla;
48:
49: #endregion
50:
51:
52: #region Fields
53:
54: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; const double NEGONE = - 1.0E+0; int CTEMP = 0; int I = 0;
55: int J = 0;int JC = 0; int KTEMP = 0; int M = 0; int N = 0; int NLP1 = 0; int NLP2 = 0; int NRP1 = 0; double RHO = 0;
56: double TEMP = 0;
57:
58: #endregion
59:
60: public DLASD3(DLAMC3 dlamc3, DNRM2 dnrm2, DCOPY dcopy, DGEMM dgemm, DLACPY dlacpy, DLASCL dlascl, DLASD4 dlasd4, XERBLA xerbla)
61: {
62:
63:
64: #region Set Dependencies
65:
66: this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy;
67: this._dlascl = dlascl;this._dlasd4 = dlasd4; this._xerbla = xerbla;
68:
69: #endregion
70:
71: }
72:
73: public DLASD3()
74: {
75:
76:
77: #region Dependencies (Initialization)
78:
79: DLAMC3 dlamc3 = new DLAMC3();
80: DNRM2 dnrm2 = new DNRM2();
81: DCOPY dcopy = new DCOPY();
82: LSAME lsame = new LSAME();
83: XERBLA xerbla = new XERBLA();
84: DLASD5 dlasd5 = new DLASD5();
85: DGEMM dgemm = new DGEMM(lsame, xerbla);
86: DLACPY dlacpy = new DLACPY(lsame);
87: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
88: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
89: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
90: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
91: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
92: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
93: DLAED6 dlaed6 = new DLAED6(dlamch);
94: DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
95:
96: #endregion
97:
98:
99: #region Set Dependencies
100:
101: this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy;
102: this._dlascl = dlascl;this._dlasd4 = dlasd4; this._xerbla = xerbla;
103:
104: #endregion
105:
106: }
107: /// <summary>
108: /// Purpose
109: /// =======
110: ///
111: /// DLASD3 finds all the square roots of the roots of the secular
112: /// equation, as defined by the values in D and Z. It makes the
113: /// appropriate calls to DLASD4 and then updates the singular
114: /// vectors by matrix multiplication.
115: ///
116: /// This code makes very mild assumptions about floating point
117: /// arithmetic. It will work on machines with a guard digit in
118: /// add/subtract, or on those binary machines without guard digits
119: /// which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
120: /// It could conceivably fail on hexadecimal or decimal machines
121: /// without guard digits, but we know of none.
122: ///
123: /// DLASD3 is called from DLASD1.
124: ///
125: ///</summary>
126: /// <param name="NL">
127: /// (input) INTEGER
128: /// The row dimension of the upper block. NL .GE. 1.
129: ///</param>
130: /// <param name="NR">
131: /// (input) INTEGER
132: /// The row dimension of the lower block. NR .GE. 1.
133: ///</param>
134: /// <param name="SQRE">
135: /// (input) INTEGER
136: /// = 0: the lower block is an NR-by-NR square matrix.
137: /// = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
138: ///
139: /// The bidiagonal matrix has N = NL + NR + 1 rows and
140: /// M = N + SQRE .GE. N columns.
141: ///</param>
142: /// <param name="K">
143: /// (input) INTEGER
144: /// The size of the secular equation, 1 =.LT. K = .LT. N.
145: ///</param>
146: /// <param name="D">
147: /// (output) DOUBLE PRECISION array, dimension(K)
148: /// On exit the square roots of the roots of the secular equation,
149: /// in ascending order.
150: ///</param>
151: /// <param name="Q">
152: /// (workspace) DOUBLE PRECISION array,
153: /// dimension at least (LDQ,K).
154: ///</param>
155: /// <param name="LDQ">
156: /// (input) INTEGER
157: /// The leading dimension of the array Q. LDQ .GE. K.
158: ///</param>
159: /// <param name="DSIGMA">
160: /// (input) DOUBLE PRECISION array, dimension(K)
161: /// The first K elements of this array contain the old roots
162: /// of the deflated updating problem. These are the poles
163: /// of the secular equation.
164: ///</param>
165: /// <param name="U">
166: /// (output) DOUBLE PRECISION array, dimension (LDU, N)
167: /// The last N - K columns of this matrix contain the deflated
168: /// left singular vectors.
169: ///</param>
170: /// <param name="LDU">
171: /// (input) INTEGER
172: /// The leading dimension of the array U. LDU .GE. N.
173: ///</param>
174: /// <param name="U2">
175: /// (input/output) DOUBLE PRECISION array, dimension (LDU2, N)
176: /// The first K columns of this matrix contain the non-deflated
177: /// left singular vectors for the split problem.
178: ///</param>
179: /// <param name="LDU2">
180: /// (input) INTEGER
181: /// The leading dimension of the array U2. LDU2 .GE. N.
182: ///</param>
183: /// <param name="VT">
184: /// (output) DOUBLE PRECISION array, dimension (LDVT, M)
185: /// The last M - K columns of VT' contain the deflated
186: /// right singular vectors.
187: ///</param>
188: /// <param name="LDVT">
189: /// (input) INTEGER
190: /// The leading dimension of the array VT. LDVT .GE. N.
191: ///</param>
192: /// <param name="VT2">
193: /// (input/output) DOUBLE PRECISION array, dimension (LDVT2, N)
194: /// The first K columns of VT2' contain the non-deflated
195: /// right singular vectors for the split problem.
196: ///</param>
197: /// <param name="LDVT2">
198: /// (input) INTEGER
199: /// The leading dimension of the array VT2. LDVT2 .GE. N.
200: ///</param>
201: /// <param name="IDXC">
202: /// (input) INTEGER array, dimension ( N )
203: /// The permutation used to arrange the columns of U (and rows of
204: /// VT) into three groups: the first group contains non-zero
205: /// entries only at and above (or before) NL +1; the second
206: /// contains non-zero entries only at and below (or after) NL+2;
207: /// and the third is dense. The first column of U and the row of
208: /// VT are treated separately, however.
209: ///
210: /// The rows of the singular vectors found by DLASD4
211: /// must be likewise permuted before the matrix multiplies can
212: /// take place.
213: ///</param>
214: /// <param name="CTOT">
215: /// (input) INTEGER array, dimension ( 4 )
216: /// A count of the total number of the various types of columns
217: /// in U (or rows in VT), as described in IDXC. The fourth column
218: /// type is any column which has been deflated.
219: ///</param>
220: /// <param name="Z">
221: /// (input) DOUBLE PRECISION array, dimension (K)
222: /// The first K elements of this array contain the components
223: /// of the deflation-adjusted updating row vector.
224: ///</param>
225: /// <param name="INFO">
226: /// (output) INTEGER
227: /// = 0: successful exit.
228: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
229: /// .GT. 0: if INFO = 1, an singular value did not converge
230: ///</param>
231: public void Run(int NL, int NR, int SQRE, int K, ref double[] D, int offset_d, ref double[] Q, int offset_q
232: , int LDQ, ref double[] DSIGMA, int offset_dsigma, ref double[] U, int offset_u, int LDU, double[] U2, int offset_u2, int LDU2
233: , ref double[] VT, int offset_vt, int LDVT, ref double[] VT2, int offset_vt2, int LDVT2, int[] IDXC, int offset_idxc, int[] CTOT, int offset_ctot
234: , ref double[] Z, int offset_z, ref int INFO)
235: {
236:
237: #region Array Index Correction
238:
239: int o_d = -1 + offset_d; int o_q = -1 - LDQ + offset_q; int o_dsigma = -1 + offset_dsigma;
240: int o_u = -1 - LDU + offset_u; int o_u2 = -1 - LDU2 + offset_u2; int o_vt = -1 - LDVT + offset_vt;
241: int o_vt2 = -1 - LDVT2 + offset_vt2; int o_idxc = -1 + offset_idxc; int o_ctot = -1 + offset_ctot;
242: int o_z = -1 + offset_z;
243:
244: #endregion
245:
246:
247: #region Prolog
248:
249: // *
250: // * -- LAPACK auxiliary routine (version 3.1) --
251: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
252: // * November 2006
253: // *
254: // * .. Scalar Arguments ..
255: // * ..
256: // * .. Array Arguments ..
257: // * ..
258: // *
259: // * Purpose
260: // * =======
261: // *
262: // * DLASD3 finds all the square roots of the roots of the secular
263: // * equation, as defined by the values in D and Z. It makes the
264: // * appropriate calls to DLASD4 and then updates the singular
265: // * vectors by matrix multiplication.
266: // *
267: // * This code makes very mild assumptions about floating point
268: // * arithmetic. It will work on machines with a guard digit in
269: // * add/subtract, or on those binary machines without guard digits
270: // * which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
271: // * It could conceivably fail on hexadecimal or decimal machines
272: // * without guard digits, but we know of none.
273: // *
274: // * DLASD3 is called from DLASD1.
275: // *
276: // * Arguments
277: // * =========
278: // *
279: // * NL (input) INTEGER
280: // * The row dimension of the upper block. NL >= 1.
281: // *
282: // * NR (input) INTEGER
283: // * The row dimension of the lower block. NR >= 1.
284: // *
285: // * SQRE (input) INTEGER
286: // * = 0: the lower block is an NR-by-NR square matrix.
287: // * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
288: // *
289: // * The bidiagonal matrix has N = NL + NR + 1 rows and
290: // * M = N + SQRE >= N columns.
291: // *
292: // * K (input) INTEGER
293: // * The size of the secular equation, 1 =< K = < N.
294: // *
295: // * D (output) DOUBLE PRECISION array, dimension(K)
296: // * On exit the square roots of the roots of the secular equation,
297: // * in ascending order.
298: // *
299: // * Q (workspace) DOUBLE PRECISION array,
300: // * dimension at least (LDQ,K).
301: // *
302: // * LDQ (input) INTEGER
303: // * The leading dimension of the array Q. LDQ >= K.
304: // *
305: // * DSIGMA (input) DOUBLE PRECISION array, dimension(K)
306: // * The first K elements of this array contain the old roots
307: // * of the deflated updating problem. These are the poles
308: // * of the secular equation.
309: // *
310: // * U (output) DOUBLE PRECISION array, dimension (LDU, N)
311: // * The last N - K columns of this matrix contain the deflated
312: // * left singular vectors.
313: // *
314: // * LDU (input) INTEGER
315: // * The leading dimension of the array U. LDU >= N.
316: // *
317: // * U2 (input/output) DOUBLE PRECISION array, dimension (LDU2, N)
318: // * The first K columns of this matrix contain the non-deflated
319: // * left singular vectors for the split problem.
320: // *
321: // * LDU2 (input) INTEGER
322: // * The leading dimension of the array U2. LDU2 >= N.
323: // *
324: // * VT (output) DOUBLE PRECISION array, dimension (LDVT, M)
325: // * The last M - K columns of VT' contain the deflated
326: // * right singular vectors.
327: // *
328: // * LDVT (input) INTEGER
329: // * The leading dimension of the array VT. LDVT >= N.
330: // *
331: // * VT2 (input/output) DOUBLE PRECISION array, dimension (LDVT2, N)
332: // * The first K columns of VT2' contain the non-deflated
333: // * right singular vectors for the split problem.
334: // *
335: // * LDVT2 (input) INTEGER
336: // * The leading dimension of the array VT2. LDVT2 >= N.
337: // *
338: // * IDXC (input) INTEGER array, dimension ( N )
339: // * The permutation used to arrange the columns of U (and rows of
340: // * VT) into three groups: the first group contains non-zero
341: // * entries only at and above (or before) NL +1; the second
342: // * contains non-zero entries only at and below (or after) NL+2;
343: // * and the third is dense. The first column of U and the row of
344: // * VT are treated separately, however.
345: // *
346: // * The rows of the singular vectors found by DLASD4
347: // * must be likewise permuted before the matrix multiplies can
348: // * take place.
349: // *
350: // * CTOT (input) INTEGER array, dimension ( 4 )
351: // * A count of the total number of the various types of columns
352: // * in U (or rows in VT), as described in IDXC. The fourth column
353: // * type is any column which has been deflated.
354: // *
355: // * Z (input) DOUBLE PRECISION array, dimension (K)
356: // * The first K elements of this array contain the components
357: // * of the deflation-adjusted updating row vector.
358: // *
359: // * INFO (output) INTEGER
360: // * = 0: successful exit.
361: // * < 0: if INFO = -i, the i-th argument had an illegal value.
362: // * > 0: if INFO = 1, an singular value did not converge
363: // *
364: // * Further Details
365: // * ===============
366: // *
367: // * Based on contributions by
368: // * Ming Gu and Huan Ren, Computer Science Division, University of
369: // * California at Berkeley, USA
370: // *
371: // * =====================================================================
372: // *
373: // * .. Parameters ..
374: // * ..
375: // * .. Local Scalars ..
376: // * ..
377: // * .. External Functions ..
378: // * ..
379: // * .. External Subroutines ..
380: // * ..
381: // * .. Intrinsic Functions ..
382: // INTRINSIC ABS, SIGN, SQRT;
383: // * ..
384: // * .. Executable Statements ..
385: // *
386: // * Test the input parameters.
387: // *
388:
389: #endregion
390:
391:
392: #region Body
393:
394: INFO = 0;
395: // *
396: if (NL < 1)
397: {
398: INFO = - 1;
399: }
400: else
401: {
402: if (NR < 1)
403: {
404: INFO = - 2;
405: }
406: else
407: {
408: if ((SQRE != 1) && (SQRE != 0))
409: {
410: INFO = - 3;
411: }
412: }
413: }
414: // *
415: N = NL + NR + 1;
416: M = N + SQRE;
417: NLP1 = NL + 1;
418: NLP2 = NL + 2;
419: // *
420: if ((K < 1) || (K > N))
421: {
422: INFO = - 4;
423: }
424: else
425: {
426: if (LDQ < K)
427: {
428: INFO = - 7;
429: }
430: else
431: {
432: if (LDU < N)
433: {
434: INFO = - 10;
435: }
436: else
437: {
438: if (LDU2 < N)
439: {
440: INFO = - 12;
441: }
442: else
443: {
444: if (LDVT < M)
445: {
446: INFO = - 14;
447: }
448: else
449: {
450: if (LDVT2 < M)
451: {
452: INFO = - 16;
453: }
454: }
455: }
456: }
457: }
458: }
459: if (INFO != 0)
460: {
461: this._xerbla.Run("DLASD3", - INFO);
462: return;
463: }
464: // *
465: // * Quick return if possible
466: // *
467: if (K == 1)
468: {
469: D[1 + o_d] = Math.Abs(Z[1 + o_z]);
470: this._dcopy.Run(M, VT2, 1+1 * LDVT2 + o_vt2, LDVT2, ref VT, 1+1 * LDVT + o_vt, LDVT);
471: if (Z[1 + o_z] > ZERO)
472: {
473: this._dcopy.Run(N, U2, 1+1 * LDU2 + o_u2, 1, ref U, 1+1 * LDU + o_u, 1);
474: }
475: else
476: {
477: for (I = 1; I <= N; I++)
478: {
479: U[I+1 * LDU + o_u] = - U2[I+1 * LDU2 + o_u2];
480: }
481: }
482: return;
483: }
484: // *
485: // * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
486: // * be computed with high relative accuracy (barring over/underflow).
487: // * This is a problem on machines without a guard digit in
488: // * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
489: // * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
490: // * which on any of these machines zeros out the bottommost
491: // * bit of DSIGMA(I) if it is 1; this makes the subsequent
492: // * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
493: // * occurs. On binary machines with a guard digit (almost all
494: // * machines) it does not change DSIGMA(I) at all. On hexadecimal
495: // * and decimal machines with a guard digit, it slightly
496: // * changes the bottommost bits of DSIGMA(I). It does not account
497: // * for hexadecimal or decimal machines without guard digits
498: // * (we know of none). We use a subroutine call to compute
499: // * 2*DSIGMA(I) to prevent optimizing compilers from eliminating
500: // * this code.
501: // *
502: for (I = 1; I <= K; I++)
503: {
504: DSIGMA[I + o_dsigma] = this._dlamc3.Run(DSIGMA[I + o_dsigma], DSIGMA[I + o_dsigma]) - DSIGMA[I + o_dsigma];
505: }
506: // *
507: // * Keep a copy of Z.
508: // *
509: this._dcopy.Run(K, Z, offset_z, 1, ref Q, offset_q, 1);
510: // *
511: // * Normalize Z.
512: // *
513: RHO = this._dnrm2.Run(K, Z, offset_z, 1);
514: this._dlascl.Run("G", 0, 0, RHO, ONE, K
515: , 1, ref Z, offset_z, K, ref INFO);
516: RHO = RHO * RHO;
517: // *
518: // * Find the new singular values.
519: // *
520: for (J = 1; J <= K; J++)
521: {
522: this._dlasd4.Run(K, J, DSIGMA, offset_dsigma, Z, offset_z, ref U, 1+J * LDU + o_u, RHO
523: , ref D[J + o_d], ref VT, 1+J * LDVT + o_vt, ref INFO);
524: // *
525: // * If the zero finder fails, the computation is terminated.
526: // *
527: if (INFO != 0)
528: {
529: return;
530: }
531: }
532: // *
533: // * Compute updated Z.
534: // *
535: for (I = 1; I <= K; I++)
536: {
537: Z[I + o_z] = U[I+K * LDU + o_u] * VT[I+K * LDVT + o_vt];
538: for (J = 1; J <= I - 1; J++)
539: {
540: Z[I + o_z] = Z[I + o_z] * (U[I+J * LDU + o_u] * VT[I+J * LDVT + o_vt] / (DSIGMA[I + o_dsigma] - DSIGMA[J + o_dsigma]) / (DSIGMA[I + o_dsigma] + DSIGMA[J + o_dsigma]));
541: }
542: for (J = I; J <= K - 1; J++)
543: {
544: Z[I + o_z] = Z[I + o_z] * (U[I+J * LDU + o_u] * VT[I+J * LDVT + o_vt] / (DSIGMA[I + o_dsigma] - DSIGMA[J + 1 + o_dsigma]) / (DSIGMA[I + o_dsigma] + DSIGMA[J + 1 + o_dsigma]));
545: }
546: Z[I + o_z] = FortranLib.Sign(Math.Sqrt(Math.Abs(Z[I + o_z])),Q[I+1 * LDQ + o_q]);
547: }
548: // *
549: // * Compute left singular vectors of the modified diagonal matrix,
550: // * and store related information for the right singular vectors.
551: // *
552: for (I = 1; I <= K; I++)
553: {
554: VT[1+I * LDVT + o_vt] = Z[1 + o_z] / U[1+I * LDU + o_u] / VT[1+I * LDVT + o_vt];
555: U[1+I * LDU + o_u] = NEGONE;
556: for (J = 2; J <= K; J++)
557: {
558: VT[J+I * LDVT + o_vt] = Z[J + o_z] / U[J+I * LDU + o_u] / VT[J+I * LDVT + o_vt];
559: U[J+I * LDU + o_u] = DSIGMA[J + o_dsigma] * VT[J+I * LDVT + o_vt];
560: }
561: TEMP = this._dnrm2.Run(K, U, 1+I * LDU + o_u, 1);
562: Q[1+I * LDQ + o_q] = U[1+I * LDU + o_u] / TEMP;
563: for (J = 2; J <= K; J++)
564: {
565: JC = IDXC[J + o_idxc];
566: Q[J+I * LDQ + o_q] = U[JC+I * LDU + o_u] / TEMP;
567: }
568: }
569: // *
570: // * Update the left singular vector matrix.
571: // *
572: if (K == 2)
573: {
574: this._dgemm.Run("N", "N", N, K, K, ONE
575: , U2, offset_u2, LDU2, Q, offset_q, LDQ, ZERO, ref U, offset_u
576: , LDU);
577: goto LABEL100;
578: }
579: if (CTOT[1 + o_ctot] > 0)
580: {
581: this._dgemm.Run("N", "N", NL, K, CTOT[1 + o_ctot], ONE
582: , U2, 1+2 * LDU2 + o_u2, LDU2, Q, 2+1 * LDQ + o_q, LDQ, ZERO, ref U, 1+1 * LDU + o_u
583: , LDU);
584: if (CTOT[3 + o_ctot] > 0)
585: {
586: KTEMP = 2 + CTOT[1 + o_ctot] + CTOT[2 + o_ctot];
587: this._dgemm.Run("N", "N", NL, K, CTOT[3 + o_ctot], ONE
588: , U2, 1+KTEMP * LDU2 + o_u2, LDU2, Q, KTEMP+1 * LDQ + o_q, LDQ, ONE, ref U, 1+1 * LDU + o_u
589: , LDU);
590: }
591: }
592: else
593: {
594: if (CTOT[3 + o_ctot] > 0)
595: {
596: KTEMP = 2 + CTOT[1 + o_ctot] + CTOT[2 + o_ctot];
597: this._dgemm.Run("N", "N", NL, K, CTOT[3 + o_ctot], ONE
598: , U2, 1+KTEMP * LDU2 + o_u2, LDU2, Q, KTEMP+1 * LDQ + o_q, LDQ, ZERO, ref U, 1+1 * LDU + o_u
599: , LDU);
600: }
601: else
602: {
603: this._dlacpy.Run("F", NL, K, U2, offset_u2, LDU2, ref U, offset_u
604: , LDU);
605: }
606: }
607: this._dcopy.Run(K, Q, 1+1 * LDQ + o_q, LDQ, ref U, NLP1+1 * LDU + o_u, LDU);
608: KTEMP = 2 + CTOT[1 + o_ctot];
609: CTEMP = CTOT[2 + o_ctot] + CTOT[3 + o_ctot];
610: this._dgemm.Run("N", "N", NR, K, CTEMP, ONE
611: , U2, NLP2+KTEMP * LDU2 + o_u2, LDU2, Q, KTEMP+1 * LDQ + o_q, LDQ, ZERO, ref U, NLP2+1 * LDU + o_u
612: , LDU);
613: // *
614: // * Generate the right singular vectors.
615: // *
616: LABEL100:;
617: for (I = 1; I <= K; I++)
618: {
619: TEMP = this._dnrm2.Run(K, VT, 1+I * LDVT + o_vt, 1);
620: Q[I+1 * LDQ + o_q] = VT[1+I * LDVT + o_vt] / TEMP;
621: for (J = 2; J <= K; J++)
622: {
623: JC = IDXC[J + o_idxc];
624: Q[I+J * LDQ + o_q] = VT[JC+I * LDVT + o_vt] / TEMP;
625: }
626: }
627: // *
628: // * Update the right singular vector matrix.
629: // *
630: if (K == 2)
631: {
632: this._dgemm.Run("N", "N", K, M, K, ONE
633: , Q, offset_q, LDQ, VT2, offset_vt2, LDVT2, ZERO, ref VT, offset_vt
634: , LDVT);
635: return;
636: }
637: KTEMP = 1 + CTOT[1 + o_ctot];
638: this._dgemm.Run("N", "N", K, NLP1, KTEMP, ONE
639: , Q, 1+1 * LDQ + o_q, LDQ, VT2, 1+1 * LDVT2 + o_vt2, LDVT2, ZERO, ref VT, 1+1 * LDVT + o_vt
640: , LDVT);
641: KTEMP = 2 + CTOT[1 + o_ctot] + CTOT[2 + o_ctot];
642: if (KTEMP <= LDVT2)
643: {
644: this._dgemm.Run("N", "N", K, NLP1, CTOT[3 + o_ctot], ONE
645: , Q, 1+KTEMP * LDQ + o_q, LDQ, VT2, KTEMP+1 * LDVT2 + o_vt2, LDVT2, ONE, ref VT, 1+1 * LDVT + o_vt
646: , LDVT);
647: }
648: // *
649: KTEMP = CTOT[1 + o_ctot] + 1;
650: NRP1 = NR + SQRE;
651: if (KTEMP > 1)
652: {
653: for (I = 1; I <= K; I++)
654: {
655: Q[I+KTEMP * LDQ + o_q] = Q[I+1 * LDQ + o_q];
656: }
657: for (I = NLP2; I <= M; I++)
658: {
659: VT2[KTEMP+I * LDVT2 + o_vt2] = VT2[1+I * LDVT2 + o_vt2];
660: }
661: }
662: CTEMP = 1 + CTOT[2 + o_ctot] + CTOT[3 + o_ctot];
663: this._dgemm.Run("N", "N", K, NRP1, CTEMP, ONE
664: , Q, 1+KTEMP * LDQ + o_q, LDQ, VT2, KTEMP+NLP2 * LDVT2 + o_vt2, LDVT2, ZERO, ref VT, 1+NLP2 * LDVT + o_vt
665: , LDVT);
666: // *
667: return;
668: // *
669: // * End of DLASD3
670: // *
671:
672: #endregion
673:
674: }
675: }
676: }