1: #region Translated by Jose Antonio De Santiago-Castillo.
2:
3: //Translated by Jose Antonio De Santiago-Castillo.
4: //E-mail:JAntonioDeSantiago@gmail.com
5: //Web: www.DotNumerics.com
6: //
7: //Fortran to C# Translation.
8: //Translated by:
9: //F2CSharp Version 0.71 (November 10, 2009)
10: //Code Optimizations: None
11: //
12: #endregion
13:
14: using System;
15: using DotNumerics.FortranLibrary;
16:
17: namespace DotNumerics.CSLapack
18: {
19: /// <summary>
20: /// -- LAPACK routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DLALS0 applies back the multiplying factors of either the left or the
27: /// right singular vector matrix of a diagonal matrix appended by a row
28: /// to the right hand side matrix B in solving the least squares problem
29: /// using the divide-and-conquer SVD approach.
30: ///
31: /// For the left singular vector matrix, three types of orthogonal
32: /// matrices are involved:
33: ///
34: /// (1L) Givens rotations: the number of such rotations is GIVPTR; the
35: /// pairs of columns/rows they were applied to are stored in GIVCOL;
36: /// and the C- and S-values of these rotations are stored in GIVNUM.
37: ///
38: /// (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
39: /// row, and for J=2:N, PERM(J)-th row of B is to be moved to the
40: /// J-th row.
41: ///
42: /// (3L) The left singular vector matrix of the remaining matrix.
43: ///
44: /// For the right singular vector matrix, four types of orthogonal
45: /// matrices are involved:
46: ///
47: /// (1R) The right singular vector matrix of the remaining matrix.
48: ///
49: /// (2R) If SQRE = 1, one extra Givens rotation to generate the right
50: /// null space.
51: ///
52: /// (3R) The inverse transformation of (2L).
53: ///
54: /// (4R) The inverse transformation of (1L).
55: ///
56: ///</summary>
57: public class DLALS0
58: {
59:
60:
61: #region Dependencies
62:
63: DCOPY _dcopy; DGEMV _dgemv; DLACPY _dlacpy; DLASCL _dlascl; DROT _drot; DSCAL _dscal; XERBLA _xerbla; DLAMC3 _dlamc3;
64: DNRM2 _dnrm2;
65:
66: #endregion
67:
68:
69: #region Fields
70:
71: const double ONE = 1.0E0; const double ZERO = 0.0E0; const double NEGONE = - 1.0E0; int I = 0; int J = 0; int M = 0;
72: int N = 0;int NLP1 = 0; double DIFLJ = 0; double DIFRJ = 0; double DJ = 0; double DSIGJ = 0; double DSIGJP = 0;
73: double TEMP = 0;
74:
75: #endregion
76:
77: public DLALS0(DCOPY dcopy, DGEMV dgemv, DLACPY dlacpy, DLASCL dlascl, DROT drot, DSCAL dscal, XERBLA xerbla, DLAMC3 dlamc3, DNRM2 dnrm2)
78: {
79:
80:
81: #region Set Dependencies
82:
83: this._dcopy = dcopy; this._dgemv = dgemv; this._dlacpy = dlacpy; this._dlascl = dlascl; this._drot = drot;
84: this._dscal = dscal;this._xerbla = xerbla; this._dlamc3 = dlamc3; this._dnrm2 = dnrm2;
85:
86: #endregion
87:
88: }
89:
90: public DLALS0()
91: {
92:
93:
94: #region Dependencies (Initialization)
95:
96: DCOPY dcopy = new DCOPY();
97: LSAME lsame = new LSAME();
98: XERBLA xerbla = new XERBLA();
99: DLAMC3 dlamc3 = new DLAMC3();
100: DROT drot = new DROT();
101: DSCAL dscal = new DSCAL();
102: DNRM2 dnrm2 = new DNRM2();
103: DGEMV dgemv = new DGEMV(lsame, xerbla);
104: DLACPY dlacpy = new DLACPY(lsame);
105: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
106: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
107: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
108: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
109: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
110: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
111:
112: #endregion
113:
114:
115: #region Set Dependencies
116:
117: this._dcopy = dcopy; this._dgemv = dgemv; this._dlacpy = dlacpy; this._dlascl = dlascl; this._drot = drot;
118: this._dscal = dscal;this._xerbla = xerbla; this._dlamc3 = dlamc3; this._dnrm2 = dnrm2;
119:
120: #endregion
121:
122: }
123: /// <summary>
124: /// Purpose
125: /// =======
126: ///
127: /// DLALS0 applies back the multiplying factors of either the left or the
128: /// right singular vector matrix of a diagonal matrix appended by a row
129: /// to the right hand side matrix B in solving the least squares problem
130: /// using the divide-and-conquer SVD approach.
131: ///
132: /// For the left singular vector matrix, three types of orthogonal
133: /// matrices are involved:
134: ///
135: /// (1L) Givens rotations: the number of such rotations is GIVPTR; the
136: /// pairs of columns/rows they were applied to are stored in GIVCOL;
137: /// and the C- and S-values of these rotations are stored in GIVNUM.
138: ///
139: /// (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
140: /// row, and for J=2:N, PERM(J)-th row of B is to be moved to the
141: /// J-th row.
142: ///
143: /// (3L) The left singular vector matrix of the remaining matrix.
144: ///
145: /// For the right singular vector matrix, four types of orthogonal
146: /// matrices are involved:
147: ///
148: /// (1R) The right singular vector matrix of the remaining matrix.
149: ///
150: /// (2R) If SQRE = 1, one extra Givens rotation to generate the right
151: /// null space.
152: ///
153: /// (3R) The inverse transformation of (2L).
154: ///
155: /// (4R) The inverse transformation of (1L).
156: ///
157: ///</summary>
158: /// <param name="ICOMPQ">
159: /// (input) INTEGER
160: /// Specifies whether singular vectors are to be computed in
161: /// factored form:
162: /// = 0: Left singular vector matrix.
163: /// = 1: Right singular vector matrix.
164: ///</param>
165: /// <param name="NL">
166: /// (input) INTEGER
167: /// The row dimension of the upper block. NL .GE. 1.
168: ///</param>
169: /// <param name="NR">
170: /// (input) INTEGER
171: /// The row dimension of the lower block. NR .GE. 1.
172: ///</param>
173: /// <param name="SQRE">
174: /// (input) INTEGER
175: /// = 0: the lower block is an NR-by-NR square matrix.
176: /// = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
177: ///
178: /// The bidiagonal matrix has row dimension N = NL + NR + 1,
179: /// and column dimension M = N + SQRE.
180: ///</param>
181: /// <param name="NRHS">
182: /// (input) INTEGER
183: /// The number of columns of B and BX. NRHS must be at least 1.
184: ///</param>
185: /// <param name="B">
186: /// (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS )
187: /// On input, B contains the right hand sides of the least
188: /// squares problem in rows 1 through M. On output, B contains
189: /// the solution X in rows 1 through N.
190: ///</param>
191: /// <param name="LDB">
192: /// (input) INTEGER
193: /// The leading dimension of B. LDB must be at least
194: /// max(1,MAX( M, N ) ).
195: ///</param>
196: /// <param name="BX">
197: /// (workspace) DOUBLE PRECISION array, dimension ( LDBX, NRHS )
198: ///</param>
199: /// <param name="LDBX">
200: /// (input) INTEGER
201: /// The leading dimension of BX.
202: ///</param>
203: /// <param name="PERM">
204: /// (input) INTEGER array, dimension ( N )
205: /// The permutations (from deflation and sorting) applied
206: /// to the two blocks.
207: ///</param>
208: /// <param name="GIVPTR">
209: /// (input) INTEGER
210: /// The number of Givens rotations which took place in this
211: /// subproblem.
212: ///</param>
213: /// <param name="GIVCOL">
214: /// (input) INTEGER array, dimension ( LDGCOL, 2 )
215: /// Each pair of numbers indicates a pair of rows/columns
216: /// involved in a Givens rotation.
217: ///</param>
218: /// <param name="LDGCOL">
219: /// (input) INTEGER
220: /// The leading dimension of GIVCOL, must be at least N.
221: ///</param>
222: /// <param name="GIVNUM">
223: /// (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
224: /// Each number indicates the C or S value used in the
225: /// corresponding Givens rotation.
226: ///</param>
227: /// <param name="LDGNUM">
228: /// (input) INTEGER
229: /// The leading dimension of arrays DIFR, POLES and
230: /// GIVNUM, must be at least K.
231: ///</param>
232: /// <param name="POLES">
233: /// (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
234: /// On entry, POLES(1:K, 1) contains the new singular
235: /// values obtained from solving the secular equation, and
236: /// POLES(1:K, 2) is an array containing the poles in the secular
237: /// equation.
238: ///</param>
239: /// <param name="DIFL">
240: /// (input) DOUBLE PRECISION array, dimension ( K ).
241: /// On entry, DIFL(I) is the distance between I-th updated
242: /// (undeflated) singular value and the I-th (undeflated) old
243: /// singular value.
244: ///</param>
245: /// <param name="DIFR">
246: /// (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
247: /// On entry, DIFR(I, 1) contains the distances between I-th
248: /// updated (undeflated) singular value and the I+1-th
249: /// (undeflated) old singular value. And DIFR(I, 2) is the
250: /// normalizing factor for the I-th right singular vector.
251: ///</param>
252: /// <param name="Z">
253: /// (input) DOUBLE PRECISION array, dimension ( K )
254: /// Contain the components of the deflation-adjusted updating row
255: /// vector.
256: ///</param>
257: /// <param name="K">
258: /// (input) INTEGER
259: /// Contains the dimension of the non-deflated matrix,
260: /// This is the order of the related secular equation. 1 .LE. K .LE.N.
261: ///</param>
262: /// <param name="C">
263: /// (input) DOUBLE PRECISION
264: /// C contains garbage if SQRE =0 and the C-value of a Givens
265: /// rotation related to the right null space if SQRE = 1.
266: ///</param>
267: /// <param name="S">
268: /// (input) DOUBLE PRECISION
269: /// S contains garbage if SQRE =0 and the S-value of a Givens
270: /// rotation related to the right null space if SQRE = 1.
271: ///</param>
272: /// <param name="WORK">
273: /// (workspace) DOUBLE PRECISION array, dimension ( K )
274: ///</param>
275: /// <param name="INFO">
276: /// (output) INTEGER
277: /// = 0: successful exit.
278: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
279: ///</param>
280: public void Run(int ICOMPQ, int NL, int NR, int SQRE, int NRHS, ref double[] B, int offset_b
281: , int LDB, ref double[] BX, int offset_bx, int LDBX, int[] PERM, int offset_perm, int GIVPTR, int[] GIVCOL, int offset_givcol
282: , int LDGCOL, double[] GIVNUM, int offset_givnum, int LDGNUM, double[] POLES, int offset_poles, double[] DIFL, int offset_difl, double[] DIFR, int offset_difr
283: , double[] Z, int offset_z, int K, double C, double S, ref double[] WORK, int offset_work, ref int INFO)
284: {
285:
286: #region Array Index Correction
287:
288: int o_b = -1 - LDB + offset_b; int o_bx = -1 - LDBX + offset_bx; int o_perm = -1 + offset_perm;
289: int o_givcol = -1 - LDGCOL + offset_givcol; int o_givnum = -1 - LDGNUM + offset_givnum;
290: int o_poles = -1 - LDGNUM + offset_poles; int o_difl = -1 + offset_difl; int o_difr = -1 - LDGNUM + offset_difr;
291: int o_z = -1 + offset_z; int o_work = -1 + offset_work;
292:
293: #endregion
294:
295:
296: #region Prolog
297:
298: // *
299: // * -- LAPACK routine (version 3.1) --
300: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
301: // * November 2006
302: // *
303: // * .. Scalar Arguments ..
304: // * ..
305: // * .. Array Arguments ..
306: // * ..
307: // *
308: // * Purpose
309: // * =======
310: // *
311: // * DLALS0 applies back the multiplying factors of either the left or the
312: // * right singular vector matrix of a diagonal matrix appended by a row
313: // * to the right hand side matrix B in solving the least squares problem
314: // * using the divide-and-conquer SVD approach.
315: // *
316: // * For the left singular vector matrix, three types of orthogonal
317: // * matrices are involved:
318: // *
319: // * (1L) Givens rotations: the number of such rotations is GIVPTR; the
320: // * pairs of columns/rows they were applied to are stored in GIVCOL;
321: // * and the C- and S-values of these rotations are stored in GIVNUM.
322: // *
323: // * (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
324: // * row, and for J=2:N, PERM(J)-th row of B is to be moved to the
325: // * J-th row.
326: // *
327: // * (3L) The left singular vector matrix of the remaining matrix.
328: // *
329: // * For the right singular vector matrix, four types of orthogonal
330: // * matrices are involved:
331: // *
332: // * (1R) The right singular vector matrix of the remaining matrix.
333: // *
334: // * (2R) If SQRE = 1, one extra Givens rotation to generate the right
335: // * null space.
336: // *
337: // * (3R) The inverse transformation of (2L).
338: // *
339: // * (4R) The inverse transformation of (1L).
340: // *
341: // * Arguments
342: // * =========
343: // *
344: // * ICOMPQ (input) INTEGER
345: // * Specifies whether singular vectors are to be computed in
346: // * factored form:
347: // * = 0: Left singular vector matrix.
348: // * = 1: Right singular vector matrix.
349: // *
350: // * NL (input) INTEGER
351: // * The row dimension of the upper block. NL >= 1.
352: // *
353: // * NR (input) INTEGER
354: // * The row dimension of the lower block. NR >= 1.
355: // *
356: // * SQRE (input) INTEGER
357: // * = 0: the lower block is an NR-by-NR square matrix.
358: // * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
359: // *
360: // * The bidiagonal matrix has row dimension N = NL + NR + 1,
361: // * and column dimension M = N + SQRE.
362: // *
363: // * NRHS (input) INTEGER
364: // * The number of columns of B and BX. NRHS must be at least 1.
365: // *
366: // * B (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS )
367: // * On input, B contains the right hand sides of the least
368: // * squares problem in rows 1 through M. On output, B contains
369: // * the solution X in rows 1 through N.
370: // *
371: // * LDB (input) INTEGER
372: // * The leading dimension of B. LDB must be at least
373: // * max(1,MAX( M, N ) ).
374: // *
375: // * BX (workspace) DOUBLE PRECISION array, dimension ( LDBX, NRHS )
376: // *
377: // * LDBX (input) INTEGER
378: // * The leading dimension of BX.
379: // *
380: // * PERM (input) INTEGER array, dimension ( N )
381: // * The permutations (from deflation and sorting) applied
382: // * to the two blocks.
383: // *
384: // * GIVPTR (input) INTEGER
385: // * The number of Givens rotations which took place in this
386: // * subproblem.
387: // *
388: // * GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 )
389: // * Each pair of numbers indicates a pair of rows/columns
390: // * involved in a Givens rotation.
391: // *
392: // * LDGCOL (input) INTEGER
393: // * The leading dimension of GIVCOL, must be at least N.
394: // *
395: // * GIVNUM (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
396: // * Each number indicates the C or S value used in the
397: // * corresponding Givens rotation.
398: // *
399: // * LDGNUM (input) INTEGER
400: // * The leading dimension of arrays DIFR, POLES and
401: // * GIVNUM, must be at least K.
402: // *
403: // * POLES (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
404: // * On entry, POLES(1:K, 1) contains the new singular
405: // * values obtained from solving the secular equation, and
406: // * POLES(1:K, 2) is an array containing the poles in the secular
407: // * equation.
408: // *
409: // * DIFL (input) DOUBLE PRECISION array, dimension ( K ).
410: // * On entry, DIFL(I) is the distance between I-th updated
411: // * (undeflated) singular value and the I-th (undeflated) old
412: // * singular value.
413: // *
414: // * DIFR (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
415: // * On entry, DIFR(I, 1) contains the distances between I-th
416: // * updated (undeflated) singular value and the I+1-th
417: // * (undeflated) old singular value. And DIFR(I, 2) is the
418: // * normalizing factor for the I-th right singular vector.
419: // *
420: // * Z (input) DOUBLE PRECISION array, dimension ( K )
421: // * Contain the components of the deflation-adjusted updating row
422: // * vector.
423: // *
424: // * K (input) INTEGER
425: // * Contains the dimension of the non-deflated matrix,
426: // * This is the order of the related secular equation. 1 <= K <=N.
427: // *
428: // * C (input) DOUBLE PRECISION
429: // * C contains garbage if SQRE =0 and the C-value of a Givens
430: // * rotation related to the right null space if SQRE = 1.
431: // *
432: // * S (input) DOUBLE PRECISION
433: // * S contains garbage if SQRE =0 and the S-value of a Givens
434: // * rotation related to the right null space if SQRE = 1.
435: // *
436: // * WORK (workspace) DOUBLE PRECISION array, dimension ( K )
437: // *
438: // * INFO (output) INTEGER
439: // * = 0: successful exit.
440: // * < 0: if INFO = -i, the i-th argument had an illegal value.
441: // *
442: // * Further Details
443: // * ===============
444: // *
445: // * Based on contributions by
446: // * Ming Gu and Ren-Cang Li, Computer Science Division, University of
447: // * California at Berkeley, USA
448: // * Osni Marques, LBNL/NERSC, USA
449: // *
450: // * =====================================================================
451: // *
452: // * .. Parameters ..
453: // * ..
454: // * .. Local Scalars ..
455: // * ..
456: // * .. External Subroutines ..
457: // * ..
458: // * .. External Functions ..
459: // * ..
460: // * .. Intrinsic Functions ..
461: // INTRINSIC MAX;
462: // * ..
463: // * .. Executable Statements ..
464: // *
465: // * Test the input parameters.
466: // *
467:
468: #endregion
469:
470:
471: #region Body
472:
473: INFO = 0;
474: // *
475: if ((ICOMPQ < 0) || (ICOMPQ > 1))
476: {
477: INFO = - 1;
478: }
479: else
480: {
481: if (NL < 1)
482: {
483: INFO = - 2;
484: }
485: else
486: {
487: if (NR < 1)
488: {
489: INFO = - 3;
490: }
491: else
492: {
493: if ((SQRE < 0) || (SQRE > 1))
494: {
495: INFO = - 4;
496: }
497: }
498: }
499: }
500: // *
501: N = NL + NR + 1;
502: // *
503: if (NRHS < 1)
504: {
505: INFO = - 5;
506: }
507: else
508: {
509: if (LDB < N)
510: {
511: INFO = - 7;
512: }
513: else
514: {
515: if (LDBX < N)
516: {
517: INFO = - 9;
518: }
519: else
520: {
521: if (GIVPTR < 0)
522: {
523: INFO = - 11;
524: }
525: else
526: {
527: if (LDGCOL < N)
528: {
529: INFO = - 13;
530: }
531: else
532: {
533: if (LDGNUM < N)
534: {
535: INFO = - 15;
536: }
537: else
538: {
539: if (K < 1)
540: {
541: INFO = - 20;
542: }
543: }
544: }
545: }
546: }
547: }
548: }
549: if (INFO != 0)
550: {
551: this._xerbla.Run("DLALS0", - INFO);
552: return;
553: }
554: // *
555: M = N + SQRE;
556: NLP1 = NL + 1;
557: // *
558: if (ICOMPQ == 0)
559: {
560: // *
561: // * Apply back orthogonal transformations from the left.
562: // *
563: // * Step (1L): apply back the Givens rotations performed.
564: // *
565: for (I = 1; I <= GIVPTR; I++)
566: {
567: this._drot.Run(NRHS, ref B, GIVCOL[I+2 * LDGCOL + o_givcol]+1 * LDB + o_b, LDB, ref B, GIVCOL[I+1 * LDGCOL + o_givcol]+1 * LDB + o_b, LDB, GIVNUM[I+2 * LDGNUM + o_givnum]
568: , GIVNUM[I+1 * LDGNUM + o_givnum]);
569: }
570: // *
571: // * Step (2L): permute rows of B.
572: // *
573: this._dcopy.Run(NRHS, B, NLP1+1 * LDB + o_b, LDB, ref BX, 1+1 * LDBX + o_bx, LDBX);
574: for (I = 2; I <= N; I++)
575: {
576: this._dcopy.Run(NRHS, B, PERM[I + o_perm]+1 * LDB + o_b, LDB, ref BX, I+1 * LDBX + o_bx, LDBX);
577: }
578: // *
579: // * Step (3L): apply the inverse of the left singular vector
580: // * matrix to BX.
581: // *
582: if (K == 1)
583: {
584: this._dcopy.Run(NRHS, BX, offset_bx, LDBX, ref B, offset_b, LDB);
585: if (Z[1 + o_z] < ZERO)
586: {
587: this._dscal.Run(NRHS, NEGONE, ref B, offset_b, LDB);
588: }
589: }
590: else
591: {
592: for (J = 1; J <= K; J++)
593: {
594: DIFLJ = DIFL[J + o_difl];
595: DJ = POLES[J+1 * LDGNUM + o_poles];
596: DSIGJ = - POLES[J+2 * LDGNUM + o_poles];
597: if (J < K)
598: {
599: DIFRJ = - DIFR[J+1 * LDGNUM + o_difr];
600: DSIGJP = - POLES[J + 1+2 * LDGNUM + o_poles];
601: }
602: if ((Z[J + o_z] == ZERO) || (POLES[J+2 * LDGNUM + o_poles] == ZERO))
603: {
604: WORK[J + o_work] = ZERO;
605: }
606: else
607: {
608: WORK[J + o_work] = - POLES[J+2 * LDGNUM + o_poles] * Z[J + o_z] / DIFLJ / (POLES[J+2 * LDGNUM + o_poles] + DJ);
609: }
610: for (I = 1; I <= J - 1; I++)
611: {
612: if ((Z[I + o_z] == ZERO) || (POLES[I+2 * LDGNUM + o_poles] == ZERO))
613: {
614: WORK[I + o_work] = ZERO;
615: }
616: else
617: {
618: WORK[I + o_work] = POLES[I+2 * LDGNUM + o_poles] * Z[I + o_z] / (this._dlamc3.Run(POLES[I+2 * LDGNUM + o_poles], DSIGJ) - DIFLJ) / (POLES[I+2 * LDGNUM + o_poles] + DJ);
619: }
620: }
621: for (I = J + 1; I <= K; I++)
622: {
623: if ((Z[I + o_z] == ZERO) || (POLES[I+2 * LDGNUM + o_poles] == ZERO))
624: {
625: WORK[I + o_work] = ZERO;
626: }
627: else
628: {
629: WORK[I + o_work] = POLES[I+2 * LDGNUM + o_poles] * Z[I + o_z] / (this._dlamc3.Run(POLES[I+2 * LDGNUM + o_poles], DSIGJP) + DIFRJ) / (POLES[I+2 * LDGNUM + o_poles] + DJ);
630: }
631: }
632: WORK[1 + o_work] = NEGONE;
633: TEMP = this._dnrm2.Run(K, WORK, offset_work, 1);
634: this._dgemv.Run("T", K, NRHS, ONE, BX, offset_bx, LDBX
635: , WORK, offset_work, 1, ZERO, ref B, J+1 * LDB + o_b, LDB);
636: this._dlascl.Run("G", 0, 0, TEMP, ONE, 1
637: , NRHS, ref B, J+1 * LDB + o_b, LDB, ref INFO);
638: }
639: }
640: // *
641: // * Move the deflated rows of BX to B also.
642: // *
643: if (K < Math.Max(M, N))
644: {
645: this._dlacpy.Run("A", N - K, NRHS, BX, K + 1+1 * LDBX + o_bx, LDBX, ref B, K + 1+1 * LDB + o_b
646: , LDB);
647: }
648: }
649: else
650: {
651: // *
652: // * Apply back the right orthogonal transformations.
653: // *
654: // * Step (1R): apply back the new right singular vector matrix
655: // * to B.
656: // *
657: if (K == 1)
658: {
659: this._dcopy.Run(NRHS, B, offset_b, LDB, ref BX, offset_bx, LDBX);
660: }
661: else
662: {
663: for (J = 1; J <= K; J++)
664: {
665: DSIGJ = POLES[J+2 * LDGNUM + o_poles];
666: if (Z[J + o_z] == ZERO)
667: {
668: WORK[J + o_work] = ZERO;
669: }
670: else
671: {
672: WORK[J + o_work] = - Z[J + o_z] / DIFL[J + o_difl] / (DSIGJ + POLES[J+1 * LDGNUM + o_poles]) / DIFR[J+2 * LDGNUM + o_difr];
673: }
674: for (I = 1; I <= J - 1; I++)
675: {
676: if (Z[J + o_z] == ZERO)
677: {
678: WORK[I + o_work] = ZERO;
679: }
680: else
681: {
682: WORK[I + o_work] = Z[J + o_z] / (this._dlamc3.Run(DSIGJ, - POLES[I + 1+2 * LDGNUM + o_poles]) - DIFR[I+1 * LDGNUM + o_difr]) / (DSIGJ + POLES[I+1 * LDGNUM + o_poles]) / DIFR[I+2 * LDGNUM + o_difr];
683: }
684: }
685: for (I = J + 1; I <= K; I++)
686: {
687: if (Z[J + o_z] == ZERO)
688: {
689: WORK[I + o_work] = ZERO;
690: }
691: else
692: {
693: WORK[I + o_work] = Z[J + o_z] / (this._dlamc3.Run(DSIGJ, - POLES[I+2 * LDGNUM + o_poles]) - DIFL[I + o_difl]) / (DSIGJ + POLES[I+1 * LDGNUM + o_poles]) / DIFR[I+2 * LDGNUM + o_difr];
694: }
695: }
696: this._dgemv.Run("T", K, NRHS, ONE, B, offset_b, LDB
697: , WORK, offset_work, 1, ZERO, ref BX, J+1 * LDBX + o_bx, LDBX);
698: }
699: }
700: // *
701: // * Step (2R): if SQRE = 1, apply back the rotation that is
702: // * related to the right null space of the subproblem.
703: // *
704: if (SQRE == 1)
705: {
706: this._dcopy.Run(NRHS, B, M+1 * LDB + o_b, LDB, ref BX, M+1 * LDBX + o_bx, LDBX);
707: this._drot.Run(NRHS, ref BX, 1+1 * LDBX + o_bx, LDBX, ref BX, M+1 * LDBX + o_bx, LDBX, C
708: , S);
709: }
710: if (K < Math.Max(M, N))
711: {
712: this._dlacpy.Run("A", N - K, NRHS, B, K + 1+1 * LDB + o_b, LDB, ref BX, K + 1+1 * LDBX + o_bx
713: , LDBX);
714: }
715: // *
716: // * Step (3R): permute rows of B.
717: // *
718: this._dcopy.Run(NRHS, BX, 1+1 * LDBX + o_bx, LDBX, ref B, NLP1+1 * LDB + o_b, LDB);
719: if (SQRE == 1)
720: {
721: this._dcopy.Run(NRHS, BX, M+1 * LDBX + o_bx, LDBX, ref B, M+1 * LDB + o_b, LDB);
722: }
723: for (I = 2; I <= N; I++)
724: {
725: this._dcopy.Run(NRHS, BX, I+1 * LDBX + o_bx, LDBX, ref B, PERM[I + o_perm]+1 * LDB + o_b, LDB);
726: }
727: // *
728: // * Step (4R): apply back the Givens rotations performed.
729: // *
730: for (I = GIVPTR; I >= 1; I += - 1)
731: {
732: this._drot.Run(NRHS, ref B, GIVCOL[I+2 * LDGCOL + o_givcol]+1 * LDB + o_b, LDB, ref B, GIVCOL[I+1 * LDGCOL + o_givcol]+1 * LDB + o_b, LDB, GIVNUM[I+2 * LDGNUM + o_givnum]
733: , - GIVNUM[I+1 * LDGNUM + o_givnum]);
734: }
735: }
736: // *
737: return;
738: // *
739: // * End of DLALS0
740: // *
741:
742: #endregion
743:
744: }
745: }
746: }