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: /// DLASD8 finds the square roots of the roots of the secular equation,
27: /// as defined by the values in DSIGMA and Z. It makes the appropriate
28: /// calls to DLASD4, and stores, for each element in D, the distance
29: /// to its two nearest poles (elements in DSIGMA). It also updates
30: /// the arrays VF and VL, the first and last components of all the
31: /// right singular vectors of the original bidiagonal matrix.
32: ///
33: /// DLASD8 is called from DLASD6.
34: ///
35: ///</summary>
36: public class DLASD8
37: {
38:
39:
40: #region Dependencies
41:
42: DCOPY _dcopy; DLASCL _dlascl; DLASD4 _dlasd4; DLASET _dlaset; XERBLA _xerbla; DDOT _ddot; DLAMC3 _dlamc3; DNRM2 _dnrm2;
43:
44: #endregion
45:
46:
47: #region Fields
48:
49: const double ONE = 1.0E+0; int I = 0; int IWK1 = 0; int IWK2 = 0; int IWK2I = 0; int IWK3 = 0; int IWK3I = 0; int J = 0;
50: double DIFLJ = 0;double DIFRJ = 0; double DJ = 0; double DSIGJ = 0; double DSIGJP = 0; double RHO = 0; double TEMP = 0;
51:
52: #endregion
53:
54: public DLASD8(DCOPY dcopy, DLASCL dlascl, DLASD4 dlasd4, DLASET dlaset, XERBLA xerbla, DDOT ddot, DLAMC3 dlamc3, DNRM2 dnrm2)
55: {
56:
57:
58: #region Set Dependencies
59:
60: this._dcopy = dcopy; this._dlascl = dlascl; this._dlasd4 = dlasd4; this._dlaset = dlaset; this._xerbla = xerbla;
61: this._ddot = ddot;this._dlamc3 = dlamc3; this._dnrm2 = dnrm2;
62:
63: #endregion
64:
65: }
66:
67: public DLASD8()
68: {
69:
70:
71: #region Dependencies (Initialization)
72:
73: DCOPY dcopy = new DCOPY();
74: LSAME lsame = new LSAME();
75: DLAMC3 dlamc3 = new DLAMC3();
76: XERBLA xerbla = new XERBLA();
77: DLASD5 dlasd5 = new DLASD5();
78: DDOT ddot = new DDOT();
79: DNRM2 dnrm2 = new DNRM2();
80: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
81: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
82: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
83: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
84: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
85: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
86: DLAED6 dlaed6 = new DLAED6(dlamch);
87: DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
88: DLASET dlaset = new DLASET(lsame);
89:
90: #endregion
91:
92:
93: #region Set Dependencies
94:
95: this._dcopy = dcopy; this._dlascl = dlascl; this._dlasd4 = dlasd4; this._dlaset = dlaset; this._xerbla = xerbla;
96: this._ddot = ddot;this._dlamc3 = dlamc3; this._dnrm2 = dnrm2;
97:
98: #endregion
99:
100: }
101: /// <summary>
102: /// Purpose
103: /// =======
104: ///
105: /// DLASD8 finds the square roots of the roots of the secular equation,
106: /// as defined by the values in DSIGMA and Z. It makes the appropriate
107: /// calls to DLASD4, and stores, for each element in D, the distance
108: /// to its two nearest poles (elements in DSIGMA). It also updates
109: /// the arrays VF and VL, the first and last components of all the
110: /// right singular vectors of the original bidiagonal matrix.
111: ///
112: /// DLASD8 is called from DLASD6.
113: ///
114: ///</summary>
115: /// <param name="ICOMPQ">
116: /// (input) INTEGER
117: /// Specifies whether singular vectors are to be computed in
118: /// factored form in the calling routine:
119: /// = 0: Compute singular values only.
120: /// = 1: Compute singular vectors in factored form as well.
121: ///</param>
122: /// <param name="K">
123: /// (input) INTEGER
124: /// The number of terms in the rational function to be solved
125: /// by DLASD4. K .GE. 1.
126: ///</param>
127: /// <param name="D">
128: /// (output) DOUBLE PRECISION array, dimension ( K )
129: /// On output, D contains the updated singular values.
130: ///</param>
131: /// <param name="Z">
132: /// (input) DOUBLE PRECISION array, dimension ( K )
133: /// The first K elements of this array contain the components
134: /// of the deflation-adjusted updating row vector.
135: ///</param>
136: /// <param name="VF">
137: /// (input/output) DOUBLE PRECISION array, dimension ( K )
138: /// On entry, VF contains information passed through DBEDE8.
139: /// On exit, VF contains the first K components of the first
140: /// components of all right singular vectors of the bidiagonal
141: /// matrix.
142: ///</param>
143: /// <param name="VL">
144: /// (input/output) DOUBLE PRECISION array, dimension ( K )
145: /// On entry, VL contains information passed through DBEDE8.
146: /// On exit, VL contains the first K components of the last
147: /// components of all right singular vectors of the bidiagonal
148: /// matrix.
149: ///</param>
150: /// <param name="DIFL">
151: /// (output) DOUBLE PRECISION array, dimension ( K )
152: /// On exit, DIFL(I) = D(I) - DSIGMA(I).
153: ///</param>
154: /// <param name="DIFR">
155: /// (output) DOUBLE PRECISION array,
156: /// dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
157: /// dimension ( K ) if ICOMPQ = 0.
158: /// On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
159: /// defined and will not be referenced.
160: ///
161: /// If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
162: /// normalizing factors for the right singular vector matrix.
163: ///</param>
164: /// <param name="LDDIFR">
165: /// (input) INTEGER
166: /// The leading dimension of DIFR, must be at least K.
167: ///</param>
168: /// <param name="DSIGMA">
169: /// (input) DOUBLE PRECISION array, dimension ( K )
170: /// The first K elements of this array contain the old roots
171: /// of the deflated updating problem. These are the poles
172: /// of the secular equation.
173: ///</param>
174: /// <param name="WORK">
175: /// (workspace) DOUBLE PRECISION array, dimension at least 3 * K
176: ///</param>
177: /// <param name="INFO">
178: /// (output) INTEGER
179: /// = 0: successful exit.
180: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
181: /// .GT. 0: if INFO = 1, an singular value did not converge
182: ///</param>
183: public void Run(int ICOMPQ, int K, ref double[] D, int offset_d, ref double[] Z, int offset_z, ref double[] VF, int offset_vf, ref double[] VL, int offset_vl
184: , ref double[] DIFL, int offset_difl, ref double[] DIFR, int offset_difr, int LDDIFR, ref double[] DSIGMA, int offset_dsigma, ref double[] WORK, int offset_work, ref int INFO)
185: {
186:
187: #region Array Index Correction
188:
189: int o_d = -1 + offset_d; int o_z = -1 + offset_z; int o_vf = -1 + offset_vf; int o_vl = -1 + offset_vl;
190: int o_difl = -1 + offset_difl; int o_difr = -1 - LDDIFR + offset_difr; int o_dsigma = -1 + offset_dsigma;
191: int o_work = -1 + offset_work;
192:
193: #endregion
194:
195:
196: #region Prolog
197:
198: // *
199: // * -- LAPACK auxiliary routine (version 3.1) --
200: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
201: // * November 2006
202: // *
203: // * .. Scalar Arguments ..
204: // * ..
205: // * .. Array Arguments ..
206: // * ..
207: // *
208: // * Purpose
209: // * =======
210: // *
211: // * DLASD8 finds the square roots of the roots of the secular equation,
212: // * as defined by the values in DSIGMA and Z. It makes the appropriate
213: // * calls to DLASD4, and stores, for each element in D, the distance
214: // * to its two nearest poles (elements in DSIGMA). It also updates
215: // * the arrays VF and VL, the first and last components of all the
216: // * right singular vectors of the original bidiagonal matrix.
217: // *
218: // * DLASD8 is called from DLASD6.
219: // *
220: // * Arguments
221: // * =========
222: // *
223: // * ICOMPQ (input) INTEGER
224: // * Specifies whether singular vectors are to be computed in
225: // * factored form in the calling routine:
226: // * = 0: Compute singular values only.
227: // * = 1: Compute singular vectors in factored form as well.
228: // *
229: // * K (input) INTEGER
230: // * The number of terms in the rational function to be solved
231: // * by DLASD4. K >= 1.
232: // *
233: // * D (output) DOUBLE PRECISION array, dimension ( K )
234: // * On output, D contains the updated singular values.
235: // *
236: // * Z (input) DOUBLE PRECISION array, dimension ( K )
237: // * The first K elements of this array contain the components
238: // * of the deflation-adjusted updating row vector.
239: // *
240: // * VF (input/output) DOUBLE PRECISION array, dimension ( K )
241: // * On entry, VF contains information passed through DBEDE8.
242: // * On exit, VF contains the first K components of the first
243: // * components of all right singular vectors of the bidiagonal
244: // * matrix.
245: // *
246: // * VL (input/output) DOUBLE PRECISION array, dimension ( K )
247: // * On entry, VL contains information passed through DBEDE8.
248: // * On exit, VL contains the first K components of the last
249: // * components of all right singular vectors of the bidiagonal
250: // * matrix.
251: // *
252: // * DIFL (output) DOUBLE PRECISION array, dimension ( K )
253: // * On exit, DIFL(I) = D(I) - DSIGMA(I).
254: // *
255: // * DIFR (output) DOUBLE PRECISION array,
256: // * dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
257: // * dimension ( K ) if ICOMPQ = 0.
258: // * On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
259: // * defined and will not be referenced.
260: // *
261: // * If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
262: // * normalizing factors for the right singular vector matrix.
263: // *
264: // * LDDIFR (input) INTEGER
265: // * The leading dimension of DIFR, must be at least K.
266: // *
267: // * DSIGMA (input) DOUBLE PRECISION array, dimension ( K )
268: // * The first K elements of this array contain the old roots
269: // * of the deflated updating problem. These are the poles
270: // * of the secular equation.
271: // *
272: // * WORK (workspace) DOUBLE PRECISION array, dimension at least 3 * K
273: // *
274: // * INFO (output) INTEGER
275: // * = 0: successful exit.
276: // * < 0: if INFO = -i, the i-th argument had an illegal value.
277: // * > 0: if INFO = 1, an singular value did not converge
278: // *
279: // * Further Details
280: // * ===============
281: // *
282: // * Based on contributions by
283: // * Ming Gu and Huan Ren, Computer Science Division, University of
284: // * California at Berkeley, USA
285: // *
286: // * =====================================================================
287: // *
288: // * .. Parameters ..
289: // * ..
290: // * .. Local Scalars ..
291: // * ..
292: // * .. External Subroutines ..
293: // * ..
294: // * .. External Functions ..
295: // * ..
296: // * .. Intrinsic Functions ..
297: // INTRINSIC ABS, SIGN, SQRT;
298: // * ..
299: // * .. Executable Statements ..
300: // *
301: // * Test the input parameters.
302: // *
303:
304: #endregion
305:
306:
307: #region Body
308:
309: INFO = 0;
310: // *
311: if ((ICOMPQ < 0) || (ICOMPQ > 1))
312: {
313: INFO = - 1;
314: }
315: else
316: {
317: if (K < 1)
318: {
319: INFO = - 2;
320: }
321: else
322: {
323: if (LDDIFR < K)
324: {
325: INFO = - 9;
326: }
327: }
328: }
329: if (INFO != 0)
330: {
331: this._xerbla.Run("DLASD8", - INFO);
332: return;
333: }
334: // *
335: // * Quick return if possible
336: // *
337: if (K == 1)
338: {
339: D[1 + o_d] = Math.Abs(Z[1 + o_z]);
340: DIFL[1 + o_difl] = D[1 + o_d];
341: if (ICOMPQ == 1)
342: {
343: DIFL[2 + o_difl] = ONE;
344: DIFR[1+2 * LDDIFR + o_difr] = ONE;
345: }
346: return;
347: }
348: // *
349: // * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
350: // * be computed with high relative accuracy (barring over/underflow).
351: // * This is a problem on machines without a guard digit in
352: // * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
353: // * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
354: // * which on any of these machines zeros out the bottommost
355: // * bit of DSIGMA(I) if it is 1; this makes the subsequent
356: // * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
357: // * occurs. On binary machines with a guard digit (almost all
358: // * machines) it does not change DSIGMA(I) at all. On hexadecimal
359: // * and decimal machines with a guard digit, it slightly
360: // * changes the bottommost bits of DSIGMA(I). It does not account
361: // * for hexadecimal or decimal machines without guard digits
362: // * (we know of none). We use a subroutine call to compute
363: // * 2*DSIGMA(I) to prevent optimizing compilers from eliminating
364: // * this code.
365: // *
366: for (I = 1; I <= K; I++)
367: {
368: DSIGMA[I + o_dsigma] = this._dlamc3.Run(DSIGMA[I + o_dsigma], DSIGMA[I + o_dsigma]) - DSIGMA[I + o_dsigma];
369: }
370: // *
371: // * Book keeping.
372: // *
373: IWK1 = 1;
374: IWK2 = IWK1 + K;
375: IWK3 = IWK2 + K;
376: IWK2I = IWK2 - 1;
377: IWK3I = IWK3 - 1;
378: // *
379: // * Normalize Z.
380: // *
381: RHO = this._dnrm2.Run(K, Z, offset_z, 1);
382: this._dlascl.Run("G", 0, 0, RHO, ONE, K
383: , 1, ref Z, offset_z, K, ref INFO);
384: RHO = RHO * RHO;
385: // *
386: // * Initialize WORK(IWK3).
387: // *
388: this._dlaset.Run("A", K, 1, ONE, ONE, ref WORK, IWK3 + o_work
389: , K);
390: // *
391: // * Compute the updated singular values, the arrays DIFL, DIFR,
392: // * and the updated Z.
393: // *
394: for (J = 1; J <= K; J++)
395: {
396: this._dlasd4.Run(K, J, DSIGMA, offset_dsigma, Z, offset_z, ref WORK, IWK1 + o_work, RHO
397: , ref D[J + o_d], ref WORK, IWK2 + o_work, ref INFO);
398: // *
399: // * If the root finder fails, the computation is terminated.
400: // *
401: if (INFO != 0)
402: {
403: return;
404: }
405: WORK[IWK3I + J + o_work] = WORK[IWK3I + J + o_work] * WORK[J + o_work] * WORK[IWK2I + J + o_work];
406: DIFL[J + o_difl] = - WORK[J + o_work];
407: DIFR[J+1 * LDDIFR + o_difr] = - WORK[J + 1 + o_work];
408: for (I = 1; I <= J - 1; I++)
409: {
410: WORK[IWK3I + I + o_work] = WORK[IWK3I + I + o_work] * WORK[I + o_work] * WORK[IWK2I + I + o_work] / (DSIGMA[I + o_dsigma] - DSIGMA[J + o_dsigma]) / (DSIGMA[I + o_dsigma] + DSIGMA[J + o_dsigma]);
411: }
412: for (I = J + 1; I <= K; I++)
413: {
414: WORK[IWK3I + I + o_work] = WORK[IWK3I + I + o_work] * WORK[I + o_work] * WORK[IWK2I + I + o_work] / (DSIGMA[I + o_dsigma] - DSIGMA[J + o_dsigma]) / (DSIGMA[I + o_dsigma] + DSIGMA[J + o_dsigma]);
415: }
416: }
417: // *
418: // * Compute updated Z.
419: // *
420: for (I = 1; I <= K; I++)
421: {
422: Z[I + o_z] = FortranLib.Sign(Math.Sqrt(Math.Abs(WORK[IWK3I + I + o_work])),Z[I + o_z]);
423: }
424: // *
425: // * Update VF and VL.
426: // *
427: for (J = 1; J <= K; J++)
428: {
429: DIFLJ = DIFL[J + o_difl];
430: DJ = D[J + o_d];
431: DSIGJ = - DSIGMA[J + o_dsigma];
432: if (J < K)
433: {
434: DIFRJ = - DIFR[J+1 * LDDIFR + o_difr];
435: DSIGJP = - DSIGMA[J + 1 + o_dsigma];
436: }
437: WORK[J + o_work] = - Z[J + o_z] / DIFLJ / (DSIGMA[J + o_dsigma] + DJ);
438: for (I = 1; I <= J - 1; I++)
439: {
440: WORK[I + o_work] = Z[I + o_z] / (this._dlamc3.Run(DSIGMA[I + o_dsigma], DSIGJ) - DIFLJ) / (DSIGMA[I + o_dsigma] + DJ);
441: }
442: for (I = J + 1; I <= K; I++)
443: {
444: WORK[I + o_work] = Z[I + o_z] / (this._dlamc3.Run(DSIGMA[I + o_dsigma], DSIGJP) + DIFRJ) / (DSIGMA[I + o_dsigma] + DJ);
445: }
446: TEMP = this._dnrm2.Run(K, WORK, offset_work, 1);
447: WORK[IWK2I + J + o_work] = this._ddot.Run(K, WORK, offset_work, 1, VF, offset_vf, 1) / TEMP;
448: WORK[IWK3I + J + o_work] = this._ddot.Run(K, WORK, offset_work, 1, VL, offset_vl, 1) / TEMP;
449: if (ICOMPQ == 1)
450: {
451: DIFR[J+2 * LDDIFR + o_difr] = TEMP;
452: }
453: }
454: // *
455: this._dcopy.Run(K, WORK, IWK2 + o_work, 1, ref VF, offset_vf, 1);
456: this._dcopy.Run(K, WORK, IWK3 + o_work, 1, ref VL, offset_vl, 1);
457: // *
458: return;
459: // *
460: // * End of DLASD8
461: // *
462:
463: #endregion
464:
465: }
466: }
467: }