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

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