00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #ifndef TEMPLATE_BLAS_TPMV_HEADER
00038 #define TEMPLATE_BLAS_TPMV_HEADER
00039
00040
00041 template<class Treal>
00042 int template_blas_tpmv(const char *uplo, const char *trans, const char *diag, const integer *n,
00043 const Treal *ap, Treal *x, const integer *incx)
00044 {
00045
00046 integer i__1, i__2;
00047
00048 integer info;
00049 Treal temp;
00050 integer i__, j, k;
00051 integer kk, ix, jx, kx;
00052 logical nounit;
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117 --x;
00118 --ap;
00119
00120 kx = 0;
00121
00122 info = 0;
00123 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
00124 info = 1;
00125 } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans,
00126 "T") && ! template_blas_lsame(trans, "C")) {
00127 info = 2;
00128 } else if (! template_blas_lsame(diag, "U") && ! template_blas_lsame(diag,
00129 "N")) {
00130 info = 3;
00131 } else if (*n < 0) {
00132 info = 4;
00133 } else if (*incx == 0) {
00134 info = 7;
00135 }
00136 if (info != 0) {
00137 template_blas_erbla("TPMV ", &info);
00138 return 0;
00139 }
00140
00141 if (*n == 0) {
00142 return 0;
00143 }
00144 nounit = template_blas_lsame(diag, "N");
00145
00146
00147 if (*incx <= 0) {
00148 kx = 1 - (*n - 1) * *incx;
00149 } else if (*incx != 1) {
00150 kx = 1;
00151 }
00152
00153
00154 if (template_blas_lsame(trans, "N")) {
00155
00156 if (template_blas_lsame(uplo, "U")) {
00157 kk = 1;
00158 if (*incx == 1) {
00159 i__1 = *n;
00160 for (j = 1; j <= i__1; ++j) {
00161 if (x[j] != 0.) {
00162 temp = x[j];
00163 k = kk;
00164 i__2 = j - 1;
00165 for (i__ = 1; i__ <= i__2; ++i__) {
00166 x[i__] += temp * ap[k];
00167 ++k;
00168
00169 }
00170 if (nounit) {
00171 x[j] *= ap[kk + j - 1];
00172 }
00173 }
00174 kk += j;
00175
00176 }
00177 } else {
00178 jx = kx;
00179 i__1 = *n;
00180 for (j = 1; j <= i__1; ++j) {
00181 if (x[jx] != 0.) {
00182 temp = x[jx];
00183 ix = kx;
00184 i__2 = kk + j - 2;
00185 for (k = kk; k <= i__2; ++k) {
00186 x[ix] += temp * ap[k];
00187 ix += *incx;
00188
00189 }
00190 if (nounit) {
00191 x[jx] *= ap[kk + j - 1];
00192 }
00193 }
00194 jx += *incx;
00195 kk += j;
00196
00197 }
00198 }
00199 } else {
00200 kk = *n * (*n + 1) / 2;
00201 if (*incx == 1) {
00202 for (j = *n; j >= 1; --j) {
00203 if (x[j] != 0.) {
00204 temp = x[j];
00205 k = kk;
00206 i__1 = j + 1;
00207 for (i__ = *n; i__ >= i__1; --i__) {
00208 x[i__] += temp * ap[k];
00209 --k;
00210
00211 }
00212 if (nounit) {
00213 x[j] *= ap[kk - *n + j];
00214 }
00215 }
00216 kk -= *n - j + 1;
00217
00218 }
00219 } else {
00220 kx += (*n - 1) * *incx;
00221 jx = kx;
00222 for (j = *n; j >= 1; --j) {
00223 if (x[jx] != 0.) {
00224 temp = x[jx];
00225 ix = kx;
00226 i__1 = kk - (*n - (j + 1));
00227 for (k = kk; k >= i__1; --k) {
00228 x[ix] += temp * ap[k];
00229 ix -= *incx;
00230
00231 }
00232 if (nounit) {
00233 x[jx] *= ap[kk - *n + j];
00234 }
00235 }
00236 jx -= *incx;
00237 kk -= *n - j + 1;
00238
00239 }
00240 }
00241 }
00242 } else {
00243
00244 if (template_blas_lsame(uplo, "U")) {
00245 kk = *n * (*n + 1) / 2;
00246 if (*incx == 1) {
00247 for (j = *n; j >= 1; --j) {
00248 temp = x[j];
00249 if (nounit) {
00250 temp *= ap[kk];
00251 }
00252 k = kk - 1;
00253 for (i__ = j - 1; i__ >= 1; --i__) {
00254 temp += ap[k] * x[i__];
00255 --k;
00256
00257 }
00258 x[j] = temp;
00259 kk -= j;
00260
00261 }
00262 } else {
00263 jx = kx + (*n - 1) * *incx;
00264 for (j = *n; j >= 1; --j) {
00265 temp = x[jx];
00266 ix = jx;
00267 if (nounit) {
00268 temp *= ap[kk];
00269 }
00270 i__1 = kk - j + 1;
00271 for (k = kk - 1; k >= i__1; --k) {
00272 ix -= *incx;
00273 temp += ap[k] * x[ix];
00274
00275 }
00276 x[jx] = temp;
00277 jx -= *incx;
00278 kk -= j;
00279
00280 }
00281 }
00282 } else {
00283 kk = 1;
00284 if (*incx == 1) {
00285 i__1 = *n;
00286 for (j = 1; j <= i__1; ++j) {
00287 temp = x[j];
00288 if (nounit) {
00289 temp *= ap[kk];
00290 }
00291 k = kk + 1;
00292 i__2 = *n;
00293 for (i__ = j + 1; i__ <= i__2; ++i__) {
00294 temp += ap[k] * x[i__];
00295 ++k;
00296
00297 }
00298 x[j] = temp;
00299 kk += *n - j + 1;
00300
00301 }
00302 } else {
00303 jx = kx;
00304 i__1 = *n;
00305 for (j = 1; j <= i__1; ++j) {
00306 temp = x[jx];
00307 ix = jx;
00308 if (nounit) {
00309 temp *= ap[kk];
00310 }
00311 i__2 = kk + *n - j;
00312 for (k = kk + 1; k <= i__2; ++k) {
00313 ix += *incx;
00314 temp += ap[k] * x[ix];
00315
00316 }
00317 x[jx] = temp;
00318 jx += *incx;
00319 kk += *n - j + 1;
00320
00321 }
00322 }
00323 }
00324 }
00325 return 0;
00326
00327 }
00328
00329 #endif