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:  }