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: /// DLALSA is an itermediate step in solving the least squares problem
27: /// by computing the SVD of the coefficient matrix in compact form (The
28: /// singular vectors are computed as products of simple orthorgonal
29: /// matrices.).
30: ///
31: /// If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector
32: /// matrix of an upper bidiagonal matrix to the right hand side; and if
33: /// ICOMPQ = 1, DLALSA applies the right singular vector matrix to the
34: /// right hand side. The singular vector matrices were generated in
35: /// compact form by DLALSA.
36: ///
37: ///</summary>
38: public class DLALSA
39: {
40:
41:
42: #region Dependencies
43:
44: DCOPY _dcopy; DGEMM _dgemm; DLALS0 _dlals0; DLASDT _dlasdt; XERBLA _xerbla;
45:
46: #endregion
47:
48:
49: #region Fields
50:
51: const double ZERO = 0.0E0; const double ONE = 1.0E0; int I = 0; int I1 = 0; int IC = 0; int IM1 = 0; int INODE = 0;
52: int J = 0;int LF = 0; int LL = 0; int LVL = 0; int LVL2 = 0; int ND = 0; int NDB1 = 0; int NDIML = 0; int NDIMR = 0;
53: int NL = 0;int NLF = 0; int NLP1 = 0; int NLVL = 0; int NR = 0; int NRF = 0; int NRP1 = 0; int SQRE = 0;
54:
55: #endregion
56:
57: public DLALSA(DCOPY dcopy, DGEMM dgemm, DLALS0 dlals0, DLASDT dlasdt, XERBLA xerbla)
58: {
59:
60:
61: #region Set Dependencies
62:
63: this._dcopy = dcopy; this._dgemm = dgemm; this._dlals0 = dlals0; this._dlasdt = dlasdt; this._xerbla = xerbla;
64:
65: #endregion
66:
67: }
68:
69: public DLALSA()
70: {
71:
72:
73: #region Dependencies (Initialization)
74:
75: DCOPY dcopy = new DCOPY();
76: LSAME lsame = new LSAME();
77: XERBLA xerbla = new XERBLA();
78: DLAMC3 dlamc3 = new DLAMC3();
79: DROT drot = new DROT();
80: DSCAL dscal = new DSCAL();
81: DNRM2 dnrm2 = new DNRM2();
82: DLASDT dlasdt = new DLASDT();
83: DGEMM dgemm = new DGEMM(lsame, xerbla);
84: DGEMV dgemv = new DGEMV(lsame, xerbla);
85: DLACPY dlacpy = new DLACPY(lsame);
86: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
87: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
88: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
89: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
90: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
91: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
92: DLALS0 dlals0 = new DLALS0(dcopy, dgemv, dlacpy, dlascl, drot, dscal, xerbla, dlamc3, dnrm2);
93:
94: #endregion
95:
96:
97: #region Set Dependencies
98:
99: this._dcopy = dcopy; this._dgemm = dgemm; this._dlals0 = dlals0; this._dlasdt = dlasdt; this._xerbla = xerbla;
100:
101: #endregion
102:
103: }
104: /// <summary>
105: /// Purpose
106: /// =======
107: ///
108: /// DLALSA is an itermediate step in solving the least squares problem
109: /// by computing the SVD of the coefficient matrix in compact form (The
110: /// singular vectors are computed as products of simple orthorgonal
111: /// matrices.).
112: ///
113: /// If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector
114: /// matrix of an upper bidiagonal matrix to the right hand side; and if
115: /// ICOMPQ = 1, DLALSA applies the right singular vector matrix to the
116: /// right hand side. The singular vector matrices were generated in
117: /// compact form by DLALSA.
118: ///
119: ///</summary>
120: /// <param name="ICOMPQ">
121: /// (input) INTEGER
122: /// Specifies whether the left or the right singular vector
123: /// matrix is involved.
124: /// = 0: Left singular vector matrix
125: /// = 1: Right singular vector matrix
126: ///</param>
127: /// <param name="SMLSIZ">
128: /// (input) INTEGER
129: /// The maximum size of the subproblems at the bottom of the
130: /// computation tree.
131: ///</param>
132: /// <param name="N">
133: /// (input) INTEGER
134: /// The row and column dimensions of the upper bidiagonal matrix.
135: ///</param>
136: /// <param name="NRHS">
137: /// (input) INTEGER
138: /// The number of columns of B and BX. NRHS must be at least 1.
139: ///</param>
140: /// <param name="B">
141: /// (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS )
142: /// On input, B contains the right hand sides of the least
143: /// squares problem in rows 1 through M.
144: /// On output, B contains the solution X in rows 1 through N.
145: ///</param>
146: /// <param name="LDB">
147: /// (input) INTEGER
148: /// The leading dimension of B in the calling subprogram.
149: /// LDB must be at least max(1,MAX( M, N ) ).
150: ///</param>
151: /// <param name="BX">
152: /// (output) DOUBLE PRECISION array, dimension ( LDBX, NRHS )
153: /// On exit, the result of applying the left or right singular
154: /// vector matrix to B.
155: ///</param>
156: /// <param name="LDBX">
157: /// (input) INTEGER
158: /// The leading dimension of BX.
159: ///</param>
160: /// <param name="U">
161: /// (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
162: /// On entry, U contains the left singular vector matrices of all
163: /// subproblems at the bottom level.
164: ///</param>
165: /// <param name="LDU">
166: /// (input) INTEGER, LDU = .GT. N.
167: /// The leading dimension of arrays U, VT, DIFL, DIFR,
168: /// POLES, GIVNUM, and Z.
169: ///</param>
170: /// <param name="VT">
171: /// (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).
172: /// On entry, VT' contains the right singular vector matrices of
173: /// all subproblems at the bottom level.
174: ///</param>
175: /// <param name="K">
176: /// (input) INTEGER array, dimension ( N ).
177: ///</param>
178: /// <param name="DIFL">
179: /// (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
180: /// where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
181: ///</param>
182: /// <param name="DIFR">
183: /// (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
184: /// On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
185: /// distances between singular values on the I-th level and
186: /// singular values on the (I -1)-th level, and DIFR(*, 2 * I)
187: /// record the normalizing factors of the right singular vectors
188: /// matrices of subproblems on I-th level.
189: ///</param>
190: /// <param name="Z">
191: /// (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
192: /// On entry, Z(1, I) contains the components of the deflation-
193: /// adjusted updating row vector for subproblems on the I-th
194: /// level.
195: ///</param>
196: /// <param name="POLES">
197: /// (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
198: /// On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
199: /// singular values involved in the secular equations on the I-th
200: /// level.
201: ///</param>
202: /// <param name="GIVPTR">
203: /// (input) INTEGER array, dimension ( N ).
204: /// On entry, GIVPTR( I ) records the number of Givens
205: /// rotations performed on the I-th problem on the computation
206: /// tree.
207: ///</param>
208: /// <param name="GIVCOL">
209: /// (input) INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
210: /// On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
211: /// locations of Givens rotations performed on the I-th level on
212: /// the computation tree.
213: ///</param>
214: /// <param name="LDGCOL">
215: /// (input) INTEGER, LDGCOL = .GT. N.
216: /// The leading dimension of arrays GIVCOL and PERM.
217: ///</param>
218: /// <param name="PERM">
219: /// (input) INTEGER array, dimension ( LDGCOL, NLVL ).
220: /// On entry, PERM(*, I) records permutations done on the I-th
221: /// level of the computation tree.
222: ///</param>
223: /// <param name="GIVNUM">
224: /// (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
225: /// On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
226: /// values of Givens rotations performed on the I-th level on the
227: /// computation tree.
228: ///</param>
229: /// <param name="C">
230: /// (input) DOUBLE PRECISION array, dimension ( N ).
231: /// On entry, if the I-th subproblem is not square,
232: /// C( I ) contains the C-value of a Givens rotation related to
233: /// the right null space of the I-th subproblem.
234: ///</param>
235: /// <param name="S">
236: /// (input) DOUBLE PRECISION array, dimension ( N ).
237: /// On entry, if the I-th subproblem is not square,
238: /// S( I ) contains the S-value of a Givens rotation related to
239: /// the right null space of the I-th subproblem.
240: ///</param>
241: /// <param name="WORK">
242: /// (workspace) DOUBLE PRECISION array.
243: /// The dimension must be at least N.
244: ///</param>
245: /// <param name="IWORK">
246: /// (workspace) INTEGER array.
247: /// The dimension must be at least 3 * N
248: ///</param>
249: /// <param name="INFO">
250: /// (output) INTEGER
251: /// = 0: successful exit.
252: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
253: ///</param>
254: public void Run(int ICOMPQ, int SMLSIZ, int N, int NRHS, ref double[] B, int offset_b, int LDB
255: , ref double[] BX, int offset_bx, int LDBX, double[] U, int offset_u, int LDU, double[] VT, int offset_vt, int[] K, int offset_k
256: , double[] DIFL, int offset_difl, double[] DIFR, int offset_difr, double[] Z, int offset_z, double[] POLES, int offset_poles, int[] GIVPTR, int offset_givptr, int[] GIVCOL, int offset_givcol
257: , int LDGCOL, int[] PERM, int offset_perm, double[] GIVNUM, int offset_givnum, double[] C, int offset_c, double[] S, int offset_s, ref double[] WORK, int offset_work
258: , ref int[] IWORK, int offset_iwork, ref int INFO)
259: {
260:
261: #region Array Index Correction
262:
263: int o_b = -1 - LDB + offset_b; int o_bx = -1 - LDBX + offset_bx; int o_u = -1 - LDU + offset_u;
264: int o_vt = -1 - LDU + offset_vt; int o_k = -1 + offset_k; int o_difl = -1 - LDU + offset_difl;
265: int o_difr = -1 - LDU + offset_difr; int o_z = -1 - LDU + offset_z; int o_poles = -1 - LDU + offset_poles;
266: int o_givptr = -1 + offset_givptr; int o_givcol = -1 - LDGCOL + offset_givcol;
267: int o_perm = -1 - LDGCOL + offset_perm; int o_givnum = -1 - LDU + offset_givnum; int o_c = -1 + offset_c;
268: int o_s = -1 + offset_s; int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
269:
270: #endregion
271:
272:
273: #region Prolog
274:
275: // *
276: // * -- LAPACK routine (version 3.1) --
277: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
278: // * November 2006
279: // *
280: // * .. Scalar Arguments ..
281: // * ..
282: // * .. Array Arguments ..
283: // * ..
284: // *
285: // * Purpose
286: // * =======
287: // *
288: // * DLALSA is an itermediate step in solving the least squares problem
289: // * by computing the SVD of the coefficient matrix in compact form (The
290: // * singular vectors are computed as products of simple orthorgonal
291: // * matrices.).
292: // *
293: // * If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector
294: // * matrix of an upper bidiagonal matrix to the right hand side; and if
295: // * ICOMPQ = 1, DLALSA applies the right singular vector matrix to the
296: // * right hand side. The singular vector matrices were generated in
297: // * compact form by DLALSA.
298: // *
299: // * Arguments
300: // * =========
301: // *
302: // *
303: // * ICOMPQ (input) INTEGER
304: // * Specifies whether the left or the right singular vector
305: // * matrix is involved.
306: // * = 0: Left singular vector matrix
307: // * = 1: Right singular vector matrix
308: // *
309: // * SMLSIZ (input) INTEGER
310: // * The maximum size of the subproblems at the bottom of the
311: // * computation tree.
312: // *
313: // * N (input) INTEGER
314: // * The row and column dimensions of the upper bidiagonal matrix.
315: // *
316: // * NRHS (input) INTEGER
317: // * The number of columns of B and BX. NRHS must be at least 1.
318: // *
319: // * B (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS )
320: // * On input, B contains the right hand sides of the least
321: // * squares problem in rows 1 through M.
322: // * On output, B contains the solution X in rows 1 through N.
323: // *
324: // * LDB (input) INTEGER
325: // * The leading dimension of B in the calling subprogram.
326: // * LDB must be at least max(1,MAX( M, N ) ).
327: // *
328: // * BX (output) DOUBLE PRECISION array, dimension ( LDBX, NRHS )
329: // * On exit, the result of applying the left or right singular
330: // * vector matrix to B.
331: // *
332: // * LDBX (input) INTEGER
333: // * The leading dimension of BX.
334: // *
335: // * U (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
336: // * On entry, U contains the left singular vector matrices of all
337: // * subproblems at the bottom level.
338: // *
339: // * LDU (input) INTEGER, LDU = > N.
340: // * The leading dimension of arrays U, VT, DIFL, DIFR,
341: // * POLES, GIVNUM, and Z.
342: // *
343: // * VT (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).
344: // * On entry, VT' contains the right singular vector matrices of
345: // * all subproblems at the bottom level.
346: // *
347: // * K (input) INTEGER array, dimension ( N ).
348: // *
349: // * DIFL (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
350: // * where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
351: // *
352: // * DIFR (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
353: // * On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
354: // * distances between singular values on the I-th level and
355: // * singular values on the (I -1)-th level, and DIFR(*, 2 * I)
356: // * record the normalizing factors of the right singular vectors
357: // * matrices of subproblems on I-th level.
358: // *
359: // * Z (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
360: // * On entry, Z(1, I) contains the components of the deflation-
361: // * adjusted updating row vector for subproblems on the I-th
362: // * level.
363: // *
364: // * POLES (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
365: // * On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
366: // * singular values involved in the secular equations on the I-th
367: // * level.
368: // *
369: // * GIVPTR (input) INTEGER array, dimension ( N ).
370: // * On entry, GIVPTR( I ) records the number of Givens
371: // * rotations performed on the I-th problem on the computation
372: // * tree.
373: // *
374: // * GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
375: // * On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
376: // * locations of Givens rotations performed on the I-th level on
377: // * the computation tree.
378: // *
379: // * LDGCOL (input) INTEGER, LDGCOL = > N.
380: // * The leading dimension of arrays GIVCOL and PERM.
381: // *
382: // * PERM (input) INTEGER array, dimension ( LDGCOL, NLVL ).
383: // * On entry, PERM(*, I) records permutations done on the I-th
384: // * level of the computation tree.
385: // *
386: // * GIVNUM (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
387: // * On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
388: // * values of Givens rotations performed on the I-th level on the
389: // * computation tree.
390: // *
391: // * C (input) DOUBLE PRECISION array, dimension ( N ).
392: // * On entry, if the I-th subproblem is not square,
393: // * C( I ) contains the C-value of a Givens rotation related to
394: // * the right null space of the I-th subproblem.
395: // *
396: // * S (input) DOUBLE PRECISION array, dimension ( N ).
397: // * On entry, if the I-th subproblem is not square,
398: // * S( I ) contains the S-value of a Givens rotation related to
399: // * the right null space of the I-th subproblem.
400: // *
401: // * WORK (workspace) DOUBLE PRECISION array.
402: // * The dimension must be at least N.
403: // *
404: // * IWORK (workspace) INTEGER array.
405: // * The dimension must be at least 3 * N
406: // *
407: // * INFO (output) INTEGER
408: // * = 0: successful exit.
409: // * < 0: if INFO = -i, the i-th argument had an illegal value.
410: // *
411: // * Further Details
412: // * ===============
413: // *
414: // * Based on contributions by
415: // * Ming Gu and Ren-Cang Li, Computer Science Division, University of
416: // * California at Berkeley, USA
417: // * Osni Marques, LBNL/NERSC, USA
418: // *
419: // * =====================================================================
420: // *
421: // * .. Parameters ..
422: // * ..
423: // * .. Local Scalars ..
424: // * ..
425: // * .. External Subroutines ..
426: // * ..
427: // * .. Executable Statements ..
428: // *
429: // * Test the input parameters.
430: // *
431:
432: #endregion
433:
434:
435: #region Body
436:
437: INFO = 0;
438: // *
439: if ((ICOMPQ < 0) || (ICOMPQ > 1))
440: {
441: INFO = - 1;
442: }
443: else
444: {
445: if (SMLSIZ < 3)
446: {
447: INFO = - 2;
448: }
449: else
450: {
451: if (N < SMLSIZ)
452: {
453: INFO = - 3;
454: }
455: else
456: {
457: if (NRHS < 1)
458: {
459: INFO = - 4;
460: }
461: else
462: {
463: if (LDB < N)
464: {
465: INFO = - 6;
466: }
467: else
468: {
469: if (LDBX < N)
470: {
471: INFO = - 8;
472: }
473: else
474: {
475: if (LDU < N)
476: {
477: INFO = - 10;
478: }
479: else
480: {
481: if (LDGCOL < N)
482: {
483: INFO = - 19;
484: }
485: }
486: }
487: }
488: }
489: }
490: }
491: }
492: if (INFO != 0)
493: {
494: this._xerbla.Run("DLALSA", - INFO);
495: return;
496: }
497: // *
498: // * Book-keeping and setting up the computation tree.
499: // *
500: INODE = 1;
501: NDIML = INODE + N;
502: NDIMR = NDIML + N;
503: // *
504: this._dlasdt.Run(N, ref NLVL, ref ND, ref IWORK, INODE + o_iwork, ref IWORK, NDIML + o_iwork, ref IWORK, NDIMR + o_iwork
505: , SMLSIZ);
506: // *
507: // * The following code applies back the left singular vector factors.
508: // * For applying back the right singular vector factors, go to 50.
509: // *
510: if (ICOMPQ == 1)
511: {
512: goto LABEL50;
513: }
514: // *
515: // * The nodes on the bottom level of the tree were solved
516: // * by DLASDQ. The corresponding left and right singular vector
517: // * matrices are in explicit form. First apply back the left
518: // * singular vector matrices.
519: // *
520: NDB1 = (ND + 1) / 2;
521: for (I = NDB1; I <= ND; I++)
522: {
523: // *
524: // * IC : center row of each node
525: // * NL : number of rows of left subproblem
526: // * NR : number of rows of right subproblem
527: // * NLF: starting row of the left subproblem
528: // * NRF: starting row of the right subproblem
529: // *
530: I1 = I - 1;
531: IC = IWORK[INODE + I1 + o_iwork];
532: NL = IWORK[NDIML + I1 + o_iwork];
533: NR = IWORK[NDIMR + I1 + o_iwork];
534: NLF = IC - NL;
535: NRF = IC + 1;
536: this._dgemm.Run("T", "N", NL, NRHS, NL, ONE
537: , U, NLF+1 * LDU + o_u, LDU, B, NLF+1 * LDB + o_b, LDB, ZERO, ref BX, NLF+1 * LDBX + o_bx
538: , LDBX);
539: this._dgemm.Run("T", "N", NR, NRHS, NR, ONE
540: , U, NRF+1 * LDU + o_u, LDU, B, NRF+1 * LDB + o_b, LDB, ZERO, ref BX, NRF+1 * LDBX + o_bx
541: , LDBX);
542: }
543: // *
544: // * Next copy the rows of B that correspond to unchanged rows
545: // * in the bidiagonal matrix to BX.
546: // *
547: for (I = 1; I <= ND; I++)
548: {
549: IC = IWORK[INODE + I - 1 + o_iwork];
550: this._dcopy.Run(NRHS, B, IC+1 * LDB + o_b, LDB, ref BX, IC+1 * LDBX + o_bx, LDBX);
551: }
552: // *
553: // * Finally go through the left singular vector matrices of all
554: // * the other subproblems bottom-up on the tree.
555: // *
556: J = (int)Math.Pow(2, NLVL);
557: SQRE = 0;
558: // *
559: for (LVL = NLVL; LVL >= 1; LVL += - 1)
560: {
561: LVL2 = 2 * LVL - 1;
562: // *
563: // * find the first node LF and last node LL on
564: // * the current level LVL
565: // *
566: if (LVL == 1)
567: {
568: LF = 1;
569: LL = 1;
570: }
571: else
572: {
573: LF = (int)Math.Pow(2, LVL - 1);
574: LL = 2 * LF - 1;
575: }
576: for (I = LF; I <= LL; I++)
577: {
578: IM1 = I - 1;
579: IC = IWORK[INODE + IM1 + o_iwork];
580: NL = IWORK[NDIML + IM1 + o_iwork];
581: NR = IWORK[NDIMR + IM1 + o_iwork];
582: NLF = IC - NL;
583: NRF = IC + 1;
584: J = J - 1;
585: this._dlals0.Run(ICOMPQ, NL, NR, SQRE, NRHS, ref BX, NLF+1 * LDBX + o_bx
586: , LDBX, ref B, NLF+1 * LDB + o_b, LDB, PERM, NLF+LVL * LDGCOL + o_perm, GIVPTR[J + o_givptr], GIVCOL, NLF+LVL2 * LDGCOL + o_givcol
587: , LDGCOL, GIVNUM, NLF+LVL2 * LDU + o_givnum, LDU, POLES, NLF+LVL2 * LDU + o_poles, DIFL, NLF+LVL * LDU + o_difl, DIFR, NLF+LVL2 * LDU + o_difr
588: , Z, NLF+LVL * LDU + o_z, K[J + o_k], C[J + o_c], S[J + o_s], ref WORK, offset_work, ref INFO);
589: }
590: }
591: goto LABEL90;
592: // *
593: // * ICOMPQ = 1: applying back the right singular vector factors.
594: // *
595: LABEL50:;
596: // *
597: // * First now go through the right singular vector matrices of all
598: // * the tree nodes top-down.
599: // *
600: J = 0;
601: for (LVL = 1; LVL <= NLVL; LVL++)
602: {
603: LVL2 = 2 * LVL - 1;
604: // *
605: // * Find the first node LF and last node LL on
606: // * the current level LVL.
607: // *
608: if (LVL == 1)
609: {
610: LF = 1;
611: LL = 1;
612: }
613: else
614: {
615: LF = (int)Math.Pow(2, LVL - 1);
616: LL = 2 * LF - 1;
617: }
618: for (I = LL; I >= LF; I += - 1)
619: {
620: IM1 = I - 1;
621: IC = IWORK[INODE + IM1 + o_iwork];
622: NL = IWORK[NDIML + IM1 + o_iwork];
623: NR = IWORK[NDIMR + IM1 + o_iwork];
624: NLF = IC - NL;
625: NRF = IC + 1;
626: if (I == LL)
627: {
628: SQRE = 0;
629: }
630: else
631: {
632: SQRE = 1;
633: }
634: J = J + 1;
635: this._dlals0.Run(ICOMPQ, NL, NR, SQRE, NRHS, ref B, NLF+1 * LDB + o_b
636: , LDB, ref BX, NLF+1 * LDBX + o_bx, LDBX, PERM, NLF+LVL * LDGCOL + o_perm, GIVPTR[J + o_givptr], GIVCOL, NLF+LVL2 * LDGCOL + o_givcol
637: , LDGCOL, GIVNUM, NLF+LVL2 * LDU + o_givnum, LDU, POLES, NLF+LVL2 * LDU + o_poles, DIFL, NLF+LVL * LDU + o_difl, DIFR, NLF+LVL2 * LDU + o_difr
638: , Z, NLF+LVL * LDU + o_z, K[J + o_k], C[J + o_c], S[J + o_s], ref WORK, offset_work, ref INFO);
639: }
640: }
641: // *
642: // * The nodes on the bottom level of the tree were solved
643: // * by DLASDQ. The corresponding right singular vector
644: // * matrices are in explicit form. Apply them back.
645: // *
646: NDB1 = (ND + 1) / 2;
647: for (I = NDB1; I <= ND; I++)
648: {
649: I1 = I - 1;
650: IC = IWORK[INODE + I1 + o_iwork];
651: NL = IWORK[NDIML + I1 + o_iwork];
652: NR = IWORK[NDIMR + I1 + o_iwork];
653: NLP1 = NL + 1;
654: if (I == ND)
655: {
656: NRP1 = NR;
657: }
658: else
659: {
660: NRP1 = NR + 1;
661: }
662: NLF = IC - NL;
663: NRF = IC + 1;
664: this._dgemm.Run("T", "N", NLP1, NRHS, NLP1, ONE
665: , VT, NLF+1 * LDU + o_vt, LDU, B, NLF+1 * LDB + o_b, LDB, ZERO, ref BX, NLF+1 * LDBX + o_bx
666: , LDBX);
667: this._dgemm.Run("T", "N", NRP1, NRHS, NRP1, ONE
668: , VT, NRF+1 * LDU + o_vt, LDU, B, NRF+1 * LDB + o_b, LDB, ZERO, ref BX, NRF+1 * LDBX + o_bx
669: , LDBX);
670: }
671: // *
672: LABEL90:;
673: // *
674: return;
675: // *
676: // * End of DLALSA
677: // *
678:
679: #endregion
680:
681: }
682: }
683: }