git.maemo.org Git - opencv/blob - 3rdparty/lapack/ssymv.c

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