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

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