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: /// DLALSD uses the singular value decomposition of A to solve the least
27: /// squares problem of finding X to minimize the Euclidean norm of each
28: /// column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
29: /// are N-by-NRHS. The solution X overwrites B.
30: ///
31: /// The singular values of A smaller than RCOND times the largest
32: /// singular value are treated as zero in solving the least squares
33: /// problem; in this case a minimum norm solution is returned.
34: /// The actual singular values are returned in D in ascending order.
35: ///
36: /// This code makes very mild assumptions about floating point
37: /// arithmetic. It will work on machines with a guard digit in
38: /// add/subtract, or on those binary machines without guard digits
39: /// which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
40: /// It could conceivably fail on hexadecimal or decimal machines
41: /// without guard digits, but we know of none.
42: ///
43: ///</summary>
44: public class DLALSD
45: {
46:
47:
48: #region Dependencies
49:
50: IDAMAX _idamax; DLAMCH _dlamch; DLANST _dlanst; DCOPY _dcopy; DGEMM _dgemm; DLACPY _dlacpy; DLALSA _dlalsa;
51: DLARTG _dlartg;DLASCL _dlascl; DLASDA _dlasda; DLASDQ _dlasdq; DLASET _dlaset; DLASRT _dlasrt; DROT _drot; XERBLA _xerbla;
52:
53: #endregion
54:
55:
56: #region Fields
57:
58: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; int BX = 0; int BXST = 0; int C = 0;
59: int DIFL = 0;int DIFR = 0; int GIVCOL = 0; int GIVNUM = 0; int GIVPTR = 0; int I = 0; int ICMPQ1 = 0; int ICMPQ2 = 0;
60: int IWK = 0;int J = 0; int K = 0; int NLVL = 0; int NM1 = 0; int NSIZE = 0; int NSUB = 0; int NWORK = 0; int PERM = 0;
61: int POLES = 0;int S = 0; int SIZEI = 0; int SMLSZP = 0; int SQRE = 0; int ST = 0; int ST1 = 0; int U = 0; int VT = 0;
62: int Z = 0;double CS = 0; double EPS = 0; double ORGNRM = 0; double R = 0; double RCND = 0; double SN = 0; double TOL = 0;
63:
64: #endregion
65:
66: public DLALSD(IDAMAX idamax, DLAMCH dlamch, DLANST dlanst, DCOPY dcopy, DGEMM dgemm, DLACPY dlacpy, DLALSA dlalsa, DLARTG dlartg, DLASCL dlascl, DLASDA dlasda
67: , DLASDQ dlasdq, DLASET dlaset, DLASRT dlasrt, DROT drot, XERBLA xerbla)
68: {
69:
70:
71: #region Set Dependencies
72:
73: this._idamax = idamax; this._dlamch = dlamch; this._dlanst = dlanst; this._dcopy = dcopy; this._dgemm = dgemm;
74: this._dlacpy = dlacpy;this._dlalsa = dlalsa; this._dlartg = dlartg; this._dlascl = dlascl; this._dlasda = dlasda;
75: this._dlasdq = dlasdq;this._dlaset = dlaset; this._dlasrt = dlasrt; this._drot = drot; this._xerbla = xerbla;
76:
77: #endregion
78:
79: }
80:
81: public DLALSD()
82: {
83:
84:
85: #region Dependencies (Initialization)
86:
87: IDAMAX idamax = new IDAMAX();
88: LSAME lsame = new LSAME();
89: DLAMC3 dlamc3 = new DLAMC3();
90: DLASSQ dlassq = new DLASSQ();
91: DCOPY dcopy = new DCOPY();
92: XERBLA xerbla = new XERBLA();
93: DROT drot = new DROT();
94: DSCAL dscal = new DSCAL();
95: DNRM2 dnrm2 = new DNRM2();
96: DLASDT dlasdt = new DLASDT();
97: DLAMRG dlamrg = new DLAMRG();
98: DLAPY2 dlapy2 = new DLAPY2();
99: DLASD5 dlasd5 = new DLASD5();
100: DDOT ddot = new DDOT();
101: DLAS2 dlas2 = new DLAS2();
102: DLASQ5 dlasq5 = new DLASQ5();
103: DLAZQ4 dlazq4 = new DLAZQ4();
104: IEEECK ieeeck = new IEEECK();
105: IPARMQ iparmq = new IPARMQ();
106: DSWAP dswap = new DSWAP();
107: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
108: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
109: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
110: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
111: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
112: DLANST dlanst = new DLANST(lsame, dlassq);
113: DGEMM dgemm = new DGEMM(lsame, xerbla);
114: DLACPY dlacpy = new DLACPY(lsame);
115: DGEMV dgemv = new DGEMV(lsame, xerbla);
116: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
117: DLALS0 dlals0 = new DLALS0(dcopy, dgemv, dlacpy, dlascl, drot, dscal, xerbla, dlamc3, dnrm2);
118: DLALSA dlalsa = new DLALSA(dcopy, dgemm, dlals0, dlasdt, xerbla);
119: DLARTG dlartg = new DLARTG(dlamch);
120: DLASD7 dlasd7 = new DLASD7(dcopy, dlamrg, drot, xerbla, dlamch, dlapy2);
121: DLAED6 dlaed6 = new DLAED6(dlamch);
122: DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
123: DLASET dlaset = new DLASET(lsame);
124: DLASD8 dlasd8 = new DLASD8(dcopy, dlascl, dlasd4, dlaset, xerbla, ddot, dlamc3, dnrm2);
125: DLASD6 dlasd6 = new DLASD6(dcopy, dlamrg, dlascl, dlasd7, dlasd8, xerbla);
126: DLASQ6 dlasq6 = new DLASQ6(dlamch);
127: DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
128: DLASRT dlasrt = new DLASRT(lsame, xerbla);
129: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
130: DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
131: DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
132: DLASR dlasr = new DLASR(lsame, xerbla);
133: DLASV2 dlasv2 = new DLASV2(dlamch);
134: DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap
135: , xerbla);
136: DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);
137: DLASDA dlasda = new DLASDA(dcopy, dlasd6, dlasdq, dlasdt, dlaset, xerbla);
138:
139: #endregion
140:
141:
142: #region Set Dependencies
143:
144: this._idamax = idamax; this._dlamch = dlamch; this._dlanst = dlanst; this._dcopy = dcopy; this._dgemm = dgemm;
145: this._dlacpy = dlacpy;this._dlalsa = dlalsa; this._dlartg = dlartg; this._dlascl = dlascl; this._dlasda = dlasda;
146: this._dlasdq = dlasdq;this._dlaset = dlaset; this._dlasrt = dlasrt; this._drot = drot; this._xerbla = xerbla;
147:
148: #endregion
149:
150: }
151: /// <summary>
152: /// Purpose
153: /// =======
154: ///
155: /// DLALSD uses the singular value decomposition of A to solve the least
156: /// squares problem of finding X to minimize the Euclidean norm of each
157: /// column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
158: /// are N-by-NRHS. The solution X overwrites B.
159: ///
160: /// The singular values of A smaller than RCOND times the largest
161: /// singular value are treated as zero in solving the least squares
162: /// problem; in this case a minimum norm solution is returned.
163: /// The actual singular values are returned in D in ascending order.
164: ///
165: /// This code makes very mild assumptions about floating point
166: /// arithmetic. It will work on machines with a guard digit in
167: /// add/subtract, or on those binary machines without guard digits
168: /// which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
169: /// It could conceivably fail on hexadecimal or decimal machines
170: /// without guard digits, but we know of none.
171: ///
172: ///</summary>
173: /// <param name="UPLO">
174: /// (input) CHARACTER*1
175: /// = 'U': D and E define an upper bidiagonal matrix.
176: /// = 'L': D and E define a lower bidiagonal matrix.
177: ///</param>
178: /// <param name="SMLSIZ">
179: /// (input) INTEGER
180: /// The maximum size of the subproblems at the bottom of the
181: /// computation tree.
182: ///</param>
183: /// <param name="N">
184: /// (input) INTEGER
185: /// The dimension of the bidiagonal matrix. N .GE. 0.
186: ///</param>
187: /// <param name="NRHS">
188: /// (input) INTEGER
189: /// The number of columns of B. NRHS must be at least 1.
190: ///</param>
191: /// <param name="D">
192: /// (input/output) DOUBLE PRECISION array, dimension (N)
193: /// On entry D contains the main diagonal of the bidiagonal
194: /// matrix. On exit, if INFO = 0, D contains its singular values.
195: ///</param>
196: /// <param name="E">
197: /// (input/output) DOUBLE PRECISION array, dimension (N-1)
198: /// Contains the super-diagonal entries of the bidiagonal matrix.
199: /// On exit, E has been destroyed.
200: ///</param>
201: /// <param name="B">
202: /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
203: /// On input, B contains the right hand sides of the least
204: /// squares problem. On output, B contains the solution X.
205: ///</param>
206: /// <param name="LDB">
207: /// (input) INTEGER
208: /// The leading dimension of B in the calling subprogram.
209: /// LDB must be at least max(1,N).
210: ///</param>
211: /// <param name="RCOND">
212: /// (input) DOUBLE PRECISION
213: /// The singular values of A less than or equal to RCOND times
214: /// the largest singular value are treated as zero in solving
215: /// the least squares problem. If RCOND is negative,
216: /// machine precision is used instead.
217: /// For example, if diag(S)*X=B were the least squares problem,
218: /// where diag(S) is a diagonal matrix of singular values, the
219: /// solution would be X(i) = B(i) / S(i) if S(i) is greater than
220: /// RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
221: /// RCOND*max(S).
222: ///</param>
223: /// <param name="RANK">
224: /// (output) INTEGER
225: /// The number of singular values of A greater than RCOND times
226: /// the largest singular value.
227: ///</param>
228: /// <param name="WORK">
229: /// (workspace) DOUBLE PRECISION array, dimension at least
230: /// (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
231: /// where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
232: ///</param>
233: /// <param name="IWORK">
234: /// (workspace) INTEGER array, dimension at least
235: /// (3*N*NLVL + 11*N)
236: ///</param>
237: /// <param name="INFO">
238: /// (output) INTEGER
239: /// = 0: successful exit.
240: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
241: /// .GT. 0: The algorithm failed to compute an singular value while
242: /// working on the submatrix lying in rows and columns
243: /// INFO/(N+1) through MOD(INFO,N+1).
244: ///</param>
245: public void Run(string UPLO, int SMLSIZ, int N, int NRHS, ref double[] D, int offset_d, ref double[] E, int offset_e
246: , ref double[] B, int offset_b, int LDB, double RCOND, ref int RANK, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork
247: , ref int INFO)
248: {
249:
250: #region Array Index Correction
251:
252: int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_b = -1 - LDB + offset_b; int o_work = -1 + offset_work;
253: int o_iwork = -1 + offset_iwork;
254:
255: #endregion
256:
257:
258: #region Strings
259:
260: UPLO = UPLO.Substring(0, 1);
261:
262: #endregion
263:
264:
265: #region Prolog
266:
267: // *
268: // * -- LAPACK routine (version 3.1) --
269: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
270: // * November 2006
271: // *
272: // * .. Scalar Arguments ..
273: // * ..
274: // * .. Array Arguments ..
275: // * ..
276: // *
277: // * Purpose
278: // * =======
279: // *
280: // * DLALSD uses the singular value decomposition of A to solve the least
281: // * squares problem of finding X to minimize the Euclidean norm of each
282: // * column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
283: // * are N-by-NRHS. The solution X overwrites B.
284: // *
285: // * The singular values of A smaller than RCOND times the largest
286: // * singular value are treated as zero in solving the least squares
287: // * problem; in this case a minimum norm solution is returned.
288: // * The actual singular values are returned in D in ascending order.
289: // *
290: // * This code makes very mild assumptions about floating point
291: // * arithmetic. It will work on machines with a guard digit in
292: // * add/subtract, or on those binary machines without guard digits
293: // * which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
294: // * It could conceivably fail on hexadecimal or decimal machines
295: // * without guard digits, but we know of none.
296: // *
297: // * Arguments
298: // * =========
299: // *
300: // * UPLO (input) CHARACTER*1
301: // * = 'U': D and E define an upper bidiagonal matrix.
302: // * = 'L': D and E define a lower bidiagonal matrix.
303: // *
304: // * SMLSIZ (input) INTEGER
305: // * The maximum size of the subproblems at the bottom of the
306: // * computation tree.
307: // *
308: // * N (input) INTEGER
309: // * The dimension of the bidiagonal matrix. N >= 0.
310: // *
311: // * NRHS (input) INTEGER
312: // * The number of columns of B. NRHS must be at least 1.
313: // *
314: // * D (input/output) DOUBLE PRECISION array, dimension (N)
315: // * On entry D contains the main diagonal of the bidiagonal
316: // * matrix. On exit, if INFO = 0, D contains its singular values.
317: // *
318: // * E (input/output) DOUBLE PRECISION array, dimension (N-1)
319: // * Contains the super-diagonal entries of the bidiagonal matrix.
320: // * On exit, E has been destroyed.
321: // *
322: // * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
323: // * On input, B contains the right hand sides of the least
324: // * squares problem. On output, B contains the solution X.
325: // *
326: // * LDB (input) INTEGER
327: // * The leading dimension of B in the calling subprogram.
328: // * LDB must be at least max(1,N).
329: // *
330: // * RCOND (input) DOUBLE PRECISION
331: // * The singular values of A less than or equal to RCOND times
332: // * the largest singular value are treated as zero in solving
333: // * the least squares problem. If RCOND is negative,
334: // * machine precision is used instead.
335: // * For example, if diag(S)*X=B were the least squares problem,
336: // * where diag(S) is a diagonal matrix of singular values, the
337: // * solution would be X(i) = B(i) / S(i) if S(i) is greater than
338: // * RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
339: // * RCOND*max(S).
340: // *
341: // * RANK (output) INTEGER
342: // * The number of singular values of A greater than RCOND times
343: // * the largest singular value.
344: // *
345: // * WORK (workspace) DOUBLE PRECISION array, dimension at least
346: // * (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
347: // * where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
348: // *
349: // * IWORK (workspace) INTEGER array, dimension at least
350: // * (3*N*NLVL + 11*N)
351: // *
352: // * INFO (output) INTEGER
353: // * = 0: successful exit.
354: // * < 0: if INFO = -i, the i-th argument had an illegal value.
355: // * > 0: The algorithm failed to compute an singular value while
356: // * working on the submatrix lying in rows and columns
357: // * INFO/(N+1) through MOD(INFO,N+1).
358: // *
359: // * Further Details
360: // * ===============
361: // *
362: // * Based on contributions by
363: // * Ming Gu and Ren-Cang Li, Computer Science Division, University of
364: // * California at Berkeley, USA
365: // * Osni Marques, LBNL/NERSC, USA
366: // *
367: // * =====================================================================
368: // *
369: // * .. Parameters ..
370: // * ..
371: // * .. Local Scalars ..
372: // * ..
373: // * .. External Functions ..
374: // * ..
375: // * .. External Subroutines ..
376: // * ..
377: // * .. Intrinsic Functions ..
378: // INTRINSIC ABS, DBLE, INT, LOG, SIGN;
379: // * ..
380: // * .. Executable Statements ..
381: // *
382: // * Test the input parameters.
383: // *
384:
385: #endregion
386:
387:
388: #region Body
389:
390: INFO = 0;
391: // *
392: if (N < 0)
393: {
394: INFO = - 3;
395: }
396: else
397: {
398: if (NRHS < 1)
399: {
400: INFO = - 4;
401: }
402: else
403: {
404: if ((LDB < 1) || (LDB < N))
405: {
406: INFO = - 8;
407: }
408: }
409: }
410: if (INFO != 0)
411: {
412: this._xerbla.Run("DLALSD", - INFO);
413: return;
414: }
415: // *
416: EPS = this._dlamch.Run("Epsilon");
417: // *
418: // * Set up the tolerance.
419: // *
420: if ((RCOND <= ZERO) || (RCOND >= ONE))
421: {
422: RCND = EPS;
423: }
424: else
425: {
426: RCND = RCOND;
427: }
428: // *
429: RANK = 0;
430: // *
431: // * Quick return if possible.
432: // *
433: if (N == 0)
434: {
435: return;
436: }
437: else
438: {
439: if (N == 1)
440: {
441: if (D[1 + o_d] == ZERO)
442: {
443: this._dlaset.Run("A", 1, NRHS, ZERO, ZERO, ref B, offset_b
444: , LDB);
445: }
446: else
447: {
448: RANK = 1;
449: this._dlascl.Run("G", 0, 0, D[1 + o_d], ONE, 1
450: , NRHS, ref B, offset_b, LDB, ref INFO);
451: D[1 + o_d] = Math.Abs(D[1 + o_d]);
452: }
453: return;
454: }
455: }
456: // *
457: // * Rotate the matrix if it is lower bidiagonal.
458: // *
459: if (UPLO == "L")
460: {
461: for (I = 1; I <= N - 1; I++)
462: {
463: this._dlartg.Run(D[I + o_d], E[I + o_e], ref CS, ref SN, ref R);
464: D[I + o_d] = R;
465: E[I + o_e] = SN * D[I + 1 + o_d];
466: D[I + 1 + o_d] = CS * D[I + 1 + o_d];
467: if (NRHS == 1)
468: {
469: this._drot.Run(1, ref B, I+1 * LDB + o_b, 1, ref B, I + 1+1 * LDB + o_b, 1, CS
470: , SN);
471: }
472: else
473: {
474: WORK[I * 2 - 1 + o_work] = CS;
475: WORK[I * 2 + o_work] = SN;
476: }
477: }
478: if (NRHS > 1)
479: {
480: for (I = 1; I <= NRHS; I++)
481: {
482: for (J = 1; J <= N - 1; J++)
483: {
484: CS = WORK[J * 2 - 1 + o_work];
485: SN = WORK[J * 2 + o_work];
486: this._drot.Run(1, ref B, J+I * LDB + o_b, 1, ref B, J + 1+I * LDB + o_b, 1, CS
487: , SN);
488: }
489: }
490: }
491: }
492: // *
493: // * Scale.
494: // *
495: NM1 = N - 1;
496: ORGNRM = this._dlanst.Run("M", N, D, offset_d, E, offset_e);
497: if (ORGNRM == ZERO)
498: {
499: this._dlaset.Run("A", N, NRHS, ZERO, ZERO, ref B, offset_b
500: , LDB);
501: return;
502: }
503: // *
504: this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N
505: , 1, ref D, offset_d, N, ref INFO);
506: this._dlascl.Run("G", 0, 0, ORGNRM, ONE, NM1
507: , 1, ref E, offset_e, NM1, ref INFO);
508: // *
509: // * If N is smaller than the minimum divide size SMLSIZ, then solve
510: // * the problem with another solver.
511: // *
512: if (N <= SMLSIZ)
513: {
514: NWORK = 1 + N * N;
515: this._dlaset.Run("A", N, N, ZERO, ONE, ref WORK, offset_work
516: , N);
517: this._dlasdq.Run("U", 0, N, N, 0, NRHS
518: , ref D, offset_d, ref E, offset_e, ref WORK, offset_work, N, ref WORK, offset_work, N
519: , ref B, offset_b, LDB, ref WORK, NWORK + o_work, ref INFO);
520: if (INFO != 0)
521: {
522: return;
523: }
524: TOL = RCND * Math.Abs(D[this._idamax.Run(N, D, offset_d, 1) + o_d]);
525: for (I = 1; I <= N; I++)
526: {
527: if (D[I + o_d] <= TOL)
528: {
529: this._dlaset.Run("A", 1, NRHS, ZERO, ZERO, ref B, I+1 * LDB + o_b
530: , LDB);
531: }
532: else
533: {
534: this._dlascl.Run("G", 0, 0, D[I + o_d], ONE, 1
535: , NRHS, ref B, I+1 * LDB + o_b, LDB, ref INFO);
536: RANK = RANK + 1;
537: }
538: }
539: this._dgemm.Run("T", "N", N, NRHS, N, ONE
540: , WORK, offset_work, N, B, offset_b, LDB, ZERO, ref WORK, NWORK + o_work
541: , N);
542: this._dlacpy.Run("A", N, NRHS, WORK, NWORK + o_work, N, ref B, offset_b
543: , LDB);
544: // *
545: // * Unscale.
546: // *
547: this._dlascl.Run("G", 0, 0, ONE, ORGNRM, N
548: , 1, ref D, offset_d, N, ref INFO);
549: this._dlasrt.Run("D", N, ref D, offset_d, ref INFO);
550: this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N
551: , NRHS, ref B, offset_b, LDB, ref INFO);
552: // *
553: return;
554: }
555: // *
556: // * Book-keeping and setting up some constants.
557: // *
558: NLVL = Convert.ToInt32(Math.Truncate(Math.Log(Convert.ToDouble(N) / Convert.ToDouble(SMLSIZ + 1)) / Math.Log(TWO))) + 1;
559: // *
560: SMLSZP = SMLSIZ + 1;
561: // *
562: U = 1;
563: VT = 1 + SMLSIZ * N;
564: DIFL = VT + SMLSZP * N;
565: DIFR = DIFL + NLVL * N;
566: Z = DIFR + NLVL * N * 2;
567: C = Z + NLVL * N;
568: S = C + N;
569: POLES = S + N;
570: GIVNUM = POLES + 2 * NLVL * N;
571: BX = GIVNUM + 2 * NLVL * N;
572: NWORK = BX + N * NRHS;
573: // *
574: SIZEI = 1 + N;
575: K = SIZEI + N;
576: GIVPTR = K + N;
577: PERM = GIVPTR + N;
578: GIVCOL = PERM + NLVL * N;
579: IWK = GIVCOL + NLVL * N * 2;
580: // *
581: ST = 1;
582: SQRE = 0;
583: ICMPQ1 = 1;
584: ICMPQ2 = 0;
585: NSUB = 0;
586: // *
587: for (I = 1; I <= N; I++)
588: {
589: if (Math.Abs(D[I + o_d]) < EPS)
590: {
591: D[I + o_d] = FortranLib.Sign(EPS,D[I + o_d]);
592: }
593: }
594: // *
595: for (I = 1; I <= NM1; I++)
596: {
597: if ((Math.Abs(E[I + o_e]) < EPS) || (I == NM1))
598: {
599: NSUB = NSUB + 1;
600: IWORK[NSUB + o_iwork] = ST;
601: // *
602: // * Subproblem found. First determine its size and then
603: // * apply divide and conquer on it.
604: // *
605: if (I < NM1)
606: {
607: // *
608: // * A subproblem with E(I) small for I < NM1.
609: // *
610: NSIZE = I - ST + 1;
611: IWORK[SIZEI + NSUB - 1 + o_iwork] = NSIZE;
612: }
613: else
614: {
615: if (Math.Abs(E[I + o_e]) >= EPS)
616: {
617: // *
618: // * A subproblem with E(NM1) not too small but I = NM1.
619: // *
620: NSIZE = N - ST + 1;
621: IWORK[SIZEI + NSUB - 1 + o_iwork] = NSIZE;
622: }
623: else
624: {
625: // *
626: // * A subproblem with E(NM1) small. This implies an
627: // * 1-by-1 subproblem at D(N), which is not solved
628: // * explicitly.
629: // *
630: NSIZE = I - ST + 1;
631: IWORK[SIZEI + NSUB - 1 + o_iwork] = NSIZE;
632: NSUB = NSUB + 1;
633: IWORK[NSUB + o_iwork] = N;
634: IWORK[SIZEI + NSUB - 1 + o_iwork] = 1;
635: this._dcopy.Run(NRHS, B, N+1 * LDB + o_b, LDB, ref WORK, BX + NM1 + o_work, N);
636: }
637: }
638: ST1 = ST - 1;
639: if (NSIZE == 1)
640: {
641: // *
642: // * This is a 1-by-1 subproblem and is not solved
643: // * explicitly.
644: // *
645: this._dcopy.Run(NRHS, B, ST+1 * LDB + o_b, LDB, ref WORK, BX + ST1 + o_work, N);
646: }
647: else
648: {
649: if (NSIZE <= SMLSIZ)
650: {
651: // *
652: // * This is a small subproblem and is solved by DLASDQ.
653: // *
654: this._dlaset.Run("A", NSIZE, NSIZE, ZERO, ONE, ref WORK, VT + ST1 + o_work
655: , N);
656: this._dlasdq.Run("U", 0, NSIZE, NSIZE, 0, NRHS
657: , ref D, ST + o_d, ref E, ST + o_e, ref WORK, VT + ST1 + o_work, N, ref WORK, NWORK + o_work, N
658: , ref B, ST+1 * LDB + o_b, LDB, ref WORK, NWORK + o_work, ref INFO);
659: if (INFO != 0)
660: {
661: return;
662: }
663: this._dlacpy.Run("A", NSIZE, NRHS, B, ST+1 * LDB + o_b, LDB, ref WORK, BX + ST1 + o_work
664: , N);
665: }
666: else
667: {
668: // *
669: // * A large problem. Solve it using divide and conquer.
670: // *
671: this._dlasda.Run(ICMPQ1, SMLSIZ, NSIZE, SQRE, ref D, ST + o_d, ref E, ST + o_e
672: , ref WORK, U + ST1 + o_work, N, ref WORK, VT + ST1 + o_work, ref IWORK, K + ST1 + o_iwork, ref WORK, DIFL + ST1 + o_work, ref WORK, DIFR + ST1 + o_work
673: , ref WORK, Z + ST1 + o_work, ref WORK, POLES + ST1 + o_work, ref IWORK, GIVPTR + ST1 + o_iwork, ref IWORK, GIVCOL + ST1 + o_iwork, N, ref IWORK, PERM + ST1 + o_iwork
674: , ref WORK, GIVNUM + ST1 + o_work, ref WORK, C + ST1 + o_work, ref WORK, S + ST1 + o_work, ref WORK, NWORK + o_work, ref IWORK, IWK + o_iwork, ref INFO);
675: if (INFO != 0)
676: {
677: return;
678: }
679: BXST = BX + ST1;
680: this._dlalsa.Run(ICMPQ2, SMLSIZ, NSIZE, NRHS, ref B, ST+1 * LDB + o_b, LDB
681: , ref WORK, BXST + o_work, N, WORK, U + ST1 + o_work, N, WORK, VT + ST1 + o_work, IWORK, K + ST1 + o_iwork
682: , WORK, DIFL + ST1 + o_work, WORK, DIFR + ST1 + o_work, WORK, Z + ST1 + o_work, WORK, POLES + ST1 + o_work, IWORK, GIVPTR + ST1 + o_iwork, IWORK, GIVCOL + ST1 + o_iwork
683: , N, IWORK, PERM + ST1 + o_iwork, WORK, GIVNUM + ST1 + o_work, WORK, C + ST1 + o_work, WORK, S + ST1 + o_work, ref WORK, NWORK + o_work
684: , ref IWORK, IWK + o_iwork, ref INFO);
685: if (INFO != 0)
686: {
687: return;
688: }
689: }
690: }
691: ST = I + 1;
692: }
693: }
694: // *
695: // * Apply the singular values and treat the tiny ones as zero.
696: // *
697: TOL = RCND * Math.Abs(D[this._idamax.Run(N, D, offset_d, 1) + o_d]);
698: // *
699: for (I = 1; I <= N; I++)
700: {
701: // *
702: // * Some of the elements in D can be negative because 1-by-1
703: // * subproblems were not solved explicitly.
704: // *
705: if (Math.Abs(D[I + o_d]) <= TOL)
706: {
707: this._dlaset.Run("A", 1, NRHS, ZERO, ZERO, ref WORK, BX + I - 1 + o_work
708: , N);
709: }
710: else
711: {
712: RANK = RANK + 1;
713: this._dlascl.Run("G", 0, 0, D[I + o_d], ONE, 1
714: , NRHS, ref WORK, BX + I - 1 + o_work, N, ref INFO);
715: }
716: D[I + o_d] = Math.Abs(D[I + o_d]);
717: }
718: // *
719: // * Now apply back the right singular vectors.
720: // *
721: ICMPQ2 = 1;
722: for (I = 1; I <= NSUB; I++)
723: {
724: ST = IWORK[I + o_iwork];
725: ST1 = ST - 1;
726: NSIZE = IWORK[SIZEI + I - 1 + o_iwork];
727: BXST = BX + ST1;
728: if (NSIZE == 1)
729: {
730: this._dcopy.Run(NRHS, WORK, BXST + o_work, N, ref B, ST+1 * LDB + o_b, LDB);
731: }
732: else
733: {
734: if (NSIZE <= SMLSIZ)
735: {
736: this._dgemm.Run("T", "N", NSIZE, NRHS, NSIZE, ONE
737: , WORK, VT + ST1 + o_work, N, WORK, BXST + o_work, N, ZERO, ref B, ST+1 * LDB + o_b
738: , LDB);
739: }
740: else
741: {
742: this._dlalsa.Run(ICMPQ2, SMLSIZ, NSIZE, NRHS, ref WORK, BXST + o_work, N
743: , ref B, ST+1 * LDB + o_b, LDB, WORK, U + ST1 + o_work, N, WORK, VT + ST1 + o_work, IWORK, K + ST1 + o_iwork
744: , WORK, DIFL + ST1 + o_work, WORK, DIFR + ST1 + o_work, WORK, Z + ST1 + o_work, WORK, POLES + ST1 + o_work, IWORK, GIVPTR + ST1 + o_iwork, IWORK, GIVCOL + ST1 + o_iwork
745: , N, IWORK, PERM + ST1 + o_iwork, WORK, GIVNUM + ST1 + o_work, WORK, C + ST1 + o_work, WORK, S + ST1 + o_work, ref WORK, NWORK + o_work
746: , ref IWORK, IWK + o_iwork, ref INFO);
747: if (INFO != 0)
748: {
749: return;
750: }
751: }
752: }
753: }
754: // *
755: // * Unscale and sort the singular values.
756: // *
757: this._dlascl.Run("G", 0, 0, ONE, ORGNRM, N
758: , 1, ref D, offset_d, N, ref INFO);
759: this._dlasrt.Run("D", N, ref D, offset_d, ref INFO);
760: this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N
761: , NRHS, ref B, offset_b, LDB, ref INFO);
762: // *
763: return;
764: // *
765: // * End of DLALSD
766: // *
767:
768: #endregion
769:
770: }
771: }
772: }