Skip Navigation Links
Numerical Libraries
Linear Algebra
Differential Equations
Optimization
Samples
Skip Navigation Links
Linear Algebra
CSLapack
CSBlas
   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:      /// DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
  27:      /// [  A   B  ]
  28:      /// [  B   C  ].
  29:      /// On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
  30:      /// eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
  31:      /// eigenvector for RT1, giving the decomposition
  32:      /// 
  33:      /// [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
  34:      /// [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
  35:      /// 
  36:      ///</summary>
  37:      public class DLAEV2
  38:      {
  39:      
  40:   
  41:          #region Fields
  42:          
  43:          const double ONE = 1.0E0; const double TWO = 2.0E0; const double ZERO = 0.0E0; const double HALF = 0.5E0; int SGN1 = 0; 
  44:          int SGN2 = 0;double AB = 0; double ACMN = 0; double ACMX = 0; double ACS = 0; double ADF = 0; double CS = 0; 
  45:          double CT = 0;double DF = 0; double RT = 0; double SM = 0; double TB = 0; double TN = 0; 
  46:   
  47:          #endregion
  48:   
  49:          public DLAEV2()
  50:          {
  51:      
  52:          }
  53:      
  54:          /// <summary>
  55:          /// Purpose
  56:          /// =======
  57:          /// 
  58:          /// DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
  59:          /// [  A   B  ]
  60:          /// [  B   C  ].
  61:          /// On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
  62:          /// eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
  63:          /// eigenvector for RT1, giving the decomposition
  64:          /// 
  65:          /// [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
  66:          /// [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
  67:          /// 
  68:          ///</summary>
  69:          /// <param name="A">
  70:          /// (input) DOUBLE PRECISION
  71:          /// The (1,1) element of the 2-by-2 matrix.
  72:          ///</param>
  73:          /// <param name="B">
  74:          /// (input) DOUBLE PRECISION
  75:          /// The (1,2) element and the conjugate of the (2,1) element of
  76:          /// the 2-by-2 matrix.
  77:          ///</param>
  78:          /// <param name="C">
  79:          /// (input) DOUBLE PRECISION
  80:          /// The (2,2) element of the 2-by-2 matrix.
  81:          ///</param>
  82:          /// <param name="RT1">
  83:          /// (output) DOUBLE PRECISION
  84:          /// The eigenvalue of larger absolute value.
  85:          ///</param>
  86:          /// <param name="RT2">
  87:          /// (output) DOUBLE PRECISION
  88:          /// The eigenvalue of smaller absolute value.
  89:          ///</param>
  90:          /// <param name="CS1">
  91:          /// (output) DOUBLE PRECISION
  92:          ///</param>
  93:          /// <param name="SN1">
  94:          /// (output) DOUBLE PRECISION
  95:          /// The vector (CS1, SN1) is a unit right eigenvector for RT1.
  96:          ///</param>
  97:          public void Run(double A, double B, double C, ref double RT1, ref double RT2, ref double CS1
  98:                           , ref double SN1)
  99:          {
 100:   
 101:              #region Prolog
 102:              
 103:              // *
 104:              // *  -- LAPACK auxiliary routine (version 3.1) --
 105:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 106:              // *     November 2006
 107:              // *
 108:              // *     .. Scalar Arguments ..
 109:              // *     ..
 110:              // *
 111:              // *  Purpose
 112:              // *  =======
 113:              // *
 114:              // *  DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
 115:              // *     [  A   B  ]
 116:              // *     [  B   C  ].
 117:              // *  On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
 118:              // *  eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
 119:              // *  eigenvector for RT1, giving the decomposition
 120:              // *
 121:              // *     [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
 122:              // *     [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
 123:              // *
 124:              // *  Arguments
 125:              // *  =========
 126:              // *
 127:              // *  A       (input) DOUBLE PRECISION
 128:              // *          The (1,1) element of the 2-by-2 matrix.
 129:              // *
 130:              // *  B       (input) DOUBLE PRECISION
 131:              // *          The (1,2) element and the conjugate of the (2,1) element of
 132:              // *          the 2-by-2 matrix.
 133:              // *
 134:              // *  C       (input) DOUBLE PRECISION
 135:              // *          The (2,2) element of the 2-by-2 matrix.
 136:              // *
 137:              // *  RT1     (output) DOUBLE PRECISION
 138:              // *          The eigenvalue of larger absolute value.
 139:              // *
 140:              // *  RT2     (output) DOUBLE PRECISION
 141:              // *          The eigenvalue of smaller absolute value.
 142:              // *
 143:              // *  CS1     (output) DOUBLE PRECISION
 144:              // *  SN1     (output) DOUBLE PRECISION
 145:              // *          The vector (CS1, SN1) is a unit right eigenvector for RT1.
 146:              // *
 147:              // *  Further Details
 148:              // *  ===============
 149:              // *
 150:              // *  RT1 is accurate to a few ulps barring over/underflow.
 151:              // *
 152:              // *  RT2 may be inaccurate if there is massive cancellation in the
 153:              // *  determinant A*C-B*B; higher precision or correctly rounded or
 154:              // *  correctly truncated arithmetic would be needed to compute RT2
 155:              // *  accurately in all cases.
 156:              // *
 157:              // *  CS1 and SN1 are accurate to a few ulps barring over/underflow.
 158:              // *
 159:              // *  Overflow is possible only if RT1 is within a factor of 5 of overflow.
 160:              // *  Underflow is harmless if the input data is 0 or exceeds
 161:              // *     underflow_threshold / macheps.
 162:              // *
 163:              // * =====================================================================
 164:              // *
 165:              // *     .. Parameters ..
 166:              // *     ..
 167:              // *     .. Local Scalars ..
 168:              // *     ..
 169:              // *     .. Intrinsic Functions ..
 170:              //      INTRINSIC          ABS, SQRT;
 171:              // *     ..
 172:              // *     .. Executable Statements ..
 173:              // *
 174:              // *     Compute the eigenvalues
 175:              // *
 176:   
 177:              #endregion
 178:   
 179:   
 180:              #region Body
 181:              
 182:              SM = A + C;
 183:              DF = A - C;
 184:              ADF = Math.Abs(DF);
 185:              TB = B + B;
 186:              AB = Math.Abs(TB);
 187:              if (Math.Abs(A) > Math.Abs(C))
 188:              {
 189:                  ACMX = A;
 190:                  ACMN = C;
 191:              }
 192:              else
 193:              {
 194:                  ACMX = C;
 195:                  ACMN = A;
 196:              }
 197:              if (ADF > AB)
 198:              {
 199:                  RT = ADF * Math.Sqrt(ONE + Math.Pow(AB / ADF,2));
 200:              }
 201:              else
 202:              {
 203:                  if (ADF < AB)
 204:                  {
 205:                      RT = AB * Math.Sqrt(ONE + Math.Pow(ADF / AB,2));
 206:                  }
 207:                  else
 208:                  {
 209:                      // *
 210:                      // *        Includes case AB=ADF=0
 211:                      // *
 212:                      RT = AB * Math.Sqrt(TWO);
 213:                  }
 214:              }
 215:              if (SM < ZERO)
 216:              {
 217:                  RT1 = HALF * (SM - RT);
 218:                  SGN1 =  - 1;
 219:                  // *
 220:                  // *        Order of execution important.
 221:                  // *        To get fully accurate smaller eigenvalue,
 222:                  // *        next line needs to be executed in higher precision.
 223:                  // *
 224:                  RT2 = (ACMX / RT1) * ACMN - (B / RT1) * B;
 225:              }
 226:              else
 227:              {
 228:                  if (SM > ZERO)
 229:                  {
 230:                      RT1 = HALF * (SM + RT);
 231:                      SGN1 = 1;
 232:                      // *
 233:                      // *        Order of execution important.
 234:                      // *        To get fully accurate smaller eigenvalue,
 235:                      // *        next line needs to be executed in higher precision.
 236:                      // *
 237:                      RT2 = (ACMX / RT1) * ACMN - (B / RT1) * B;
 238:                  }
 239:                  else
 240:                  {
 241:                      // *
 242:                      // *        Includes case RT1 = RT2 = 0
 243:                      // *
 244:                      RT1 = HALF * RT;
 245:                      RT2 =  - HALF * RT;
 246:                      SGN1 = 1;
 247:                  }
 248:              }
 249:              // *
 250:              // *     Compute the eigenvector
 251:              // *
 252:              if (DF >= ZERO)
 253:              {
 254:                  CS = DF + RT;
 255:                  SGN2 = 1;
 256:              }
 257:              else
 258:              {
 259:                  CS = DF - RT;
 260:                  SGN2 =  - 1;
 261:              }
 262:              ACS = Math.Abs(CS);
 263:              if (ACS > AB)
 264:              {
 265:                  CT =  - TB / CS;
 266:                  SN1 = ONE / Math.Sqrt(ONE + CT * CT);
 267:                  CS1 = CT * SN1;
 268:              }
 269:              else
 270:              {
 271:                  if (AB == ZERO)
 272:                  {
 273:                      CS1 = ONE;
 274:                      SN1 = ZERO;
 275:                  }
 276:                  else
 277:                  {
 278:                      TN =  - CS / TB;
 279:                      CS1 = ONE / Math.Sqrt(ONE + TN * TN);
 280:                      SN1 = TN * CS1;
 281:                  }
 282:              }
 283:              if (SGN1 == SGN2)
 284:              {
 285:                  TN = CS1;
 286:                  CS1 =  - SN1;
 287:                  SN1 = TN;
 288:              }
 289:              return;
 290:              // *
 291:              // *     End of DLAEV2
 292:              // *
 293:   
 294:              #endregion
 295:   
 296:          }
 297:      }
 298:  }