dtrmv.c
Go to the documentation of this file.
00001 /* dtrmv.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int dtrmv_(char *uplo, char *trans, char *diag, integer *n, 
00017         doublereal *a, integer *lda, doublereal *x, integer *incx)
00018 {
00019     /* System generated locals */
00020     integer a_dim1, a_offset, i__1, i__2;
00021 
00022     /* Local variables */
00023     integer i__, j, ix, jx, kx, info;
00024     doublereal temp;
00025     extern logical lsame_(char *, char *);
00026     extern /* Subroutine */ int xerbla_(char *, integer *);
00027     logical nounit;
00028 
00029 /*     .. Scalar Arguments .. */
00030 /*     .. */
00031 /*     .. Array Arguments .. */
00032 /*     .. */
00033 
00034 /*  Purpose */
00035 /*  ======= */
00036 
00037 /*  DTRMV  performs one of the matrix-vector operations */
00038 
00039 /*     x := A*x,   or   x := A'*x, */
00040 
00041 /*  where x is an n element vector and  A is an n by n unit, or non-unit, */
00042 /*  upper or lower triangular matrix. */
00043 
00044 /*  Arguments */
00045 /*  ========== */
00046 
00047 /*  UPLO   - CHARACTER*1. */
00048 /*           On entry, UPLO specifies whether the matrix is an upper or */
00049 /*           lower triangular matrix as follows: */
00050 
00051 /*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
00052 
00053 /*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
00054 
00055 /*           Unchanged on exit. */
00056 
00057 /*  TRANS  - CHARACTER*1. */
00058 /*           On entry, TRANS specifies the operation to be performed as */
00059 /*           follows: */
00060 
00061 /*              TRANS = 'N' or 'n'   x := A*x. */
00062 
00063 /*              TRANS = 'T' or 't'   x := A'*x. */
00064 
00065 /*              TRANS = 'C' or 'c'   x := A'*x. */
00066 
00067 /*           Unchanged on exit. */
00068 
00069 /*  DIAG   - CHARACTER*1. */
00070 /*           On entry, DIAG specifies whether or not A is unit */
00071 /*           triangular as follows: */
00072 
00073 /*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
00074 
00075 /*              DIAG = 'N' or 'n'   A is not assumed to be unit */
00076 /*                                  triangular. */
00077 
00078 /*           Unchanged on exit. */
00079 
00080 /*  N      - INTEGER. */
00081 /*           On entry, N specifies the order of the matrix A. */
00082 /*           N must be at least zero. */
00083 /*           Unchanged on exit. */
00084 
00085 /*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
00086 /*           Before entry with  UPLO = 'U' or 'u', the leading n by n */
00087 /*           upper triangular part of the array A must contain the upper */
00088 /*           triangular matrix and the strictly lower triangular part of */
00089 /*           A is not referenced. */
00090 /*           Before entry with UPLO = 'L' or 'l', the leading n by n */
00091 /*           lower triangular part of the array A must contain the lower */
00092 /*           triangular matrix and the strictly upper triangular part of */
00093 /*           A is not referenced. */
00094 /*           Note that when  DIAG = 'U' or 'u', the diagonal elements of */
00095 /*           A are not referenced either, but are assumed to be unity. */
00096 /*           Unchanged on exit. */
00097 
00098 /*  LDA    - INTEGER. */
00099 /*           On entry, LDA specifies the first dimension of A as declared */
00100 /*           in the calling (sub) program. LDA must be at least */
00101 /*           max( 1, n ). */
00102 /*           Unchanged on exit. */
00103 
00104 /*  X      - DOUBLE PRECISION array of dimension at least */
00105 /*           ( 1 + ( n - 1 )*abs( INCX ) ). */
00106 /*           Before entry, the incremented array X must contain the n */
00107 /*           element vector x. On exit, X is overwritten with the */
00108 /*           tranformed vector x. */
00109 
00110 /*  INCX   - INTEGER. */
00111 /*           On entry, INCX specifies the increment for the elements of */
00112 /*           X. INCX must not be zero. */
00113 /*           Unchanged on exit. */
00114 
00115 
00116 /*  Level 2 Blas routine. */
00117 
00118 /*  -- Written on 22-October-1986. */
00119 /*     Jack Dongarra, Argonne National Lab. */
00120 /*     Jeremy Du Croz, Nag Central Office. */
00121 /*     Sven Hammarling, Nag Central Office. */
00122 /*     Richard Hanson, Sandia National Labs. */
00123 
00124 
00125 /*     .. Parameters .. */
00126 /*     .. */
00127 /*     .. Local Scalars .. */
00128 /*     .. */
00129 /*     .. External Functions .. */
00130 /*     .. */
00131 /*     .. External Subroutines .. */
00132 /*     .. */
00133 /*     .. Intrinsic Functions .. */
00134 /*     .. */
00135 
00136 /*     Test the input parameters. */
00137 
00138     /* Parameter adjustments */
00139     a_dim1 = *lda;
00140     a_offset = 1 + a_dim1;
00141     a -= a_offset;
00142     --x;
00143 
00144     /* Function Body */
00145     info = 0;
00146     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00147         info = 1;
00148     } else if (! lsame_(trans, "N") && ! lsame_(trans, 
00149             "T") && ! lsame_(trans, "C")) {
00150         info = 2;
00151     } else if (! lsame_(diag, "U") && ! lsame_(diag, 
00152             "N")) {
00153         info = 3;
00154     } else if (*n < 0) {
00155         info = 4;
00156     } else if (*lda < max(1,*n)) {
00157         info = 6;
00158     } else if (*incx == 0) {
00159         info = 8;
00160     }
00161     if (info != 0) {
00162         xerbla_("DTRMV ", &info);
00163         return 0;
00164     }
00165 
00166 /*     Quick return if possible. */
00167 
00168     if (*n == 0) {
00169         return 0;
00170     }
00171 
00172     nounit = lsame_(diag, "N");
00173 
00174 /*     Set up the start point in X if the increment is not unity. This */
00175 /*     will be  ( N - 1 )*INCX  too small for descending loops. */
00176 
00177     if (*incx <= 0) {
00178         kx = 1 - (*n - 1) * *incx;
00179     } else if (*incx != 1) {
00180         kx = 1;
00181     }
00182 
00183 /*     Start the operations. In this version the elements of A are */
00184 /*     accessed sequentially with one pass through A. */
00185 
00186     if (lsame_(trans, "N")) {
00187 
00188 /*        Form  x := A*x. */
00189 
00190         if (lsame_(uplo, "U")) {
00191             if (*incx == 1) {
00192                 i__1 = *n;
00193                 for (j = 1; j <= i__1; ++j) {
00194                     if (x[j] != 0.) {
00195                         temp = x[j];
00196                         i__2 = j - 1;
00197                         for (i__ = 1; i__ <= i__2; ++i__) {
00198                             x[i__] += temp * a[i__ + j * a_dim1];
00199 /* L10: */
00200                         }
00201                         if (nounit) {
00202                             x[j] *= a[j + j * a_dim1];
00203                         }
00204                     }
00205 /* L20: */
00206                 }
00207             } else {
00208                 jx = kx;
00209                 i__1 = *n;
00210                 for (j = 1; j <= i__1; ++j) {
00211                     if (x[jx] != 0.) {
00212                         temp = x[jx];
00213                         ix = kx;
00214                         i__2 = j - 1;
00215                         for (i__ = 1; i__ <= i__2; ++i__) {
00216                             x[ix] += temp * a[i__ + j * a_dim1];
00217                             ix += *incx;
00218 /* L30: */
00219                         }
00220                         if (nounit) {
00221                             x[jx] *= a[j + j * a_dim1];
00222                         }
00223                     }
00224                     jx += *incx;
00225 /* L40: */
00226                 }
00227             }
00228         } else {
00229             if (*incx == 1) {
00230                 for (j = *n; j >= 1; --j) {
00231                     if (x[j] != 0.) {
00232                         temp = x[j];
00233                         i__1 = j + 1;
00234                         for (i__ = *n; i__ >= i__1; --i__) {
00235                             x[i__] += temp * a[i__ + j * a_dim1];
00236 /* L50: */
00237                         }
00238                         if (nounit) {
00239                             x[j] *= a[j + j * a_dim1];
00240                         }
00241                     }
00242 /* L60: */
00243                 }
00244             } else {
00245                 kx += (*n - 1) * *incx;
00246                 jx = kx;
00247                 for (j = *n; j >= 1; --j) {
00248                     if (x[jx] != 0.) {
00249                         temp = x[jx];
00250                         ix = kx;
00251                         i__1 = j + 1;
00252                         for (i__ = *n; i__ >= i__1; --i__) {
00253                             x[ix] += temp * a[i__ + j * a_dim1];
00254                             ix -= *incx;
00255 /* L70: */
00256                         }
00257                         if (nounit) {
00258                             x[jx] *= a[j + j * a_dim1];
00259                         }
00260                     }
00261                     jx -= *incx;
00262 /* L80: */
00263                 }
00264             }
00265         }
00266     } else {
00267 
00268 /*        Form  x := A'*x. */
00269 
00270         if (lsame_(uplo, "U")) {
00271             if (*incx == 1) {
00272                 for (j = *n; j >= 1; --j) {
00273                     temp = x[j];
00274                     if (nounit) {
00275                         temp *= a[j + j * a_dim1];
00276                     }
00277                     for (i__ = j - 1; i__ >= 1; --i__) {
00278                         temp += a[i__ + j * a_dim1] * x[i__];
00279 /* L90: */
00280                     }
00281                     x[j] = temp;
00282 /* L100: */
00283                 }
00284             } else {
00285                 jx = kx + (*n - 1) * *incx;
00286                 for (j = *n; j >= 1; --j) {
00287                     temp = x[jx];
00288                     ix = jx;
00289                     if (nounit) {
00290                         temp *= a[j + j * a_dim1];
00291                     }
00292                     for (i__ = j - 1; i__ >= 1; --i__) {
00293                         ix -= *incx;
00294                         temp += a[i__ + j * a_dim1] * x[ix];
00295 /* L110: */
00296                     }
00297                     x[jx] = temp;
00298                     jx -= *incx;
00299 /* L120: */
00300                 }
00301             }
00302         } else {
00303             if (*incx == 1) {
00304                 i__1 = *n;
00305                 for (j = 1; j <= i__1; ++j) {
00306                     temp = x[j];
00307                     if (nounit) {
00308                         temp *= a[j + j * a_dim1];
00309                     }
00310                     i__2 = *n;
00311                     for (i__ = j + 1; i__ <= i__2; ++i__) {
00312                         temp += a[i__ + j * a_dim1] * x[i__];
00313 /* L130: */
00314                     }
00315                     x[j] = temp;
00316 /* L140: */
00317                 }
00318             } else {
00319                 jx = kx;
00320                 i__1 = *n;
00321                 for (j = 1; j <= i__1; ++j) {
00322                     temp = x[jx];
00323                     ix = jx;
00324                     if (nounit) {
00325                         temp *= a[j + j * a_dim1];
00326                     }
00327                     i__2 = *n;
00328                     for (i__ = j + 1; i__ <= i__2; ++i__) {
00329                         ix += *incx;
00330                         temp += a[i__ + j * a_dim1] * x[ix];
00331 /* L150: */
00332                     }
00333                     x[jx] = temp;
00334                     jx += *incx;
00335 /* L160: */
00336                 }
00337             }
00338         }
00339     }
00340 
00341     return 0;
00342 
00343 /*     End of DTRMV . */
00344 
00345 } /* dtrmv_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:50