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 routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DLASQ1 computes the singular values of a real N-by-N bidiagonal
  27:      /// matrix with diagonal D and off-diagonal E. The singular values
  28:      /// are computed to high relative accuracy, in the absence of
  29:      /// denormalization, underflow and overflow. The algorithm was first
  30:      /// presented in
  31:      /// 
  32:      /// "Accurate singular values and differential qd algorithms" by K. V.
  33:      /// Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
  34:      /// 1994,
  35:      /// 
  36:      /// and the present implementation is described in "An implementation of
  37:      /// the dqds Algorithm (Positive Case)", LAPACK Working Note.
  38:      /// 
  39:      ///</summary>
  40:      public class DLASQ1
  41:      {
  42:      
  43:   
  44:          #region Dependencies
  45:          
  46:          DCOPY _dcopy; DLAS2 _dlas2; DLASCL _dlascl; DLASQ2 _dlasq2; DLASRT _dlasrt; XERBLA _xerbla; DLAMCH _dlamch; 
  47:   
  48:          #endregion
  49:   
  50:   
  51:          #region Fields
  52:          
  53:          const double ZERO = 0.0E0; int I = 0; int IINFO = 0; double EPS = 0; double SCALE = 0; double SAFMIN = 0; 
  54:          double SIGMN = 0;double SIGMX = 0; 
  55:   
  56:          #endregion
  57:   
  58:          public DLASQ1(DCOPY dcopy, DLAS2 dlas2, DLASCL dlascl, DLASQ2 dlasq2, DLASRT dlasrt, XERBLA xerbla, DLAMCH dlamch)
  59:          {
  60:      
  61:   
  62:              #region Set Dependencies
  63:              
  64:              this._dcopy = dcopy; this._dlas2 = dlas2; this._dlascl = dlascl; this._dlasq2 = dlasq2; this._dlasrt = dlasrt; 
  65:              this._xerbla = xerbla;this._dlamch = dlamch; 
  66:   
  67:              #endregion
  68:   
  69:          }
  70:      
  71:          public DLASQ1()
  72:          {
  73:      
  74:   
  75:              #region Dependencies (Initialization)
  76:              
  77:              DCOPY dcopy = new DCOPY();
  78:              DLAS2 dlas2 = new DLAS2();
  79:              LSAME lsame = new LSAME();
  80:              DLAMC3 dlamc3 = new DLAMC3();
  81:              XERBLA xerbla = new XERBLA();
  82:              DLASQ5 dlasq5 = new DLASQ5();
  83:              DLAZQ4 dlazq4 = new DLAZQ4();
  84:              IEEECK ieeeck = new IEEECK();
  85:              IPARMQ iparmq = new IPARMQ();
  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:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
  92:              DLASQ6 dlasq6 = new DLASQ6(dlamch);
  93:              DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
  94:              DLASRT dlasrt = new DLASRT(lsame, xerbla);
  95:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  96:              DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
  97:   
  98:              #endregion
  99:   
 100:   
 101:              #region Set Dependencies
 102:              
 103:              this._dcopy = dcopy; this._dlas2 = dlas2; this._dlascl = dlascl; this._dlasq2 = dlasq2; this._dlasrt = dlasrt; 
 104:              this._xerbla = xerbla;this._dlamch = dlamch; 
 105:   
 106:              #endregion
 107:   
 108:          }
 109:          /// <summary>
 110:          /// Purpose
 111:          /// =======
 112:          /// 
 113:          /// DLASQ1 computes the singular values of a real N-by-N bidiagonal
 114:          /// matrix with diagonal D and off-diagonal E. The singular values
 115:          /// are computed to high relative accuracy, in the absence of
 116:          /// denormalization, underflow and overflow. The algorithm was first
 117:          /// presented in
 118:          /// 
 119:          /// "Accurate singular values and differential qd algorithms" by K. V.
 120:          /// Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
 121:          /// 1994,
 122:          /// 
 123:          /// and the present implementation is described in "An implementation of
 124:          /// the dqds Algorithm (Positive Case)", LAPACK Working Note.
 125:          /// 
 126:          ///</summary>
 127:          /// <param name="N">
 128:          /// (input) INTEGER
 129:          /// The number of rows and columns in the matrix. N .GE. 0.
 130:          ///</param>
 131:          /// <param name="D">
 132:          /// (input/output) DOUBLE PRECISION array, dimension (N)
 133:          /// On entry, D contains the diagonal elements of the
 134:          /// bidiagonal matrix whose SVD is desired. On normal exit,
 135:          /// D contains the singular values in decreasing order.
 136:          ///</param>
 137:          /// <param name="E">
 138:          /// (input/output) DOUBLE PRECISION array, dimension (N)
 139:          /// On entry, elements E(1:N-1) contain the off-diagonal elements
 140:          /// of the bidiagonal matrix whose SVD is desired.
 141:          /// On exit, E is overwritten.
 142:          ///</param>
 143:          /// <param name="WORK">
 144:          /// (workspace) DOUBLE PRECISION array, dimension (4*N)
 145:          ///</param>
 146:          /// <param name="INFO">
 147:          /// (output) INTEGER
 148:          /// = 0: successful exit
 149:          /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
 150:          /// .GT. 0: the algorithm failed
 151:          /// = 1, a split was marked by a positive value in E
 152:          /// = 2, current block of Z not diagonalized after 30*N
 153:          /// iterations (in inner while loop)
 154:          /// = 3, termination criterion of outer while loop not met 
 155:          /// (program created more than N unreduced blocks)
 156:          ///</param>
 157:          public void Run(int N, ref double[] D, int offset_d, double[] E, int offset_e, ref double[] WORK, int offset_work, ref int INFO)
 158:          {
 159:   
 160:              #region Array Index Correction
 161:              
 162:               int o_d = -1 + offset_d;  int o_e = -1 + offset_e;  int o_work = -1 + offset_work; 
 163:   
 164:              #endregion
 165:   
 166:   
 167:              #region Prolog
 168:              
 169:              // *
 170:              // *  -- LAPACK routine (version 3.1) --
 171:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 172:              // *     November 2006
 173:              // *
 174:              // *     .. Scalar Arguments ..
 175:              // *     ..
 176:              // *     .. Array Arguments ..
 177:              // *     ..
 178:              // *
 179:              // *  Purpose
 180:              // *  =======
 181:              // *
 182:              // *  DLASQ1 computes the singular values of a real N-by-N bidiagonal
 183:              // *  matrix with diagonal D and off-diagonal E. The singular values
 184:              // *  are computed to high relative accuracy, in the absence of
 185:              // *  denormalization, underflow and overflow. The algorithm was first
 186:              // *  presented in
 187:              // *
 188:              // *  "Accurate singular values and differential qd algorithms" by K. V.
 189:              // *  Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
 190:              // *  1994,
 191:              // *
 192:              // *  and the present implementation is described in "An implementation of
 193:              // *  the dqds Algorithm (Positive Case)", LAPACK Working Note.
 194:              // *
 195:              // *  Arguments
 196:              // *  =========
 197:              // *
 198:              // *  N     (input) INTEGER
 199:              // *        The number of rows and columns in the matrix. N >= 0.
 200:              // *
 201:              // *  D     (input/output) DOUBLE PRECISION array, dimension (N)
 202:              // *        On entry, D contains the diagonal elements of the
 203:              // *        bidiagonal matrix whose SVD is desired. On normal exit,
 204:              // *        D contains the singular values in decreasing order.
 205:              // *
 206:              // *  E     (input/output) DOUBLE PRECISION array, dimension (N)
 207:              // *        On entry, elements E(1:N-1) contain the off-diagonal elements
 208:              // *        of the bidiagonal matrix whose SVD is desired.
 209:              // *        On exit, E is overwritten.
 210:              // *
 211:              // *  WORK  (workspace) DOUBLE PRECISION array, dimension (4*N)
 212:              // *
 213:              // *  INFO  (output) INTEGER
 214:              // *        = 0: successful exit
 215:              // *        < 0: if INFO = -i, the i-th argument had an illegal value
 216:              // *        > 0: the algorithm failed
 217:              // *             = 1, a split was marked by a positive value in E
 218:              // *             = 2, current block of Z not diagonalized after 30*N
 219:              // *                  iterations (in inner while loop)
 220:              // *             = 3, termination criterion of outer while loop not met 
 221:              // *                  (program created more than N unreduced blocks)
 222:              // *
 223:              // *  =====================================================================
 224:              // *
 225:              // *     .. Parameters ..
 226:              // *     ..
 227:              // *     .. Local Scalars ..
 228:              // *     ..
 229:              // *     .. External Subroutines ..
 230:              // *     ..
 231:              // *     .. External Functions ..
 232:              // *     ..
 233:              // *     .. Intrinsic Functions ..
 234:              //      INTRINSIC          ABS, MAX, SQRT;
 235:              // *     ..
 236:              // *     .. Executable Statements ..
 237:              // *
 238:   
 239:              #endregion
 240:   
 241:   
 242:              #region Body
 243:              
 244:              INFO = 0;
 245:              if (N < 0)
 246:              {
 247:                  INFO =  - 2;
 248:                  this._xerbla.Run("DLASQ1",  - INFO);
 249:                  return;
 250:              }
 251:              else
 252:              {
 253:                  if (N == 0)
 254:                  {
 255:                      return;
 256:                  }
 257:                  else
 258:                  {
 259:                      if (N == 1)
 260:                      {
 261:                          D[1 + o_d] = Math.Abs(D[1 + o_d]);
 262:                          return;
 263:                      }
 264:                      else
 265:                      {
 266:                          if (N == 2)
 267:                          {
 268:                              this._dlas2.Run(D[1 + o_d], E[1 + o_e], D[2 + o_d], ref SIGMN, ref SIGMX);
 269:                              D[1 + o_d] = SIGMX;
 270:                              D[2 + o_d] = SIGMN;
 271:                              return;
 272:                          }
 273:                      }
 274:                  }
 275:              }
 276:              // *
 277:              // *     Estimate the largest singular value.
 278:              // *
 279:              SIGMX = ZERO;
 280:              for (I = 1; I <= N - 1; I++)
 281:              {
 282:                  D[I + o_d] = Math.Abs(D[I + o_d]);
 283:                  SIGMX = Math.Max(SIGMX, Math.Abs(E[I + o_e]));
 284:              }
 285:              D[N + o_d] = Math.Abs(D[N + o_d]);
 286:              // *
 287:              // *     Early return if SIGMX is zero (matrix is already diagonal).
 288:              // *
 289:              if (SIGMX == ZERO)
 290:              {
 291:                  this._dlasrt.Run("D", N, ref D, offset_d, ref IINFO);
 292:                  return;
 293:              }
 294:              // *
 295:              for (I = 1; I <= N; I++)
 296:              {
 297:                  SIGMX = Math.Max(SIGMX, D[I + o_d]);
 298:              }
 299:              // *
 300:              // *     Copy D and E into WORK (in the Z format) and scale (squaring the
 301:              // *     input data makes scaling by a power of the radix pointless).
 302:              // *
 303:              EPS = this._dlamch.Run("Precision");
 304:              SAFMIN = this._dlamch.Run("Safe minimum");
 305:              SCALE = Math.Sqrt(EPS / SAFMIN);
 306:              this._dcopy.Run(N, D, offset_d, 1, ref WORK, 1 + o_work, 2);
 307:              this._dcopy.Run(N - 1, E, offset_e, 1, ref WORK, 2 + o_work, 2);
 308:              this._dlascl.Run("G", 0, 0, SIGMX, SCALE, 2 * N - 1
 309:                               , 1, ref WORK, offset_work, 2 * N - 1, ref IINFO);
 310:              // *         
 311:              // *     Compute the q's and e's.
 312:              // *
 313:              for (I = 1; I <= 2 * N - 1; I++)
 314:              {
 315:                  WORK[I + o_work] = Math.Pow(WORK[I + o_work],2);
 316:              }
 317:              WORK[2 * N + o_work] = ZERO;
 318:              // *
 319:              this._dlasq2.Run(N, ref WORK, offset_work, ref INFO);
 320:              // *
 321:              if (INFO == 0)
 322:              {
 323:                  for (I = 1; I <= N; I++)
 324:                  {
 325:                      D[I + o_d] = Math.Sqrt(WORK[I + o_work]);
 326:                  }
 327:                  this._dlascl.Run("G", 0, 0, SCALE, SIGMX, N
 328:                                   , 1, ref D, offset_d, N, ref IINFO);
 329:              }
 330:              // *
 331:              return;
 332:              // *
 333:              // *     End of DLASQ1
 334:              // *
 335:   
 336:              #endregion
 337:   
 338:          }
 339:      }
 340:  }