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:      /// DLAED3 finds the roots of the secular equation, as defined by the
27:      /// values in D, W, and RHO, between 1 and K.  It makes the
28:      /// appropriate calls to DLAED4 and then updates the eigenvectors by
29:      /// multiplying the matrix of eigenvectors of the pair of eigensystems
30:      /// being combined by the matrix of eigenvectors of the K-by-K system
31:      /// which is solved here.
32:      ///
33:      /// This code makes very mild assumptions about floating point
34:      /// arithmetic. It will work on machines with a guard digit in
35:      /// add/subtract, or on those binary machines without guard digits
36:      /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
37:      /// It could conceivably fail on hexadecimal or decimal machines
38:      /// without guard digits, but we know of none.
39:      ///
40:      ///</summary>
41:      public class DLAED3
42:      {
43:
44:
45:          #region Dependencies
46:
47:          DLAMC3 _dlamc3; DNRM2 _dnrm2; DCOPY _dcopy; DGEMM _dgemm; DLACPY _dlacpy; DLAED4 _dlaed4; DLASET _dlaset; XERBLA _xerbla;
48:
49:          #endregion
50:
51:
52:          #region Fields
53:
54:          const double ONE = 1.0E0; const double ZERO = 0.0E0; int I = 0; int II = 0; int IQ2 = 0; int J = 0; int N12 = 0;
55:          int N2 = 0;int N23 = 0; double TEMP = 0;
56:
57:          #endregion
58:
59:          public DLAED3(DLAMC3 dlamc3, DNRM2 dnrm2, DCOPY dcopy, DGEMM dgemm, DLACPY dlacpy, DLAED4 dlaed4, DLASET dlaset, XERBLA xerbla)
60:          {
61:
62:
63:              #region Set Dependencies
64:
65:              this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy;
66:              this._dlaed4 = dlaed4;this._dlaset = dlaset; this._xerbla = xerbla;
67:
68:              #endregion
69:
70:          }
71:
72:          public DLAED3()
73:          {
74:
75:
76:              #region Dependencies (Initialization)
77:
78:              DLAMC3 dlamc3 = new DLAMC3();
79:              DNRM2 dnrm2 = new DNRM2();
80:              DCOPY dcopy = new DCOPY();
81:              LSAME lsame = new LSAME();
82:              XERBLA xerbla = new XERBLA();
83:              DLAED5 dlaed5 = new DLAED5();
84:              DGEMM dgemm = new DGEMM(lsame, xerbla);
85:              DLACPY dlacpy = new DLACPY(lsame);
86:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
87:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
88:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
89:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
90:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
91:              DLAED6 dlaed6 = new DLAED6(dlamch);
92:              DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
93:              DLASET dlaset = new DLASET(lsame);
94:
95:              #endregion
96:
97:
98:              #region Set Dependencies
99:
100:              this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy;
101:              this._dlaed4 = dlaed4;this._dlaset = dlaset; this._xerbla = xerbla;
102:
103:              #endregion
104:
105:          }
106:          /// <summary>
107:          /// Purpose
108:          /// =======
109:          ///
110:          /// DLAED3 finds the roots of the secular equation, as defined by the
111:          /// values in D, W, and RHO, between 1 and K.  It makes the
112:          /// appropriate calls to DLAED4 and then updates the eigenvectors by
113:          /// multiplying the matrix of eigenvectors of the pair of eigensystems
114:          /// being combined by the matrix of eigenvectors of the K-by-K system
115:          /// which is solved here.
116:          ///
117:          /// This code makes very mild assumptions about floating point
118:          /// arithmetic. It will work on machines with a guard digit in
119:          /// add/subtract, or on those binary machines without guard digits
120:          /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
121:          /// It could conceivably fail on hexadecimal or decimal machines
122:          /// without guard digits, but we know of none.
123:          ///
124:          ///</summary>
125:          /// <param name="K">
126:          /// (input) INTEGER
127:          /// The number of terms in the rational function to be solved by
128:          /// DLAED4.  K .GE. 0.
129:          ///</param>
130:          /// <param name="N">
131:          /// (input) INTEGER
132:          /// The number of rows and columns in the Q matrix.
133:          /// N .GE. K (deflation may result in N.GT.K).
134:          ///</param>
135:          /// <param name="N1">
136:          /// (input) INTEGER
137:          /// The location of the last eigenvalue in the leading submatrix.
138:          /// min(1,N) .LE. N1 .LE. N/2.
139:          ///</param>
140:          /// <param name="D">
141:          /// (output) DOUBLE PRECISION array, dimension (N)
142:          /// D(I) contains the updated eigenvalues for
143:          /// 1 .LE. I .LE. K.
144:          ///</param>
145:          /// <param name="Q">
146:          /// (output) DOUBLE PRECISION array, dimension (LDQ,N)
147:          /// Initially the first K columns are used as workspace.
148:          /// On output the columns 1 to K contain
149:          /// the updated eigenvectors.
150:          ///</param>
151:          /// <param name="LDQ">
152:          /// (input) INTEGER
153:          /// The leading dimension of the array Q.  LDQ .GE. max(1,N).
154:          ///</param>
155:          /// <param name="RHO">
156:          /// (input) DOUBLE PRECISION
157:          /// The value of the parameter in the rank one update equation.
158:          /// RHO .GE. 0 required.
159:          ///</param>
160:          /// <param name="DLAMDA">
161:          /// (input/output) DOUBLE PRECISION array, dimension (K)
162:          /// The first K elements of this array contain the old roots
163:          /// of the deflated updating problem.  These are the poles
164:          /// of the secular equation. May be changed on output by
165:          /// having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
166:          /// Cray-2, or Cray C-90, as described above.
167:          ///</param>
168:          /// <param name="Q2">
169:          /// (input) DOUBLE PRECISION array, dimension (LDQ2, N)
170:          /// The first K columns of this matrix contain the non-deflated
171:          /// eigenvectors for the split problem.
172:          ///</param>
173:          /// <param name="INDX">
174:          /// (input) INTEGER array, dimension (N)
175:          /// The permutation used to arrange the columns of the deflated
176:          /// Q matrix into three groups (see DLAED2).
177:          /// The rows of the eigenvectors found by DLAED4 must be likewise
178:          /// permuted before the matrix multiply can take place.
179:          ///</param>
180:          /// <param name="CTOT">
181:          /// (input) INTEGER array, dimension (4)
182:          /// A count of the total number of the various types of columns
183:          /// in Q, as described in INDX.  The fourth column type is any
184:          /// column which has been deflated.
185:          ///</param>
186:          /// <param name="W">
187:          /// (input/output) DOUBLE PRECISION array, dimension (K)
188:          /// The first K elements of this array contain the components
189:          /// of the deflation-adjusted updating vector. Destroyed on
190:          /// output.
191:          ///</param>
192:          /// <param name="S">
193:          /// (workspace) DOUBLE PRECISION array, dimension (N1 + 1)*K
194:          /// Will contain the eigenvectors of the repaired matrix which
195:          /// will be multiplied by the previously accumulated eigenvectors
196:          /// to update the system.
197:          ///</param>
198:          /// <param name="INFO">
199:          /// (output) INTEGER
200:          /// = 0:  successful exit.
201:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
202:          /// .GT. 0:  if INFO = 1, an eigenvalue did not converge
203:          ///</param>
204:          public void Run(int K, int N, int N1, ref double[] D, int offset_d, ref double[] Q, int offset_q, int LDQ
205:                           , double RHO, ref double[] DLAMDA, int offset_dlamda, double[] Q2, int offset_q2, int[] INDX, int offset_indx, int[] CTOT, int offset_ctot, ref double[] W, int offset_w
206:                           , ref double[] S, int offset_s, ref int INFO)
207:          {
208:
209:              #region Array Index Correction
210:
211:               int o_d = -1 + offset_d;  int o_q = -1 - LDQ + offset_q;  int o_dlamda = -1 + offset_dlamda;
212:               int o_q2 = -1 + offset_q2; int o_indx = -1 + offset_indx;  int o_ctot = -1 + offset_ctot;  int o_w = -1 + offset_w;
213:               int o_s = -1 + offset_s;
214:
215:              #endregion
216:
217:
218:              #region Prolog
219:
220:              // *
221:              // *  -- LAPACK routine (version 3.1) --
222:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
223:              // *     November 2006
224:              // *
225:              // *     .. Scalar Arguments ..
226:              // *     ..
227:              // *     .. Array Arguments ..
228:              // *     ..
229:              // *
230:              // *  Purpose
231:              // *  =======
232:              // *
233:              // *  DLAED3 finds the roots of the secular equation, as defined by the
234:              // *  values in D, W, and RHO, between 1 and K.  It makes the
235:              // *  appropriate calls to DLAED4 and then updates the eigenvectors by
236:              // *  multiplying the matrix of eigenvectors of the pair of eigensystems
237:              // *  being combined by the matrix of eigenvectors of the K-by-K system
238:              // *  which is solved here.
239:              // *
240:              // *  This code makes very mild assumptions about floating point
241:              // *  arithmetic. It will work on machines with a guard digit in
242:              // *  add/subtract, or on those binary machines without guard digits
243:              // *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
244:              // *  It could conceivably fail on hexadecimal or decimal machines
245:              // *  without guard digits, but we know of none.
246:              // *
247:              // *  Arguments
248:              // *  =========
249:              // *
250:              // *  K       (input) INTEGER
251:              // *          The number of terms in the rational function to be solved by
252:              // *          DLAED4.  K >= 0.
253:              // *
254:              // *  N       (input) INTEGER
255:              // *          The number of rows and columns in the Q matrix.
256:              // *          N >= K (deflation may result in N>K).
257:              // *
258:              // *  N1      (input) INTEGER
259:              // *          The location of the last eigenvalue in the leading submatrix.
260:              // *          min(1,N) <= N1 <= N/2.
261:              // *
262:              // *  D       (output) DOUBLE PRECISION array, dimension (N)
263:              // *          D(I) contains the updated eigenvalues for
264:              // *          1 <= I <= K.
265:              // *
266:              // *  Q       (output) DOUBLE PRECISION array, dimension (LDQ,N)
267:              // *          Initially the first K columns are used as workspace.
268:              // *          On output the columns 1 to K contain
269:              // *          the updated eigenvectors.
270:              // *
271:              // *  LDQ     (input) INTEGER
272:              // *          The leading dimension of the array Q.  LDQ >= max(1,N).
273:              // *
274:              // *  RHO     (input) DOUBLE PRECISION
275:              // *          The value of the parameter in the rank one update equation.
276:              // *          RHO >= 0 required.
277:              // *
278:              // *  DLAMDA  (input/output) DOUBLE PRECISION array, dimension (K)
279:              // *          The first K elements of this array contain the old roots
280:              // *          of the deflated updating problem.  These are the poles
281:              // *          of the secular equation. May be changed on output by
282:              // *          having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
283:              // *          Cray-2, or Cray C-90, as described above.
284:              // *
285:              // *  Q2      (input) DOUBLE PRECISION array, dimension (LDQ2, N)
286:              // *          The first K columns of this matrix contain the non-deflated
287:              // *          eigenvectors for the split problem.
288:              // *
289:              // *  INDX    (input) INTEGER array, dimension (N)
290:              // *          The permutation used to arrange the columns of the deflated
291:              // *          Q matrix into three groups (see DLAED2).
292:              // *          The rows of the eigenvectors found by DLAED4 must be likewise
293:              // *          permuted before the matrix multiply can take place.
294:              // *
295:              // *  CTOT    (input) INTEGER array, dimension (4)
296:              // *          A count of the total number of the various types of columns
297:              // *          in Q, as described in INDX.  The fourth column type is any
298:              // *          column which has been deflated.
299:              // *
300:              // *  W       (input/output) DOUBLE PRECISION array, dimension (K)
301:              // *          The first K elements of this array contain the components
302:              // *          of the deflation-adjusted updating vector. Destroyed on
303:              // *          output.
304:              // *
305:              // *  S       (workspace) DOUBLE PRECISION array, dimension (N1 + 1)*K
306:              // *          Will contain the eigenvectors of the repaired matrix which
307:              // *          will be multiplied by the previously accumulated eigenvectors
308:              // *          to update the system.
309:              // *
310:              // *  LDS     (input) INTEGER
311:              // *          The leading dimension of S.  LDS >= max(1,K).
312:              // *
313:              // *  INFO    (output) INTEGER
314:              // *          = 0:  successful exit.
315:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
316:              // *          > 0:  if INFO = 1, an eigenvalue did not converge
317:              // *
318:              // *  Further Details
319:              // *  ===============
320:              // *
321:              // *  Based on contributions by
322:              // *     Jeff Rutter, Computer Science Division, University of California
323:              // *     at Berkeley, USA
324:              // *  Modified by Francoise Tisseur, University of Tennessee.
325:              // *
326:              // *  =====================================================================
327:              // *
328:              // *     .. Parameters ..
329:              // *     ..
330:              // *     .. Local Scalars ..
331:              // *     ..
332:              // *     .. External Functions ..
333:              // *     ..
334:              // *     .. External Subroutines ..
335:              // *     ..
336:              // *     .. Intrinsic Functions ..
337:              //      INTRINSIC          MAX, SIGN, SQRT;
338:              // *     ..
339:              // *     .. Executable Statements ..
340:              // *
341:              // *     Test the input parameters.
342:              // *
343:
344:              #endregion
345:
346:
347:              #region Body
348:
349:              INFO = 0;
350:              // *
351:              if (K < 0)
352:              {
353:                  INFO =  - 1;
354:              }
355:              else
356:              {
357:                  if (N < K)
358:                  {
359:                      INFO =  - 2;
360:                  }
361:                  else
362:                  {
363:                      if (LDQ < Math.Max(1, N))
364:                      {
365:                          INFO =  - 6;
366:                      }
367:                  }
368:              }
369:              if (INFO != 0)
370:              {
371:                  this._xerbla.Run("DLAED3",  - INFO);
372:                  return;
373:              }
374:              // *
375:              // *     Quick return if possible
376:              // *
377:              if (K == 0) return;
378:              // *
379:              // *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
380:              // *     be computed with high relative accuracy (barring over/underflow).
381:              // *     This is a problem on machines without a guard digit in
382:              // *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
383:              // *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
384:              // *     which on any of these machines zeros out the bottommost
385:              // *     bit of DLAMDA(I) if it is 1; this makes the subsequent
386:              // *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
387:              // *     occurs. On binary machines with a guard digit (almost all
388:              // *     machines) it does not change DLAMDA(I) at all. On hexadecimal
389:              // *     and decimal machines with a guard digit, it slightly
390:              // *     changes the bottommost bits of DLAMDA(I). It does not account
391:              // *     for hexadecimal or decimal machines without guard digits
392:              // *     (we know of none). We use a subroutine call to compute
393:              // *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
394:              // *     this code.
395:              // *
396:              for (I = 1; I <= K; I++)
397:              {
398:                  DLAMDA[I + o_dlamda] = this._dlamc3.Run(DLAMDA[I + o_dlamda], DLAMDA[I + o_dlamda]) - DLAMDA[I + o_dlamda];
399:              }
400:              // *
401:              for (J = 1; J <= K; J++)
402:              {
403:                  this._dlaed4.Run(K, J, DLAMDA, offset_dlamda, W, offset_w, ref Q, 1+J * LDQ + o_q, RHO
404:                                   , ref D[J + o_d], ref INFO);
405:                  // *
406:                  // *        If the zero finder fails, the computation is terminated.
407:                  // *
408:                  if (INFO != 0) goto LABEL120;
409:              }
410:              // *
411:              if (K == 1) goto LABEL110;
412:              if (K == 2)
413:              {
414:                  for (J = 1; J <= K; J++)
415:                  {
416:                      W[1 + o_w] = Q[1+J * LDQ + o_q];
417:                      W[2 + o_w] = Q[2+J * LDQ + o_q];
418:                      II = INDX[1 + o_indx];
419:                      Q[1+J * LDQ + o_q] = W[II + o_w];
420:                      II = INDX[2 + o_indx];
421:                      Q[2+J * LDQ + o_q] = W[II + o_w];
422:                  }
423:                  goto LABEL110;
424:              }
425:              // *
426:              // *     Compute updated W.
427:              // *
428:              this._dcopy.Run(K, W, offset_w, 1, ref S, offset_s, 1);
429:              // *
430:              // *     Initialize W(I) = Q(I,I)
431:              // *
432:              this._dcopy.Run(K, Q, offset_q, LDQ + 1, ref W, offset_w, 1);
433:              for (J = 1; J <= K; J++)
434:              {
435:                  for (I = 1; I <= J - 1; I++)
436:                  {
437:                      W[I + o_w] = W[I + o_w] * (Q[I+J * LDQ + o_q] / (DLAMDA[I + o_dlamda] - DLAMDA[J + o_dlamda]));
438:                  }
439:                  for (I = J + 1; I <= K; I++)
440:                  {
441:                      W[I + o_w] = W[I + o_w] * (Q[I+J * LDQ + o_q] / (DLAMDA[I + o_dlamda] - DLAMDA[J + o_dlamda]));
442:                  }
443:              }
444:              for (I = 1; I <= K; I++)
445:              {
446:                  W[I + o_w] = FortranLib.Sign(Math.Sqrt( - W[I + o_w]),S[I + o_s]);
447:              }
448:              // *
449:              // *     Compute eigenvectors of the modified rank-1 modification.
450:              // *
451:              for (J = 1; J <= K; J++)
452:              {
453:                  for (I = 1; I <= K; I++)
454:                  {
455:                      S[I + o_s] = W[I + o_w] / Q[I+J * LDQ + o_q];
456:                  }
457:                  TEMP = this._dnrm2.Run(K, S, offset_s, 1);
458:                  for (I = 1; I <= K; I++)
459:                  {
460:                      II = INDX[I + o_indx];
461:                      Q[I+J * LDQ + o_q] = S[II + o_s] / TEMP;
462:                  }
463:              }
464:              // *
465:              // *     Compute the updated eigenvectors.
466:              // *
467:          LABEL110:;
468:              // *
469:              N2 = N - N1;
470:              N12 = CTOT[1 + o_ctot] + CTOT[2 + o_ctot];
471:              N23 = CTOT[2 + o_ctot] + CTOT[3 + o_ctot];
472:              // *
473:              this._dlacpy.Run("A", N23, K, Q, CTOT[1 + o_ctot] + 1+1 * LDQ + o_q, LDQ, ref S, offset_s
474:                               , N23);
475:              IQ2 = N1 * N12 + 1;
476:              if (N23 != 0)
477:              {
478:                  this._dgemm.Run("N", "N", N2, K, N23, ONE
479:                                  , Q2, IQ2 + o_q2, N2, S, offset_s, N23, ZERO, ref Q, N1 + 1+1 * LDQ + o_q
480:                                  , LDQ);
481:              }
482:              else
483:              {
484:                  this._dlaset.Run("A", N2, K, ZERO, ZERO, ref Q, N1 + 1+1 * LDQ + o_q
485:                                   , LDQ);
486:              }
487:              // *
488:              this._dlacpy.Run("A", N12, K, Q, offset_q, LDQ, ref S, offset_s
489:                               , N12);
490:              if (N12 != 0)
491:              {
492:                  this._dgemm.Run("N", "N", N1, K, N12, ONE
493:                                  , Q2, offset_q2, N1, S, offset_s, N12, ZERO, ref Q, offset_q
494:                                  , LDQ);
495:              }
496:              else
497:              {
498:                  this._dlaset.Run("A", N1, K, ZERO, ZERO, ref Q, 1+1 * LDQ + o_q
499:                                   , LDQ);
500:              }
501:              // *
502:              // *
503:          LABEL120:;
504:              return;
505:              // *
506:              // *     End of DLAED3
507:              // *
508:
509:              #endregion
510:
511:          }
512:      }
513:  }