template_blas_tpsv.h
Go to the documentation of this file.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 #ifndef TEMPLATE_BLAS_TPSV_HEADER
00036 #define TEMPLATE_BLAS_TPSV_HEADER
00037
00038
00039 template<class Treal>
00040 int template_blas_tpsv(const char *uplo, const char *trans, const char *diag, const integer *n,
00041 const Treal *ap, Treal *x, const integer *incx)
00042 {
00043
00044 integer i__1, i__2;
00045
00046 integer info;
00047 Treal temp;
00048 integer i__, j, k;
00049 integer kk, ix, jx, kx;
00050 logical nounit;
00051
00052
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("DTPSV ", &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 = *n * (*n + 1) / 2;
00158 if (*incx == 1) {
00159 for (j = *n; j >= 1; --j) {
00160 if (x[j] != 0.) {
00161 if (nounit) {
00162 x[j] /= ap[kk];
00163 }
00164 temp = x[j];
00165 k = kk - 1;
00166 for (i__ = j - 1; i__ >= 1; --i__) {
00167 x[i__] -= temp * ap[k];
00168 --k;
00169
00170 }
00171 }
00172 kk -= j;
00173
00174 }
00175 } else {
00176 jx = kx + (*n - 1) * *incx;
00177 for (j = *n; j >= 1; --j) {
00178 if (x[jx] != 0.) {
00179 if (nounit) {
00180 x[jx] /= ap[kk];
00181 }
00182 temp = x[jx];
00183 ix = jx;
00184 i__1 = kk - j + 1;
00185 for (k = kk - 1; k >= i__1; --k) {
00186 ix -= *incx;
00187 x[ix] -= temp * ap[k];
00188
00189 }
00190 }
00191 jx -= *incx;
00192 kk -= j;
00193
00194 }
00195 }
00196 } else {
00197 kk = 1;
00198 if (*incx == 1) {
00199 i__1 = *n;
00200 for (j = 1; j <= i__1; ++j) {
00201 if (x[j] != 0.) {
00202 if (nounit) {
00203 x[j] /= ap[kk];
00204 }
00205 temp = x[j];
00206 k = kk + 1;
00207 i__2 = *n;
00208 for (i__ = j + 1; i__ <= i__2; ++i__) {
00209 x[i__] -= temp * ap[k];
00210 ++k;
00211
00212 }
00213 }
00214 kk += *n - j + 1;
00215
00216 }
00217 } else {
00218 jx = kx;
00219 i__1 = *n;
00220 for (j = 1; j <= i__1; ++j) {
00221 if (x[jx] != 0.) {
00222 if (nounit) {
00223 x[jx] /= ap[kk];
00224 }
00225 temp = x[jx];
00226 ix = jx;
00227 i__2 = kk + *n - j;
00228 for (k = kk + 1; k <= i__2; ++k) {
00229 ix += *incx;
00230 x[ix] -= temp * ap[k];
00231
00232 }
00233 }
00234 jx += *incx;
00235 kk += *n - j + 1;
00236
00237 }
00238 }
00239 }
00240 } else {
00241
00242 if (template_blas_lsame(uplo, "U")) {
00243 kk = 1;
00244 if (*incx == 1) {
00245 i__1 = *n;
00246 for (j = 1; j <= i__1; ++j) {
00247 temp = x[j];
00248 k = kk;
00249 i__2 = j - 1;
00250 for (i__ = 1; i__ <= i__2; ++i__) {
00251 temp -= ap[k] * x[i__];
00252 ++k;
00253
00254 }
00255 if (nounit) {
00256 temp /= ap[kk + j - 1];
00257 }
00258 x[j] = temp;
00259 kk += j;
00260
00261 }
00262 } else {
00263 jx = kx;
00264 i__1 = *n;
00265 for (j = 1; j <= i__1; ++j) {
00266 temp = x[jx];
00267 ix = kx;
00268 i__2 = kk + j - 2;
00269 for (k = kk; k <= i__2; ++k) {
00270 temp -= ap[k] * x[ix];
00271 ix += *incx;
00272
00273 }
00274 if (nounit) {
00275 temp /= ap[kk + j - 1];
00276 }
00277 x[jx] = temp;
00278 jx += *incx;
00279 kk += j;
00280
00281 }
00282 }
00283 } else {
00284 kk = *n * (*n + 1) / 2;
00285 if (*incx == 1) {
00286 for (j = *n; j >= 1; --j) {
00287 temp = x[j];
00288 k = kk;
00289 i__1 = j + 1;
00290 for (i__ = *n; i__ >= i__1; --i__) {
00291 temp -= ap[k] * x[i__];
00292 --k;
00293
00294 }
00295 if (nounit) {
00296 temp /= ap[kk - *n + j];
00297 }
00298 x[j] = temp;
00299 kk -= *n - j + 1;
00300
00301 }
00302 } else {
00303 kx += (*n - 1) * *incx;
00304 jx = kx;
00305 for (j = *n; j >= 1; --j) {
00306 temp = x[jx];
00307 ix = kx;
00308 i__1 = kk - (*n - (j + 1));
00309 for (k = kk; k >= i__1; --k) {
00310 temp -= ap[k] * x[ix];
00311 ix -= *incx;
00312
00313 }
00314 if (nounit) {
00315 temp /= ap[kk - *n + j];
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