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: /// DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
27: /// where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
28: ///
29: /// A related subroutine DLASD7 handles the case in which the singular
30: /// values (and the singular vectors in factored form) are desired.
31: ///
32: /// DLASD1 computes the SVD as follows:
33: ///
34: /// ( D1(in) 0 0 0 )
35: /// B = U(in) * ( Z1' a Z2' b ) * VT(in)
36: /// ( 0 0 D2(in) 0 )
37: ///
38: /// = U(out) * ( D(out) 0) * VT(out)
39: ///
40: /// where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
41: /// with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
42: /// elsewhere; and the entry b is empty if SQRE = 0.
43: ///
44: /// The left singular vectors of the original matrix are stored in U, and
45: /// the transpose of the right singular vectors are stored in VT, and the
46: /// singular values are in D. The algorithm consists of three stages:
47: ///
48: /// The first stage consists of deflating the size of the problem
49: /// when there are multiple singular values or when there are zeros in
50: /// the Z vector. For each such occurence the dimension of the
51: /// secular equation problem is reduced by one. This stage is
52: /// performed by the routine DLASD2.
53: ///
54: /// The second stage consists of calculating the updated
55: /// singular values. This is done by finding the square roots of the
56: /// roots of the secular equation via the routine DLASD4 (as called
57: /// by DLASD3). This routine also calculates the singular vectors of
58: /// the current problem.
59: ///
60: /// The final stage consists of computing the updated singular vectors
61: /// directly using the updated singular values. The singular vectors
62: /// for the current problem are multiplied with the singular vectors
63: /// from the overall problem.
64: ///
65: ///</summary>
66: public class DLASD1
67: {
68:
69:
70: #region Dependencies
71:
72: DLAMRG _dlamrg; DLASCL _dlascl; DLASD2 _dlasd2; DLASD3 _dlasd3; XERBLA _xerbla;
73:
74: #endregion
75:
76:
77: #region Fields
78:
79: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; int COLTYP = 0; int I = 0; int IDX = 0; int IDXC = 0; int IDXP = 0;
80: int IQ = 0;int ISIGMA = 0; int IU2 = 0; int IVT2 = 0; int IZ = 0; int K = 0; int LDQ = 0; int LDU2 = 0; int LDVT2 = 0;
81: int M = 0;int N = 0; int N1 = 0; int N2 = 0; double ORGNRM = 0;
82:
83: #endregion
84:
85: public DLASD1(DLAMRG dlamrg, DLASCL dlascl, DLASD2 dlasd2, DLASD3 dlasd3, XERBLA xerbla)
86: {
87:
88:
89: #region Set Dependencies
90:
91: this._dlamrg = dlamrg; this._dlascl = dlascl; this._dlasd2 = dlasd2; this._dlasd3 = dlasd3; this._xerbla = xerbla;
92:
93: #endregion
94:
95: }
96:
97: public DLASD1()
98: {
99:
100:
101: #region Dependencies (Initialization)
102:
103: DLAMRG dlamrg = new DLAMRG();
104: LSAME lsame = new LSAME();
105: DLAMC3 dlamc3 = new DLAMC3();
106: XERBLA xerbla = new XERBLA();
107: DLAPY2 dlapy2 = new DLAPY2();
108: DCOPY dcopy = new DCOPY();
109: DROT drot = new DROT();
110: DNRM2 dnrm2 = new DNRM2();
111: DLASD5 dlasd5 = new DLASD5();
112: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
113: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
114: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
115: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
116: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
117: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
118: DLACPY dlacpy = new DLACPY(lsame);
119: DLASET dlaset = new DLASET(lsame);
120: DLASD2 dlasd2 = new DLASD2(dlamch, dlapy2, dcopy, dlacpy, dlamrg, dlaset, drot, xerbla);
121: DGEMM dgemm = new DGEMM(lsame, xerbla);
122: DLAED6 dlaed6 = new DLAED6(dlamch);
123: DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
124: DLASD3 dlasd3 = new DLASD3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlascl, dlasd4, xerbla);
125:
126: #endregion
127:
128:
129: #region Set Dependencies
130:
131: this._dlamrg = dlamrg; this._dlascl = dlascl; this._dlasd2 = dlasd2; this._dlasd3 = dlasd3; this._xerbla = xerbla;
132:
133: #endregion
134:
135: }
136: /// <summary>
137: /// Purpose
138: /// =======
139: ///
140: /// DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
141: /// where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
142: ///
143: /// A related subroutine DLASD7 handles the case in which the singular
144: /// values (and the singular vectors in factored form) are desired.
145: ///
146: /// DLASD1 computes the SVD as follows:
147: ///
148: /// ( D1(in) 0 0 0 )
149: /// B = U(in) * ( Z1' a Z2' b ) * VT(in)
150: /// ( 0 0 D2(in) 0 )
151: ///
152: /// = U(out) * ( D(out) 0) * VT(out)
153: ///
154: /// where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
155: /// with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
156: /// elsewhere; and the entry b is empty if SQRE = 0.
157: ///
158: /// The left singular vectors of the original matrix are stored in U, and
159: /// the transpose of the right singular vectors are stored in VT, and the
160: /// singular values are in D. The algorithm consists of three stages:
161: ///
162: /// The first stage consists of deflating the size of the problem
163: /// when there are multiple singular values or when there are zeros in
164: /// the Z vector. For each such occurence the dimension of the
165: /// secular equation problem is reduced by one. This stage is
166: /// performed by the routine DLASD2.
167: ///
168: /// The second stage consists of calculating the updated
169: /// singular values. This is done by finding the square roots of the
170: /// roots of the secular equation via the routine DLASD4 (as called
171: /// by DLASD3). This routine also calculates the singular vectors of
172: /// the current problem.
173: ///
174: /// The final stage consists of computing the updated singular vectors
175: /// directly using the updated singular values. The singular vectors
176: /// for the current problem are multiplied with the singular vectors
177: /// from the overall problem.
178: ///
179: ///</summary>
180: /// <param name="NL">
181: /// (input) INTEGER
182: /// The row dimension of the upper block. NL .GE. 1.
183: ///</param>
184: /// <param name="NR">
185: /// (input) INTEGER
186: /// The row dimension of the lower block. NR .GE. 1.
187: ///</param>
188: /// <param name="SQRE">
189: /// (input) INTEGER
190: /// = 0: the lower block is an NR-by-NR square matrix.
191: /// = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
192: ///
193: /// The bidiagonal matrix has row dimension N = NL + NR + 1,
194: /// and column dimension M = N + SQRE.
195: ///</param>
196: /// <param name="D">
197: /// (input/output) DOUBLE PRECISION array,
198: /// dimension (N = NL+NR+1).
199: /// On entry D(1:NL,1:NL) contains the singular values of the
200: /// upper block; and D(NL+2:N) contains the singular values of
201: /// the lower block. On exit D(1:N) contains the singular values
202: /// of the modified matrix.
203: ///</param>
204: /// <param name="ALPHA">
205: /// (input/output) DOUBLE PRECISION
206: /// Contains the diagonal element associated with the added row.
207: ///</param>
208: /// <param name="BETA">
209: /// (input/output) DOUBLE PRECISION
210: /// Contains the off-diagonal element associated with the added
211: /// row.
212: ///</param>
213: /// <param name="U">
214: /// (input/output) DOUBLE PRECISION array, dimension(LDU,N)
215: /// On entry U(1:NL, 1:NL) contains the left singular vectors of
216: /// the upper block; U(NL+2:N, NL+2:N) contains the left singular
217: /// vectors of the lower block. On exit U contains the left
218: /// singular vectors of the bidiagonal matrix.
219: ///</param>
220: /// <param name="LDU">
221: /// (input) INTEGER
222: /// The leading dimension of the array U. LDU .GE. max( 1, N ).
223: ///</param>
224: /// <param name="VT">
225: /// (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
226: /// where M = N + SQRE.
227: /// On entry VT(1:NL+1, 1:NL+1)' contains the right singular
228: /// vectors of the upper block; VT(NL+2:M, NL+2:M)' contains
229: /// the right singular vectors of the lower block. On exit
230: /// VT' contains the right singular vectors of the
231: /// bidiagonal matrix.
232: ///</param>
233: /// <param name="LDVT">
234: /// (input) INTEGER
235: /// The leading dimension of the array VT. LDVT .GE. max( 1, M ).
236: ///</param>
237: /// <param name="IDXQ">
238: /// (output) INTEGER array, dimension(N)
239: /// This contains the permutation which will reintegrate the
240: /// subproblem just solved back into sorted order, i.e.
241: /// D( IDXQ( I = 1, N ) ) will be in ascending order.
242: ///</param>
243: /// <param name="IWORK">
244: /// (workspace) INTEGER array, dimension( 4 * N )
245: ///</param>
246: /// <param name="WORK">
247: /// (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
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: /// .GT. 0: if INFO = 1, an singular value did not converge
254: ///</param>
255: public void Run(int NL, int NR, int SQRE, ref double[] D, int offset_d, ref double ALPHA, ref double BETA
256: , ref double[] U, int offset_u, int LDU, ref double[] VT, int offset_vt, int LDVT, ref int[] IDXQ, int offset_idxq, ref int[] IWORK, int offset_iwork
257: , ref double[] WORK, int offset_work, ref int INFO)
258: {
259:
260: #region Array Index Correction
261:
262: int o_d = -1 + offset_d; int o_u = -1 - LDU + offset_u; int o_vt = -1 - LDVT + offset_vt;
263: int o_idxq = -1 + offset_idxq; int o_iwork = -1 + offset_iwork; int o_work = -1 + offset_work;
264:
265: #endregion
266:
267:
268: #region Prolog
269:
270: // *
271: // * -- LAPACK auxiliary routine (version 3.1) --
272: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
273: // * November 2006
274: // *
275: // * .. Scalar Arguments ..
276: // * ..
277: // * .. Array Arguments ..
278: // * ..
279: // *
280: // * Purpose
281: // * =======
282: // *
283: // * DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
284: // * where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
285: // *
286: // * A related subroutine DLASD7 handles the case in which the singular
287: // * values (and the singular vectors in factored form) are desired.
288: // *
289: // * DLASD1 computes the SVD as follows:
290: // *
291: // * ( D1(in) 0 0 0 )
292: // * B = U(in) * ( Z1' a Z2' b ) * VT(in)
293: // * ( 0 0 D2(in) 0 )
294: // *
295: // * = U(out) * ( D(out) 0) * VT(out)
296: // *
297: // * where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
298: // * with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
299: // * elsewhere; and the entry b is empty if SQRE = 0.
300: // *
301: // * The left singular vectors of the original matrix are stored in U, and
302: // * the transpose of the right singular vectors are stored in VT, and the
303: // * singular values are in D. The algorithm consists of three stages:
304: // *
305: // * The first stage consists of deflating the size of the problem
306: // * when there are multiple singular values or when there are zeros in
307: // * the Z vector. For each such occurence the dimension of the
308: // * secular equation problem is reduced by one. This stage is
309: // * performed by the routine DLASD2.
310: // *
311: // * The second stage consists of calculating the updated
312: // * singular values. This is done by finding the square roots of the
313: // * roots of the secular equation via the routine DLASD4 (as called
314: // * by DLASD3). This routine also calculates the singular vectors of
315: // * the current problem.
316: // *
317: // * The final stage consists of computing the updated singular vectors
318: // * directly using the updated singular values. The singular vectors
319: // * for the current problem are multiplied with the singular vectors
320: // * from the overall problem.
321: // *
322: // * Arguments
323: // * =========
324: // *
325: // * NL (input) INTEGER
326: // * The row dimension of the upper block. NL >= 1.
327: // *
328: // * NR (input) INTEGER
329: // * The row dimension of the lower block. NR >= 1.
330: // *
331: // * SQRE (input) INTEGER
332: // * = 0: the lower block is an NR-by-NR square matrix.
333: // * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
334: // *
335: // * The bidiagonal matrix has row dimension N = NL + NR + 1,
336: // * and column dimension M = N + SQRE.
337: // *
338: // * D (input/output) DOUBLE PRECISION array,
339: // * dimension (N = NL+NR+1).
340: // * On entry D(1:NL,1:NL) contains the singular values of the
341: // * upper block; and D(NL+2:N) contains the singular values of
342: // * the lower block. On exit D(1:N) contains the singular values
343: // * of the modified matrix.
344: // *
345: // * ALPHA (input/output) DOUBLE PRECISION
346: // * Contains the diagonal element associated with the added row.
347: // *
348: // * BETA (input/output) DOUBLE PRECISION
349: // * Contains the off-diagonal element associated with the added
350: // * row.
351: // *
352: // * U (input/output) DOUBLE PRECISION array, dimension(LDU,N)
353: // * On entry U(1:NL, 1:NL) contains the left singular vectors of
354: // * the upper block; U(NL+2:N, NL+2:N) contains the left singular
355: // * vectors of the lower block. On exit U contains the left
356: // * singular vectors of the bidiagonal matrix.
357: // *
358: // * LDU (input) INTEGER
359: // * The leading dimension of the array U. LDU >= max( 1, N ).
360: // *
361: // * VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
362: // * where M = N + SQRE.
363: // * On entry VT(1:NL+1, 1:NL+1)' contains the right singular
364: // * vectors of the upper block; VT(NL+2:M, NL+2:M)' contains
365: // * the right singular vectors of the lower block. On exit
366: // * VT' contains the right singular vectors of the
367: // * bidiagonal matrix.
368: // *
369: // * LDVT (input) INTEGER
370: // * The leading dimension of the array VT. LDVT >= max( 1, M ).
371: // *
372: // * IDXQ (output) INTEGER array, dimension(N)
373: // * This contains the permutation which will reintegrate the
374: // * subproblem just solved back into sorted order, i.e.
375: // * D( IDXQ( I = 1, N ) ) will be in ascending order.
376: // *
377: // * IWORK (workspace) INTEGER array, dimension( 4 * N )
378: // *
379: // * WORK (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
380: // *
381: // * INFO (output) INTEGER
382: // * = 0: successful exit.
383: // * < 0: if INFO = -i, the i-th argument had an illegal value.
384: // * > 0: if INFO = 1, an singular value did not converge
385: // *
386: // * Further Details
387: // * ===============
388: // *
389: // * Based on contributions by
390: // * Ming Gu and Huan Ren, Computer Science Division, University of
391: // * California at Berkeley, USA
392: // *
393: // * =====================================================================
394: // *
395: // * .. Parameters ..
396: // *
397: // * ..
398: // * .. Local Scalars ..
399: // * ..
400: // * .. External Subroutines ..
401: // * ..
402: // * .. Intrinsic Functions ..
403: // INTRINSIC ABS, MAX;
404: // * ..
405: // * .. Executable Statements ..
406: // *
407: // * Test the input parameters.
408: // *
409:
410: #endregion
411:
412:
413: #region Body
414:
415: INFO = 0;
416: // *
417: if (NL < 1)
418: {
419: INFO = - 1;
420: }
421: else
422: {
423: if (NR < 1)
424: {
425: INFO = - 2;
426: }
427: else
428: {
429: if ((SQRE < 0) || (SQRE > 1))
430: {
431: INFO = - 3;
432: }
433: }
434: }
435: if (INFO != 0)
436: {
437: this._xerbla.Run("DLASD1", - INFO);
438: return;
439: }
440: // *
441: N = NL + NR + 1;
442: M = N + SQRE;
443: // *
444: // * The following values are for bookkeeping purposes only. They are
445: // * integer pointers which indicate the portion of the workspace
446: // * used by a particular array in DLASD2 and DLASD3.
447: // *
448: LDU2 = N;
449: LDVT2 = M;
450: // *
451: IZ = 1;
452: ISIGMA = IZ + M;
453: IU2 = ISIGMA + N;
454: IVT2 = IU2 + LDU2 * N;
455: IQ = IVT2 + LDVT2 * M;
456: // *
457: IDX = 1;
458: IDXC = IDX + N;
459: COLTYP = IDXC + N;
460: IDXP = COLTYP + N;
461: // *
462: // * Scale.
463: // *
464: ORGNRM = Math.Max(Math.Abs(ALPHA), Math.Abs(BETA));
465: D[NL + 1 + o_d] = ZERO;
466: for (I = 1; I <= N; I++)
467: {
468: if (Math.Abs(D[I + o_d]) > ORGNRM)
469: {
470: ORGNRM = Math.Abs(D[I + o_d]);
471: }
472: }
473: this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N
474: , 1, ref D, offset_d, N, ref INFO);
475: ALPHA = ALPHA / ORGNRM;
476: BETA = BETA / ORGNRM;
477: // *
478: // * Deflate singular values.
479: // *
480: this._dlasd2.Run(NL, NR, SQRE, ref K, ref D, offset_d, ref WORK, IZ + o_work
481: , ALPHA, BETA, ref U, offset_u, LDU, ref VT, offset_vt, LDVT
482: , ref WORK, ISIGMA + o_work, ref WORK, IU2 + o_work, LDU2, ref WORK, IVT2 + o_work, LDVT2, ref IWORK, IDXP + o_iwork
483: , ref IWORK, IDX + o_iwork, ref IWORK, IDXC + o_iwork, ref IDXQ, offset_idxq, ref IWORK, COLTYP + o_iwork, ref INFO);
484: // *
485: // * Solve Secular Equation and update singular vectors.
486: // *
487: LDQ = K;
488: this._dlasd3.Run(NL, NR, SQRE, K, ref D, offset_d, ref WORK, IQ + o_work
489: , LDQ, ref WORK, ISIGMA + o_work, ref U, offset_u, LDU, WORK, IU2 + o_work, LDU2
490: , ref VT, offset_vt, LDVT, ref WORK, IVT2 + o_work, LDVT2, IWORK, IDXC + o_iwork, IWORK, COLTYP + o_iwork
491: , ref WORK, IZ + o_work, ref INFO);
492: if (INFO != 0)
493: {
494: return;
495: }
496: // *
497: // * Unscale.
498: // *
499: this._dlascl.Run("G", 0, 0, ONE, ORGNRM, N
500: , 1, ref D, offset_d, N, ref INFO);
501: // *
502: // * Prepare the IDXQ sorting permutation.
503: // *
504: N1 = K;
505: N2 = N - K;
506: this._dlamrg.Run(N1, N2, D, offset_d, 1, - 1, ref IDXQ, offset_idxq);
507: // *
508: return;
509: // *
510: // * End of DLASD1
511: // *
512:
513: #endregion
514:
515: }
516: }
517: }