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: /// DLAIC1 applies one step of incremental condition estimation in
27: /// its simplest version:
28: ///
29: /// Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
30: /// lower triangular matrix L, such that
31: /// twonorm(L*x) = sest
32: /// Then DLAIC1 computes sestpr, s, c such that
33: /// the vector
34: /// [ s*x ]
35: /// xhat = [ c ]
36: /// is an approximate singular vector of
37: /// [ L 0 ]
38: /// Lhat = [ w' gamma ]
39: /// in the sense that
40: /// twonorm(Lhat*xhat) = sestpr.
41: ///
42: /// Depending on JOB, an estimate for the largest or smallest singular
43: /// value is computed.
44: ///
45: /// Note that [s c]' and sestpr**2 is an eigenpair of the system
46: ///
47: /// diag(sest*sest, 0) + [alpha gamma] * [ alpha ]
48: /// [ gamma ]
49: ///
50: /// where alpha = x'*w.
51: ///
52: ///</summary>
53: public class DLAIC1
54: {
55:
56:
57: #region Dependencies
58:
59: DDOT _ddot; DLAMCH _dlamch;
60:
61: #endregion
62:
63:
64: #region Fields
65:
66: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; const double HALF = 0.5E0;
67: const double FOUR = 4.0E0;double ABSALP = 0; double ABSEST = 0; double ABSGAM = 0; double ALPHA = 0; double B = 0;
68: double COSINE = 0;double EPS = 0; double NORMA = 0; double S1 = 0; double S2 = 0; double SINE = 0; double T = 0;
69: double TEST = 0;double TMP = 0; double ZETA1 = 0; double ZETA2 = 0;
70:
71: #endregion
72:
73: public DLAIC1(DDOT ddot, DLAMCH dlamch)
74: {
75:
76:
77: #region Set Dependencies
78:
79: this._ddot = ddot; this._dlamch = dlamch;
80:
81: #endregion
82:
83: }
84:
85: public DLAIC1()
86: {
87:
88:
89: #region Dependencies (Initialization)
90:
91: DDOT ddot = new DDOT();
92: LSAME lsame = new LSAME();
93: DLAMC3 dlamc3 = new DLAMC3();
94: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
95: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
96: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
97: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
98: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
99:
100: #endregion
101:
102:
103: #region Set Dependencies
104:
105: this._ddot = ddot; this._dlamch = dlamch;
106:
107: #endregion
108:
109: }
110: /// <summary>
111: /// Purpose
112: /// =======
113: ///
114: /// DLAIC1 applies one step of incremental condition estimation in
115: /// its simplest version:
116: ///
117: /// Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
118: /// lower triangular matrix L, such that
119: /// twonorm(L*x) = sest
120: /// Then DLAIC1 computes sestpr, s, c such that
121: /// the vector
122: /// [ s*x ]
123: /// xhat = [ c ]
124: /// is an approximate singular vector of
125: /// [ L 0 ]
126: /// Lhat = [ w' gamma ]
127: /// in the sense that
128: /// twonorm(Lhat*xhat) = sestpr.
129: ///
130: /// Depending on JOB, an estimate for the largest or smallest singular
131: /// value is computed.
132: ///
133: /// Note that [s c]' and sestpr**2 is an eigenpair of the system
134: ///
135: /// diag(sest*sest, 0) + [alpha gamma] * [ alpha ]
136: /// [ gamma ]
137: ///
138: /// where alpha = x'*w.
139: ///
140: ///</summary>
141: /// <param name="JOB">
142: /// (input) INTEGER
143: /// = 1: an estimate for the largest singular value is computed.
144: /// = 2: an estimate for the smallest singular value is computed.
145: ///</param>
146: /// <param name="J">
147: /// (input) INTEGER
148: /// Length of X and W
149: ///</param>
150: /// <param name="X">
151: /// (input) DOUBLE PRECISION array, dimension (J)
152: /// The j-vector x.
153: ///</param>
154: /// <param name="SEST">
155: /// (input) DOUBLE PRECISION
156: /// Estimated singular value of j by j matrix L
157: ///</param>
158: /// <param name="W">
159: /// (input) DOUBLE PRECISION array, dimension (J)
160: /// The j-vector w.
161: ///</param>
162: /// <param name="GAMMA">
163: /// (input) DOUBLE PRECISION
164: /// The diagonal element gamma.
165: ///</param>
166: /// <param name="SESTPR">
167: /// (output) DOUBLE PRECISION
168: /// Estimated singular value of (j+1) by (j+1) matrix Lhat.
169: ///</param>
170: /// <param name="S">
171: /// (output) DOUBLE PRECISION
172: /// Sine needed in forming xhat.
173: ///</param>
174: /// <param name="C">
175: /// (output) DOUBLE PRECISION
176: /// Cosine needed in forming xhat.
177: ///</param>
178: public void Run(int JOB, int J, double[] X, int offset_x, double SEST, double[] W, int offset_w, double GAMMA
179: , ref double SESTPR, ref double S, ref double C)
180: {
181:
182: #region Array Index Correction
183:
184: int o_x = -1 + offset_x; int o_w = -1 + offset_w;
185:
186: #endregion
187:
188:
189: #region Prolog
190:
191: // *
192: // * -- LAPACK auxiliary routine (version 3.1) --
193: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
194: // * November 2006
195: // *
196: // * .. Scalar Arguments ..
197: // * ..
198: // * .. Array Arguments ..
199: // * ..
200: // *
201: // * Purpose
202: // * =======
203: // *
204: // * DLAIC1 applies one step of incremental condition estimation in
205: // * its simplest version:
206: // *
207: // * Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
208: // * lower triangular matrix L, such that
209: // * twonorm(L*x) = sest
210: // * Then DLAIC1 computes sestpr, s, c such that
211: // * the vector
212: // * [ s*x ]
213: // * xhat = [ c ]
214: // * is an approximate singular vector of
215: // * [ L 0 ]
216: // * Lhat = [ w' gamma ]
217: // * in the sense that
218: // * twonorm(Lhat*xhat) = sestpr.
219: // *
220: // * Depending on JOB, an estimate for the largest or smallest singular
221: // * value is computed.
222: // *
223: // * Note that [s c]' and sestpr**2 is an eigenpair of the system
224: // *
225: // * diag(sest*sest, 0) + [alpha gamma] * [ alpha ]
226: // * [ gamma ]
227: // *
228: // * where alpha = x'*w.
229: // *
230: // * Arguments
231: // * =========
232: // *
233: // * JOB (input) INTEGER
234: // * = 1: an estimate for the largest singular value is computed.
235: // * = 2: an estimate for the smallest singular value is computed.
236: // *
237: // * J (input) INTEGER
238: // * Length of X and W
239: // *
240: // * X (input) DOUBLE PRECISION array, dimension (J)
241: // * The j-vector x.
242: // *
243: // * SEST (input) DOUBLE PRECISION
244: // * Estimated singular value of j by j matrix L
245: // *
246: // * W (input) DOUBLE PRECISION array, dimension (J)
247: // * The j-vector w.
248: // *
249: // * GAMMA (input) DOUBLE PRECISION
250: // * The diagonal element gamma.
251: // *
252: // * SESTPR (output) DOUBLE PRECISION
253: // * Estimated singular value of (j+1) by (j+1) matrix Lhat.
254: // *
255: // * S (output) DOUBLE PRECISION
256: // * Sine needed in forming xhat.
257: // *
258: // * C (output) DOUBLE PRECISION
259: // * Cosine needed in forming xhat.
260: // *
261: // * =====================================================================
262: // *
263: // * .. Parameters ..
264: // * ..
265: // * .. Local Scalars ..
266: // * ..
267: // * .. Intrinsic Functions ..
268: // INTRINSIC ABS, MAX, SIGN, SQRT;
269: // * ..
270: // * .. External Functions ..
271: // * ..
272: // * .. Executable Statements ..
273: // *
274:
275: #endregion
276:
277:
278: #region Body
279:
280: EPS = this._dlamch.Run("Epsilon");
281: ALPHA = this._ddot.Run(J, X, offset_x, 1, W, offset_w, 1);
282: // *
283: ABSALP = Math.Abs(ALPHA);
284: ABSGAM = Math.Abs(GAMMA);
285: ABSEST = Math.Abs(SEST);
286: // *
287: if (JOB == 1)
288: {
289: // *
290: // * Estimating largest singular value
291: // *
292: // * special cases
293: // *
294: if (SEST == ZERO)
295: {
296: S1 = Math.Max(ABSGAM, ABSALP);
297: if (S1 == ZERO)
298: {
299: S = ZERO;
300: C = ONE;
301: SESTPR = ZERO;
302: }
303: else
304: {
305: S = ALPHA / S1;
306: C = GAMMA / S1;
307: TMP = Math.Sqrt(S * S + C * C);
308: S = S / TMP;
309: C = C / TMP;
310: SESTPR = S1 * TMP;
311: }
312: return;
313: }
314: else
315: {
316: if (ABSGAM <= EPS * ABSEST)
317: {
318: S = ONE;
319: C = ZERO;
320: TMP = Math.Max(ABSEST, ABSALP);
321: S1 = ABSEST / TMP;
322: S2 = ABSALP / TMP;
323: SESTPR = TMP * Math.Sqrt(S1 * S1 + S2 * S2);
324: return;
325: }
326: else
327: {
328: if (ABSALP <= EPS * ABSEST)
329: {
330: S1 = ABSGAM;
331: S2 = ABSEST;
332: if (S1 <= S2)
333: {
334: S = ONE;
335: C = ZERO;
336: SESTPR = S2;
337: }
338: else
339: {
340: S = ZERO;
341: C = ONE;
342: SESTPR = S1;
343: }
344: return;
345: }
346: else
347: {
348: if (ABSEST <= EPS * ABSALP || ABSEST <= EPS * ABSGAM)
349: {
350: S1 = ABSGAM;
351: S2 = ABSALP;
352: if (S1 <= S2)
353: {
354: TMP = S1 / S2;
355: S = Math.Sqrt(ONE + TMP * TMP);
356: SESTPR = S2 * S;
357: C = (GAMMA / S2) / S;
358: S = FortranLib.Sign(ONE,ALPHA) / S;
359: }
360: else
361: {
362: TMP = S2 / S1;
363: C = Math.Sqrt(ONE + TMP * TMP);
364: SESTPR = S1 * C;
365: S = (ALPHA / S1) / C;
366: C = FortranLib.Sign(ONE,GAMMA) / C;
367: }
368: return;
369: }
370: else
371: {
372: // *
373: // * normal case
374: // *
375: ZETA1 = ALPHA / ABSEST;
376: ZETA2 = GAMMA / ABSEST;
377: // *
378: B = (ONE - ZETA1 * ZETA1 - ZETA2 * ZETA2) * HALF;
379: C = ZETA1 * ZETA1;
380: if (B > ZERO)
381: {
382: T = C / (B + Math.Sqrt(B * B + C));
383: }
384: else
385: {
386: T = Math.Sqrt(B * B + C) - B;
387: }
388: // *
389: SINE = - ZETA1 / T;
390: COSINE = - ZETA2 / (ONE + T);
391: TMP = Math.Sqrt(SINE * SINE + COSINE * COSINE);
392: S = SINE / TMP;
393: C = COSINE / TMP;
394: SESTPR = Math.Sqrt(T + ONE) * ABSEST;
395: return;
396: }
397: }
398: }
399: }
400: // *
401: }
402: else
403: {
404: if (JOB == 2)
405: {
406: // *
407: // * Estimating smallest singular value
408: // *
409: // * special cases
410: // *
411: if (SEST == ZERO)
412: {
413: SESTPR = ZERO;
414: if (Math.Max(ABSGAM, ABSALP) == ZERO)
415: {
416: SINE = ONE;
417: COSINE = ZERO;
418: }
419: else
420: {
421: SINE = - GAMMA;
422: COSINE = ALPHA;
423: }
424: S1 = Math.Max(Math.Abs(SINE), Math.Abs(COSINE));
425: S = SINE / S1;
426: C = COSINE / S1;
427: TMP = Math.Sqrt(S * S + C * C);
428: S = S / TMP;
429: C = C / TMP;
430: return;
431: }
432: else
433: {
434: if (ABSGAM <= EPS * ABSEST)
435: {
436: S = ZERO;
437: C = ONE;
438: SESTPR = ABSGAM;
439: return;
440: }
441: else
442: {
443: if (ABSALP <= EPS * ABSEST)
444: {
445: S1 = ABSGAM;
446: S2 = ABSEST;
447: if (S1 <= S2)
448: {
449: S = ZERO;
450: C = ONE;
451: SESTPR = S1;
452: }
453: else
454: {
455: S = ONE;
456: C = ZERO;
457: SESTPR = S2;
458: }
459: return;
460: }
461: else
462: {
463: if (ABSEST <= EPS * ABSALP || ABSEST <= EPS * ABSGAM)
464: {
465: S1 = ABSGAM;
466: S2 = ABSALP;
467: if (S1 <= S2)
468: {
469: TMP = S1 / S2;
470: C = Math.Sqrt(ONE + TMP * TMP);
471: SESTPR = ABSEST * (TMP / C);
472: S = - (GAMMA / S2) / C;
473: C = FortranLib.Sign(ONE,ALPHA) / C;
474: }
475: else
476: {
477: TMP = S2 / S1;
478: S = Math.Sqrt(ONE + TMP * TMP);
479: SESTPR = ABSEST / S;
480: C = (ALPHA / S1) / S;
481: S = - FortranLib.Sign(ONE,GAMMA) / S;
482: }
483: return;
484: }
485: else
486: {
487: // *
488: // * normal case
489: // *
490: ZETA1 = ALPHA / ABSEST;
491: ZETA2 = GAMMA / ABSEST;
492: // *
493: NORMA = Math.Max(ONE + ZETA1 * ZETA1 + Math.Abs(ZETA1 * ZETA2), Math.Abs(ZETA1 * ZETA2) + ZETA2 * ZETA2);
494: // *
495: // * See if root is closer to zero or to ONE
496: // *
497: TEST = ONE + TWO * (ZETA1 - ZETA2) * (ZETA1 + ZETA2);
498: if (TEST >= ZERO)
499: {
500: // *
501: // * root is close to zero, compute directly
502: // *
503: B = (ZETA1 * ZETA1 + ZETA2 * ZETA2 + ONE) * HALF;
504: C = ZETA2 * ZETA2;
505: T = C / (B + Math.Sqrt(Math.Abs(B * B - C)));
506: SINE = ZETA1 / (ONE - T);
507: COSINE = - ZETA2 / T;
508: SESTPR = Math.Sqrt(T + FOUR * EPS * EPS * NORMA) * ABSEST;
509: }
510: else
511: {
512: // *
513: // * root is closer to ONE, shift by that amount
514: // *
515: B = (ZETA2 * ZETA2 + ZETA1 * ZETA1 - ONE) * HALF;
516: C = ZETA1 * ZETA1;
517: if (B >= ZERO)
518: {
519: T = - C / (B + Math.Sqrt(B * B + C));
520: }
521: else
522: {
523: T = B - Math.Sqrt(B * B + C);
524: }
525: SINE = - ZETA1 / T;
526: COSINE = - ZETA2 / (ONE + T);
527: SESTPR = Math.Sqrt(ONE + T + FOUR * EPS * EPS * NORMA) * ABSEST;
528: }
529: TMP = Math.Sqrt(SINE * SINE + COSINE * COSINE);
530: S = SINE / TMP;
531: C = COSINE / TMP;
532: return;
533: // *
534: }
535: }
536: }
537: }
538: }
539: }
540: return;
541: // *
542: // * End of DLAIC1
543: // *
544:
545: #endregion
546:
547: }
548: }
549: }