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:      /// Purpose
  21:      /// =======
  22:      /// 
  23:      /// DSYR2  performs the symmetric rank 2 operation
  24:      /// 
  25:      /// A := alpha*x*y' + alpha*y*x' + A,
  26:      /// 
  27:      /// where alpha is a scalar, x and y are n element vectors and A is an n
  28:      /// by n symmetric matrix.
  29:      /// 
  30:      ///</summary>
  31:      public class DSYR2
  32:      {
  33:      
  34:   
  35:          #region Dependencies
  36:          
  37:          LSAME _lsame; XERBLA _xerbla; 
  38:   
  39:          #endregion
  40:   
  41:   
  42:          #region Fields
  43:          
  44:          const double ZERO = 0.0E+0; double TEMP1 = 0; double TEMP2 = 0; int I = 0; int INFO = 0; int IX = 0; int IY = 0; 
  45:          int J = 0;int JX = 0; int JY = 0; int KX = 0; int KY = 0; 
  46:   
  47:          #endregion
  48:   
  49:          public DSYR2(LSAME lsame, XERBLA xerbla)
  50:          {
  51:      
  52:   
  53:              #region Set Dependencies
  54:              
  55:              this._lsame = lsame; this._xerbla = xerbla; 
  56:   
  57:              #endregion
  58:   
  59:          }
  60:      
  61:          public DSYR2()
  62:          {
  63:      
  64:   
  65:              #region Dependencies (Initialization)
  66:              
  67:              LSAME lsame = new LSAME();
  68:              XERBLA xerbla = new XERBLA();
  69:   
  70:              #endregion
  71:   
  72:   
  73:              #region Set Dependencies
  74:              
  75:              this._lsame = lsame; this._xerbla = xerbla; 
  76:   
  77:              #endregion
  78:   
  79:          }
  80:          /// <summary>
  81:          /// Purpose
  82:          /// =======
  83:          /// 
  84:          /// DSYR2  performs the symmetric rank 2 operation
  85:          /// 
  86:          /// A := alpha*x*y' + alpha*y*x' + A,
  87:          /// 
  88:          /// where alpha is a scalar, x and y are n element vectors and A is an n
  89:          /// by n symmetric matrix.
  90:          /// 
  91:          ///</summary>
  92:          /// <param name="UPLO">
  93:          /// - CHARACTER*1.
  94:          /// On entry, UPLO specifies whether the upper or lower
  95:          /// triangular part of the array A is to be referenced as
  96:          /// follows:
  97:          /// 
  98:          /// UPLO = 'U' or 'u'   Only the upper triangular part of A
  99:          /// is to be referenced.
 100:          /// 
 101:          /// UPLO = 'L' or 'l'   Only the lower triangular part of A
 102:          /// is to be referenced.
 103:          /// 
 104:          /// Unchanged on exit.
 105:          ///</param>
 106:          /// <param name="N">
 107:          /// - INTEGER.
 108:          /// On entry, N specifies the order of the matrix A.
 109:          /// N must be at least zero.
 110:          /// Unchanged on exit.
 111:          ///</param>
 112:          /// <param name="ALPHA">
 113:          /// - DOUBLE PRECISION.
 114:          /// On entry, ALPHA specifies the scalar alpha.
 115:          /// Unchanged on exit.
 116:          ///</param>
 117:          /// <param name="X">
 118:          /// - DOUBLE PRECISION array of dimension at least
 119:          /// ( 1 + ( n - 1 )*abs( INCX ) ).
 120:          /// Before entry, the incremented array X must contain the n
 121:          /// element vector x.
 122:          /// Unchanged on exit.
 123:          ///</param>
 124:          /// <param name="INCX">
 125:          /// - INTEGER.
 126:          /// On entry, INCX specifies the increment for the elements of
 127:          /// X. INCX must not be zero.
 128:          /// Unchanged on exit.
 129:          ///</param>
 130:          /// <param name="Y">
 131:          /// - DOUBLE PRECISION array of dimension at least
 132:          /// ( 1 + ( n - 1 )*abs( INCY ) ).
 133:          /// Before entry, the incremented array Y must contain the n
 134:          /// element vector y.
 135:          /// Unchanged on exit.
 136:          ///</param>
 137:          /// <param name="INCY">
 138:          /// - INTEGER.
 139:          /// On entry, INCY specifies the increment for the elements of
 140:          /// Y. INCY must not be zero.
 141:          /// Unchanged on exit.
 142:          ///</param>
 143:          /// <param name="A">
 144:          /// := alpha*x*y' + alpha*y*x' + A,
 145:          ///</param>
 146:          /// <param name="LDA">
 147:          /// - INTEGER.
 148:          /// On entry, LDA specifies the first dimension of A as declared
 149:          /// in the calling (sub) program. LDA must be at least
 150:          /// max( 1, n ).
 151:          /// Unchanged on exit.
 152:          /// 
 153:          ///</param>
 154:          public void Run(string UPLO, int N, double ALPHA, double[] X, int offset_x, int INCX, double[] Y, int offset_y
 155:                           , int INCY, ref double[] A, int offset_a, int LDA)
 156:          {
 157:   
 158:              #region Array Index Correction
 159:              
 160:               int o_x = -1 + offset_x;  int o_y = -1 + offset_y;  int o_a = -1 - LDA + offset_a; 
 161:   
 162:              #endregion
 163:   
 164:   
 165:              #region Strings
 166:              
 167:              UPLO = UPLO.Substring(0, 1);  
 168:   
 169:              #endregion
 170:   
 171:   
 172:              #region Prolog
 173:              
 174:              // *     .. Scalar Arguments ..
 175:              // *     ..
 176:              // *     .. Array Arguments ..
 177:              // *     ..
 178:              // *
 179:              // *  Purpose
 180:              // *  =======
 181:              // *
 182:              // *  DSYR2  performs the symmetric rank 2 operation
 183:              // *
 184:              // *     A := alpha*x*y' + alpha*y*x' + A,
 185:              // *
 186:              // *  where alpha is a scalar, x and y are n element vectors and A is an n
 187:              // *  by n symmetric matrix.
 188:              // *
 189:              // *  Arguments
 190:              // *  ==========
 191:              // *
 192:              // *  UPLO   - CHARACTER*1.
 193:              // *           On entry, UPLO specifies whether the upper or lower
 194:              // *           triangular part of the array A is to be referenced as
 195:              // *           follows:
 196:              // *
 197:              // *              UPLO = 'U' or 'u'   Only the upper triangular part of A
 198:              // *                                  is to be referenced.
 199:              // *
 200:              // *              UPLO = 'L' or 'l'   Only the lower triangular part of A
 201:              // *                                  is to be referenced.
 202:              // *
 203:              // *           Unchanged on exit.
 204:              // *
 205:              // *  N      - INTEGER.
 206:              // *           On entry, N specifies the order of the matrix A.
 207:              // *           N must be at least zero.
 208:              // *           Unchanged on exit.
 209:              // *
 210:              // *  ALPHA  - DOUBLE PRECISION.
 211:              // *           On entry, ALPHA specifies the scalar alpha.
 212:              // *           Unchanged on exit.
 213:              // *
 214:              // *  X      - DOUBLE PRECISION array of dimension at least
 215:              // *           ( 1 + ( n - 1 )*abs( INCX ) ).
 216:              // *           Before entry, the incremented array X must contain the n
 217:              // *           element vector x.
 218:              // *           Unchanged on exit.
 219:              // *
 220:              // *  INCX   - INTEGER.
 221:              // *           On entry, INCX specifies the increment for the elements of
 222:              // *           X. INCX must not be zero.
 223:              // *           Unchanged on exit.
 224:              // *
 225:              // *  Y      - DOUBLE PRECISION array of dimension at least
 226:              // *           ( 1 + ( n - 1 )*abs( INCY ) ).
 227:              // *           Before entry, the incremented array Y must contain the n
 228:              // *           element vector y.
 229:              // *           Unchanged on exit.
 230:              // *
 231:              // *  INCY   - INTEGER.
 232:              // *           On entry, INCY specifies the increment for the elements of
 233:              // *           Y. INCY must not be zero.
 234:              // *           Unchanged on exit.
 235:              // *
 236:              // *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
 237:              // *           Before entry with  UPLO = 'U' or 'u', the leading n by n
 238:              // *           upper triangular part of the array A must contain the upper
 239:              // *           triangular part of the symmetric matrix and the strictly
 240:              // *           lower triangular part of A is not referenced. On exit, the
 241:              // *           upper triangular part of the array A is overwritten by the
 242:              // *           upper triangular part of the updated matrix.
 243:              // *           Before entry with UPLO = 'L' or 'l', the leading n by n
 244:              // *           lower triangular part of the array A must contain the lower
 245:              // *           triangular part of the symmetric matrix and the strictly
 246:              // *           upper triangular part of A is not referenced. On exit, the
 247:              // *           lower triangular part of the array A is overwritten by the
 248:              // *           lower triangular part of the updated matrix.
 249:              // *
 250:              // *  LDA    - INTEGER.
 251:              // *           On entry, LDA specifies the first dimension of A as declared
 252:              // *           in the calling (sub) program. LDA must be at least
 253:              // *           max( 1, n ).
 254:              // *           Unchanged on exit.
 255:              // *
 256:              // *
 257:              // *  Level 2 Blas routine.
 258:              // *
 259:              // *  -- Written on 22-October-1986.
 260:              // *     Jack Dongarra, Argonne National Lab.
 261:              // *     Jeremy Du Croz, Nag Central Office.
 262:              // *     Sven Hammarling, Nag Central Office.
 263:              // *     Richard Hanson, Sandia National Labs.
 264:              // *
 265:              // *
 266:              // *     .. Parameters ..
 267:              // *     ..
 268:              // *     .. Local Scalars ..
 269:              // *     ..
 270:              // *     .. External Functions ..
 271:              // *     ..
 272:              // *     .. External Subroutines ..
 273:              // *     ..
 274:              // *     .. Intrinsic Functions ..
 275:              //      INTRINSIC MAX;
 276:              // *     ..
 277:              // *
 278:              // *     Test the input parameters.
 279:              // *
 280:   
 281:              #endregion
 282:   
 283:   
 284:              #region Body
 285:              
 286:              INFO = 0;
 287:              if (!this._lsame.Run(UPLO, "U") && !this._lsame.Run(UPLO, "L"))
 288:              {
 289:                  INFO = 1;
 290:              }
 291:              else
 292:              {
 293:                  if (N < 0)
 294:                  {
 295:                      INFO = 2;
 296:                  }
 297:                  else
 298:                  {
 299:                      if (INCX == 0)
 300:                      {
 301:                          INFO = 5;
 302:                      }
 303:                      else
 304:                      {
 305:                          if (INCY == 0)
 306:                          {
 307:                              INFO = 7;
 308:                          }
 309:                          else
 310:                          {
 311:                              if (LDA < Math.Max(1, N))
 312:                              {
 313:                                  INFO = 9;
 314:                              }
 315:                          }
 316:                      }
 317:                  }
 318:              }
 319:              if (INFO != 0)
 320:              {
 321:                  this._xerbla.Run("DSYR2 ", INFO);
 322:                  return;
 323:              }
 324:              // *
 325:              // *     Quick return if possible.
 326:              // *
 327:              if ((N == 0) || (ALPHA == ZERO)) return;
 328:              // *
 329:              // *     Set up the start points in X and Y if the increments are not both
 330:              // *     unity.
 331:              // *
 332:              if ((INCX != 1) || (INCY != 1))
 333:              {
 334:                  if (INCX > 0)
 335:                  {
 336:                      KX = 1;
 337:                  }
 338:                  else
 339:                  {
 340:                      KX = 1 - (N - 1) * INCX;
 341:                  }
 342:                  if (INCY > 0)
 343:                  {
 344:                      KY = 1;
 345:                  }
 346:                  else
 347:                  {
 348:                      KY = 1 - (N - 1) * INCY;
 349:                  }
 350:                  JX = KX;
 351:                  JY = KY;
 352:              }
 353:              // *
 354:              // *     Start the operations. In this version the elements of A are
 355:              // *     accessed sequentially with one pass through the triangular part
 356:              // *     of A.
 357:              // *
 358:              if (this._lsame.Run(UPLO, "U"))
 359:              {
 360:                  // *
 361:                  // *        Form  A  when A is stored in the upper triangle.
 362:                  // *
 363:                  if ((INCX == 1) && (INCY == 1))
 364:                  {
 365:                      for (J = 1; J <= N; J++)
 366:                      {
 367:                          if ((X[J + o_x] != ZERO) || (Y[J + o_y] != ZERO))
 368:                          {
 369:                              TEMP1 = ALPHA * Y[J + o_y];
 370:                              TEMP2 = ALPHA * X[J + o_x];
 371:                              for (I = 1; I <= J; I++)
 372:                              {
 373:                                  A[I+J * LDA + o_a] = A[I+J * LDA + o_a] + X[I + o_x] * TEMP1 + Y[I + o_y] * TEMP2;
 374:                              }
 375:                          }
 376:                      }
 377:                  }
 378:                  else
 379:                  {
 380:                      for (J = 1; J <= N; J++)
 381:                      {
 382:                          if ((X[JX + o_x] != ZERO) || (Y[JY + o_y] != ZERO))
 383:                          {
 384:                              TEMP1 = ALPHA * Y[JY + o_y];
 385:                              TEMP2 = ALPHA * X[JX + o_x];
 386:                              IX = KX;
 387:                              IY = KY;
 388:                              for (I = 1; I <= J; I++)
 389:                              {
 390:                                  A[I+J * LDA + o_a] = A[I+J * LDA + o_a] + X[IX + o_x] * TEMP1 + Y[IY + o_y] * TEMP2;
 391:                                  IX = IX + INCX;
 392:                                  IY = IY + INCY;
 393:                              }
 394:                          }
 395:                          JX = JX + INCX;
 396:                          JY = JY + INCY;
 397:                      }
 398:                  }
 399:              }
 400:              else
 401:              {
 402:                  // *
 403:                  // *        Form  A  when A is stored in the lower triangle.
 404:                  // *
 405:                  if ((INCX == 1) && (INCY == 1))
 406:                  {
 407:                      for (J = 1; J <= N; J++)
 408:                      {
 409:                          if ((X[J + o_x] != ZERO) || (Y[J + o_y] != ZERO))
 410:                          {
 411:                              TEMP1 = ALPHA * Y[J + o_y];
 412:                              TEMP2 = ALPHA * X[J + o_x];
 413:                              for (I = J; I <= N; I++)
 414:                              {
 415:                                  A[I+J * LDA + o_a] = A[I+J * LDA + o_a] + X[I + o_x] * TEMP1 + Y[I + o_y] * TEMP2;
 416:                              }
 417:                          }
 418:                      }
 419:                  }
 420:                  else
 421:                  {
 422:                      for (J = 1; J <= N; J++)
 423:                      {
 424:                          if ((X[JX + o_x] != ZERO) || (Y[JY + o_y] != ZERO))
 425:                          {
 426:                              TEMP1 = ALPHA * Y[JY + o_y];
 427:                              TEMP2 = ALPHA * X[JX + o_x];
 428:                              IX = JX;
 429:                              IY = JY;
 430:                              for (I = J; I <= N; I++)
 431:                              {
 432:                                  A[I+J * LDA + o_a] = A[I+J * LDA + o_a] + X[IX + o_x] * TEMP1 + Y[IY + o_y] * TEMP2;
 433:                                  IX = IX + INCX;
 434:                                  IY = IY + INCY;
 435:                              }
 436:                          }
 437:                          JX = JX + INCX;
 438:                          JY = JY + INCY;
 439:                      }
 440:                  }
 441:              }
 442:              // *
 443:              return;
 444:              // *
 445:              // *     End of DSYR2 .
 446:              // *
 447:   
 448:              #endregion
 449:   
 450:          }
 451:      }
 452:  }