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_TPSV_HEADER
00038 #define TEMPLATE_BLAS_TPSV_HEADER
00039
00040
00041 template<class Treal>
00042 int template_blas_tpsv(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
00118
00119 --x;
00120 --ap;
00121
00122 kx = 0;
00123
00124 info = 0;
00125 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
00126 info = 1;
00127 } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans,
00128 "T") && ! template_blas_lsame(trans, "C")) {
00129 info = 2;
00130 } else if (! template_blas_lsame(diag, "U") && ! template_blas_lsame(diag,
00131 "N")) {
00132 info = 3;
00133 } else if (*n < 0) {
00134 info = 4;
00135 } else if (*incx == 0) {
00136 info = 7;
00137 }
00138 if (info != 0) {
00139 template_blas_erbla("DTPSV ", &info);
00140 return 0;
00141 }
00142
00143 if (*n == 0) {
00144 return 0;
00145 }
00146 nounit = template_blas_lsame(diag, "N");
00147
00148
00149 if (*incx <= 0) {
00150 kx = 1 - (*n - 1) * *incx;
00151 } else if (*incx != 1) {
00152 kx = 1;
00153 }
00154
00155
00156 if (template_blas_lsame(trans, "N")) {
00157
00158 if (template_blas_lsame(uplo, "U")) {
00159 kk = *n * (*n + 1) / 2;
00160 if (*incx == 1) {
00161 for (j = *n; j >= 1; --j) {
00162 if (x[j] != 0.) {
00163 if (nounit) {
00164 x[j] /= ap[kk];
00165 }
00166 temp = x[j];
00167 k = kk - 1;
00168 for (i__ = j - 1; i__ >= 1; --i__) {
00169 x[i__] -= temp * ap[k];
00170 --k;
00171
00172 }
00173 }
00174 kk -= j;
00175
00176 }
00177 } else {
00178 jx = kx + (*n - 1) * *incx;
00179 for (j = *n; j >= 1; --j) {
00180 if (x[jx] != 0.) {
00181 if (nounit) {
00182 x[jx] /= ap[kk];
00183 }
00184 temp = x[jx];
00185 ix = jx;
00186 i__1 = kk - j + 1;
00187 for (k = kk - 1; k >= i__1; --k) {
00188 ix -= *incx;
00189 x[ix] -= temp * ap[k];
00190
00191 }
00192 }
00193 jx -= *incx;
00194 kk -= j;
00195
00196 }
00197 }
00198 } else {
00199 kk = 1;
00200 if (*incx == 1) {
00201 i__1 = *n;
00202 for (j = 1; j <= i__1; ++j) {
00203 if (x[j] != 0.) {
00204 if (nounit) {
00205 x[j] /= ap[kk];
00206 }
00207 temp = x[j];
00208 k = kk + 1;
00209 i__2 = *n;
00210 for (i__ = j + 1; i__ <= i__2; ++i__) {
00211 x[i__] -= temp * ap[k];
00212 ++k;
00213
00214 }
00215 }
00216 kk += *n - j + 1;
00217
00218 }
00219 } else {
00220 jx = kx;
00221 i__1 = *n;
00222 for (j = 1; j <= i__1; ++j) {
00223 if (x[jx] != 0.) {
00224 if (nounit) {
00225 x[jx] /= ap[kk];
00226 }
00227 temp = x[jx];
00228 ix = jx;
00229 i__2 = kk + *n - j;
00230 for (k = kk + 1; k <= i__2; ++k) {
00231 ix += *incx;
00232 x[ix] -= temp * ap[k];
00233
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 = 1;
00246 if (*incx == 1) {
00247 i__1 = *n;
00248 for (j = 1; j <= i__1; ++j) {
00249 temp = x[j];
00250 k = kk;
00251 i__2 = j - 1;
00252 for (i__ = 1; i__ <= i__2; ++i__) {
00253 temp -= ap[k] * x[i__];
00254 ++k;
00255
00256 }
00257 if (nounit) {
00258 temp /= ap[kk + j - 1];
00259 }
00260 x[j] = temp;
00261 kk += j;
00262
00263 }
00264 } else {
00265 jx = kx;
00266 i__1 = *n;
00267 for (j = 1; j <= i__1; ++j) {
00268 temp = x[jx];
00269 ix = kx;
00270 i__2 = kk + j - 2;
00271 for (k = kk; k <= i__2; ++k) {
00272 temp -= ap[k] * x[ix];
00273 ix += *incx;
00274
00275 }
00276 if (nounit) {
00277 temp /= ap[kk + j - 1];
00278 }
00279 x[jx] = temp;
00280 jx += *incx;
00281 kk += j;
00282
00283 }
00284 }
00285 } else {
00286 kk = *n * (*n + 1) / 2;
00287 if (*incx == 1) {
00288 for (j = *n; j >= 1; --j) {
00289 temp = x[j];
00290 k = kk;
00291 i__1 = j + 1;
00292 for (i__ = *n; i__ >= i__1; --i__) {
00293 temp -= ap[k] * x[i__];
00294 --k;
00295
00296 }
00297 if (nounit) {
00298 temp /= ap[kk - *n + j];
00299 }
00300 x[j] = temp;
00301 kk -= *n - j + 1;
00302
00303 }
00304 } else {
00305 kx += (*n - 1) * *incx;
00306 jx = kx;
00307 for (j = *n; j >= 1; --j) {
00308 temp = x[jx];
00309 ix = kx;
00310 i__1 = kk - (*n - (j + 1));
00311 for (k = kk; k >= i__1; --k) {
00312 temp -= ap[k] * x[ix];
00313 ix -= *incx;
00314
00315 }
00316 if (nounit) {
00317 temp /= ap[kk - *n + j];
00318 }
00319 x[jx] = temp;
00320 jx -= *incx;
00321 kk -= *n - j + 1;
00322
00323 }
00324 }
00325 }
00326 }
00327 return 0;
00328
00329 }
00330
00331 #endif