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_SYMV_HEADER
00036 #define TEMPLATE_BLAS_SYMV_HEADER
00037
00038
00039 template<class Treal>
00040 int template_blas_symv(const char *uplo, const integer *n, const Treal *alpha,
00041 const Treal *a, const integer *lda, const Treal *x, const integer *incx, const Treal
00042 *beta, Treal *y, const integer *incy)
00043 {
00044
00045 integer a_dim1, a_offset, i__1, i__2;
00046
00047 integer info;
00048 Treal temp1, temp2;
00049 integer i__, j;
00050 integer ix, iy, jx, jy, kx, ky;
00051 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
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
00118
00119
00120
00121 a_dim1 = *lda;
00122 a_offset = 1 + a_dim1 * 1;
00123 a -= a_offset;
00124 --x;
00125 --y;
00126
00127 info = 0;
00128 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
00129 info = 1;
00130 } else if (*n < 0) {
00131 info = 2;
00132 } else if (*lda < maxMACRO(1,*n)) {
00133 info = 5;
00134 } else if (*incx == 0) {
00135 info = 7;
00136 } else if (*incy == 0) {
00137 info = 10;
00138 }
00139 if (info != 0) {
00140 template_blas_erbla("SYMV ", &info);
00141 return 0;
00142 }
00143
00144 if (*n == 0 || (*alpha == 0. && *beta == 1.) ) {
00145 return 0;
00146 }
00147
00148 if (*incx > 0) {
00149 kx = 1;
00150 } else {
00151 kx = 1 - (*n - 1) * *incx;
00152 }
00153 if (*incy > 0) {
00154 ky = 1;
00155 } else {
00156 ky = 1 - (*n - 1) * *incy;
00157 }
00158
00159
00160
00161
00162 if (*beta != 1.) {
00163 if (*incy == 1) {
00164 if (*beta == 0.) {
00165 i__1 = *n;
00166 for (i__ = 1; i__ <= i__1; ++i__) {
00167 y[i__] = 0.;
00168
00169 }
00170 } else {
00171 i__1 = *n;
00172 for (i__ = 1; i__ <= i__1; ++i__) {
00173 y[i__] = *beta * y[i__];
00174
00175 }
00176 }
00177 } else {
00178 iy = ky;
00179 if (*beta == 0.) {
00180 i__1 = *n;
00181 for (i__ = 1; i__ <= i__1; ++i__) {
00182 y[iy] = 0.;
00183 iy += *incy;
00184
00185 }
00186 } else {
00187 i__1 = *n;
00188 for (i__ = 1; i__ <= i__1; ++i__) {
00189 y[iy] = *beta * y[iy];
00190 iy += *incy;
00191
00192 }
00193 }
00194 }
00195 }
00196 if (*alpha == 0.) {
00197 return 0;
00198 }
00199 if (template_blas_lsame(uplo, "U")) {
00200
00201 if (*incx == 1 && *incy == 1) {
00202 i__1 = *n;
00203 for (j = 1; j <= i__1; ++j) {
00204 temp1 = *alpha * x[j];
00205 temp2 = 0.;
00206 i__2 = j - 1;
00207 for (i__ = 1; i__ <= i__2; ++i__) {
00208 y[i__] += temp1 * a_ref(i__, j);
00209 temp2 += a_ref(i__, j) * x[i__];
00210
00211 }
00212 y[j] = y[j] + temp1 * a_ref(j, j) + *alpha * temp2;
00213
00214 }
00215 } else {
00216 jx = kx;
00217 jy = ky;
00218 i__1 = *n;
00219 for (j = 1; j <= i__1; ++j) {
00220 temp1 = *alpha * x[jx];
00221 temp2 = 0.;
00222 ix = kx;
00223 iy = ky;
00224 i__2 = j - 1;
00225 for (i__ = 1; i__ <= i__2; ++i__) {
00226 y[iy] += temp1 * a_ref(i__, j);
00227 temp2 += a_ref(i__, j) * x[ix];
00228 ix += *incx;
00229 iy += *incy;
00230
00231 }
00232 y[jy] = y[jy] + temp1 * a_ref(j, j) + *alpha * temp2;
00233 jx += *incx;
00234 jy += *incy;
00235
00236 }
00237 }
00238 } else {
00239
00240 if (*incx == 1 && *incy == 1) {
00241 i__1 = *n;
00242 for (j = 1; j <= i__1; ++j) {
00243 temp1 = *alpha * x[j];
00244 temp2 = 0.;
00245 y[j] += temp1 * a_ref(j, j);
00246 i__2 = *n;
00247 for (i__ = j + 1; i__ <= i__2; ++i__) {
00248 y[i__] += temp1 * a_ref(i__, j);
00249 temp2 += a_ref(i__, j) * x[i__];
00250
00251 }
00252 y[j] += *alpha * temp2;
00253
00254 }
00255 } else {
00256 jx = kx;
00257 jy = ky;
00258 i__1 = *n;
00259 for (j = 1; j <= i__1; ++j) {
00260 temp1 = *alpha * x[jx];
00261 temp2 = 0.;
00262 y[jy] += temp1 * a_ref(j, j);
00263 ix = jx;
00264 iy = jy;
00265 i__2 = *n;
00266 for (i__ = j + 1; i__ <= i__2; ++i__) {
00267 ix += *incx;
00268 iy += *incy;
00269 y[iy] += temp1 * a_ref(i__, j);
00270 temp2 += a_ref(i__, j) * x[ix];
00271
00272 }
00273 y[jy] += *alpha * temp2;
00274 jx += *incx;
00275 jy += *incy;
00276
00277 }
00278 }
00279 }
00280 return 0;
00281
00282 }
00283 #undef a_ref
00284
00285 #endif