3rdparty/lapack/ssyr2.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int ssyr2_(char *uplo, integer *n, real *alpha, real *x,
   4         integer *incx, real *y, integer *incy, real *a, integer *lda)
   5 {
   6     /* System generated locals */
   7     integer a_dim1, a_offset, i__1, i__2;
   8
   9     /* Local variables */
  10     integer i__, j, ix, iy, jx, jy, kx, ky, info;
  11     real temp1, temp2;
  12     extern logical lsame_(char *, char *);
  13     extern /* Subroutine */ int xerbla_(char *, integer *);
  14
  15 /*     .. Scalar Arguments .. */
  16 /*     .. */
  17 /*     .. Array Arguments .. */
  18 /*     .. */
  19
  20 /*  Purpose */
  21 /*  ======= */
  22
  23 /*  SSYR2  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 /*  Arguments */
  31 /*  ========== */
  32
  33 /*  UPLO   - CHARACTER*1. */
  34 /*           On entry, UPLO specifies whether the upper or lower */
  35 /*           triangular part of the array A is to be referenced as */
  36 /*           follows: */
  37
  38 /*              UPLO = 'U' or 'u'   Only the upper triangular part of A */
  39 /*                                  is to be referenced. */
  40
  41 /*              UPLO = 'L' or 'l'   Only the lower triangular part of A */
  42 /*                                  is to be referenced. */
  43
  44 /*           Unchanged on exit. */
  45
  46 /*  N      - INTEGER. */
  47 /*           On entry, N specifies the order of the matrix A. */
  48 /*           N must be at least zero. */
  49 /*           Unchanged on exit. */
  50
  51 /*  ALPHA  - REAL            . */
  52 /*           On entry, ALPHA specifies the scalar alpha. */
  53 /*           Unchanged on exit. */
  54
  55 /*  X      - REAL             array of dimension at least */
  56 /*           ( 1 + ( n - 1 )*abs( INCX ) ). */
  57 /*           Before entry, the incremented array X must contain the n */
  58 /*           element vector x. */
  59 /*           Unchanged on exit. */
  60
  61 /*  INCX   - INTEGER. */
  62 /*           On entry, INCX specifies the increment for the elements of */
  63 /*           X. INCX must not be zero. */
  64 /*           Unchanged on exit. */
  65
  66 /*  Y      - REAL             array of dimension at least */
  67 /*           ( 1 + ( n - 1 )*abs( INCY ) ). */
  68 /*           Before entry, the incremented array Y must contain the n */
  69 /*           element vector y. */
  70 /*           Unchanged on exit. */
  71
  72 /*  INCY   - INTEGER. */
  73 /*           On entry, INCY specifies the increment for the elements of */
  74 /*           Y. INCY must not be zero. */
  75 /*           Unchanged on exit. */
  76
  77 /*  A      - REAL             array of DIMENSION ( LDA, n ). */
  78 /*           Before entry with  UPLO = 'U' or 'u', the leading n by n */
  79 /*           upper triangular part of the array A must contain the upper */
  80 /*           triangular part of the symmetric matrix and the strictly */
  81 /*           lower triangular part of A is not referenced. On exit, the */
  82 /*           upper triangular part of the array A is overwritten by the */
  83 /*           upper triangular part of the updated matrix. */
  84 /*           Before entry with UPLO = 'L' or 'l', the leading n by n */
  85 /*           lower triangular part of the array A must contain the lower */
  86 /*           triangular part of the symmetric matrix and the strictly */
  87 /*           upper triangular part of A is not referenced. On exit, the */
  88 /*           lower triangular part of the array A is overwritten by the */
  89 /*           lower triangular part of the updated matrix. */
  90
  91 /*  LDA    - INTEGER. */
  92 /*           On entry, LDA specifies the first dimension of A as declared */
  93 /*           in the calling (sub) program. LDA must be at least */
  94 /*           max( 1, n ). */
  95 /*           Unchanged on exit. */
  96
  97
  98 /*  Level 2 Blas routine. */
  99
 100 /*  -- Written on 22-October-1986. */
 101 /*     Jack Dongarra, Argonne National Lab. */
 102 /*     Jeremy Du Croz, Nag Central Office. */
 103 /*     Sven Hammarling, Nag Central Office. */
 104 /*     Richard Hanson, Sandia National Labs. */
 105
 106
 107 /*     .. Parameters .. */
 108 /*     .. */
 109 /*     .. Local Scalars .. */
 110 /*     .. */
 111 /*     .. External Functions .. */
 112 /*     .. */
 113 /*     .. External Subroutines .. */
 114 /*     .. */
 115 /*     .. Intrinsic Functions .. */
 116 /*     .. */
 117
 118 /*     Test the input parameters. */
 119
 120     /* Parameter adjustments */
 121     --x;
 122     --y;
 123     a_dim1 = *lda;
 124     a_offset = 1 + a_dim1;
 125     a -= a_offset;
 126
 127     /* Function Body */
 128     info = 0;
 129     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
 130         info = 1;
 131     } else if (*n < 0) {
 132         info = 2;
 133     } else if (*incx == 0) {
 134         info = 5;
 135     } else if (*incy == 0) {
 136         info = 7;
 137     } else if (*lda < max(1,*n)) {
 138         info = 9;
 139     }
 140     if (info != 0) {
 141         xerbla_("SSYR2 ", &info);
 142         return 0;
 143     }
 144
 145 /*     Quick return if possible. */
 146
 147     if (*n == 0 || *alpha == 0.f) {
 148         return 0;
 149     }
 150
 151 /*     Set up the start points in X and Y if the increments are not both */
 152 /*     unity. */
 153
 154     if (*incx != 1 || *incy != 1) {
 155         if (*incx > 0) {
 156             kx = 1;
 157         } else {
 158             kx = 1 - (*n - 1) * *incx;
 159         }
 160         if (*incy > 0) {
 161             ky = 1;
 162         } else {
 163             ky = 1 - (*n - 1) * *incy;
 164         }
 165         jx = kx;
 166         jy = ky;
 167     }
 168
 169 /*     Start the operations. In this version the elements of A are */
 170 /*     accessed sequentially with one pass through the triangular part */
 171 /*     of A. */
 172
 173     if (lsame_(uplo, "U")) {
 174
 175 /*        Form  A  when A is stored in the upper triangle. */
 176
 177         if (*incx == 1 && *incy == 1) {
 178             i__1 = *n;
 179             for (j = 1; j <= i__1; ++j) {
 180                 if (x[j] != 0.f || y[j] != 0.f) {
 181                     temp1 = *alpha * y[j];
 182                     temp2 = *alpha * x[j];
 183                     i__2 = j;
 184                     for (i__ = 1; i__ <= i__2; ++i__) {
 185                         a[i__ + j * a_dim1] = a[i__ + j * a_dim1] + x[i__] *
 186                                 temp1 + y[i__] * temp2;
 187 /* L10: */
 188                     }
 189                 }
 190 /* L20: */
 191             }
 192         } else {
 193             i__1 = *n;
 194             for (j = 1; j <= i__1; ++j) {
 195                 if (x[jx] != 0.f || y[jy] != 0.f) {
 196                     temp1 = *alpha * y[jy];
 197                     temp2 = *alpha * x[jx];
 198                     ix = kx;
 199                     iy = ky;
 200                     i__2 = j;
 201                     for (i__ = 1; i__ <= i__2; ++i__) {
 202                         a[i__ + j * a_dim1] = a[i__ + j * a_dim1] + x[ix] *
 203                                 temp1 + y[iy] * temp2;
 204                         ix += *incx;
 205                         iy += *incy;
 206 /* L30: */
 207                     }
 208                 }
 209                 jx += *incx;
 210                 jy += *incy;
 211 /* L40: */
 212             }
 213         }
 214     } else {
 215
 216 /*        Form  A  when A is stored in the lower triangle. */
 217
 218         if (*incx == 1 && *incy == 1) {
 219             i__1 = *n;
 220             for (j = 1; j <= i__1; ++j) {
 221                 if (x[j] != 0.f || y[j] != 0.f) {
 222                     temp1 = *alpha * y[j];
 223                     temp2 = *alpha * x[j];
 224                     i__2 = *n;
 225                     for (i__ = j; i__ <= i__2; ++i__) {
 226                         a[i__ + j * a_dim1] = a[i__ + j * a_dim1] + x[i__] *
 227                                 temp1 + y[i__] * temp2;
 228 /* L50: */
 229                     }
 230                 }
 231 /* L60: */
 232             }
 233         } else {
 234             i__1 = *n;
 235             for (j = 1; j <= i__1; ++j) {
 236                 if (x[jx] != 0.f || y[jy] != 0.f) {
 237                     temp1 = *alpha * y[jy];
 238                     temp2 = *alpha * x[jx];
 239                     ix = jx;
 240                     iy = jy;
 241                     i__2 = *n;
 242                     for (i__ = j; i__ <= i__2; ++i__) {
 243                         a[i__ + j * a_dim1] = a[i__ + j * a_dim1] + x[ix] *
 244                                 temp1 + y[iy] * temp2;
 245                         ix += *incx;
 246                         iy += *incy;
 247 /* L70: */
 248                     }
 249                 }
 250                 jx += *incx;
 251                 jy += *incy;
 252 /* L80: */
 253             }
 254         }
 255     }
 256
 257     return 0;
 258
 259 /*     End of SSYR2 . */
 260
 261 } /* ssyr2_ */