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_LAPACK_LATRD_HEADER
00038 #define TEMPLATE_LAPACK_LATRD_HEADER
00039
00040
00041 template<class Treal>
00042 int template_lapack_latrd(const char *uplo, const integer *n, const integer *nb, Treal *
00043 a, const integer *lda, Treal *e, Treal *tau, Treal *w,
00044 const integer *ldw)
00045 {
00046
00047
00048
00049
00050
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
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180 Treal c_b5 = -1.;
00181 Treal c_b6 = 1.;
00182 integer c__1 = 1;
00183 Treal c_b16 = 0.;
00184
00185
00186 integer a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3;
00187
00188 integer i__;
00189 Treal alpha;
00190 integer iw;
00191 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00192 #define w_ref(a_1,a_2) w[(a_2)*w_dim1 + a_1]
00193
00194
00195 a_dim1 = *lda;
00196 a_offset = 1 + a_dim1 * 1;
00197 a -= a_offset;
00198 --e;
00199 --tau;
00200 w_dim1 = *ldw;
00201 w_offset = 1 + w_dim1 * 1;
00202 w -= w_offset;
00203
00204
00205 if (*n <= 0) {
00206 return 0;
00207 }
00208
00209 if (template_blas_lsame(uplo, "U")) {
00210
00211
00212
00213 i__1 = *n - *nb + 1;
00214 for (i__ = *n; i__ >= i__1; --i__) {
00215 iw = i__ - *n + *nb;
00216 if (i__ < *n) {
00217
00218
00219
00220 i__2 = *n - i__;
00221 template_blas_gemv("No transpose", &i__, &i__2, &c_b5, &a_ref(1, i__ + 1),
00222 lda, &w_ref(i__, iw + 1), ldw, &c_b6, &a_ref(1, i__),
00223 &c__1);
00224 i__2 = *n - i__;
00225 template_blas_gemv("No transpose", &i__, &i__2, &c_b5, &w_ref(1, iw + 1),
00226 ldw, &a_ref(i__, i__ + 1), lda, &c_b6, &a_ref(1, i__),
00227 &c__1);
00228 }
00229 if (i__ > 1) {
00230
00231
00232
00233
00234 i__2 = i__ - 1;
00235 template_lapack_larfg(&i__2, &a_ref(i__ - 1, i__), &a_ref(1, i__), &c__1, &
00236 tau[i__ - 1]);
00237 e[i__ - 1] = a_ref(i__ - 1, i__);
00238 a_ref(i__ - 1, i__) = 1.;
00239
00240
00241
00242 i__2 = i__ - 1;
00243 template_blas_symv("Upper", &i__2, &c_b6, &a[a_offset], lda, &a_ref(1,
00244 i__), &c__1, &c_b16, &w_ref(1, iw), &c__1);
00245 if (i__ < *n) {
00246 i__2 = i__ - 1;
00247 i__3 = *n - i__;
00248 template_blas_gemv("Transpose", &i__2, &i__3, &c_b6, &w_ref(1, iw + 1)
00249 , ldw, &a_ref(1, i__), &c__1, &c_b16, &w_ref(i__
00250 + 1, iw), &c__1);
00251 i__2 = i__ - 1;
00252 i__3 = *n - i__;
00253 template_blas_gemv("No transpose", &i__2, &i__3, &c_b5, &a_ref(1, i__
00254 + 1), lda, &w_ref(i__ + 1, iw), &c__1, &c_b6, &
00255 w_ref(1, iw), &c__1);
00256 i__2 = i__ - 1;
00257 i__3 = *n - i__;
00258 template_blas_gemv("Transpose", &i__2, &i__3, &c_b6, &a_ref(1, i__ +
00259 1), lda, &a_ref(1, i__), &c__1, &c_b16, &w_ref(
00260 i__ + 1, iw), &c__1);
00261 i__2 = i__ - 1;
00262 i__3 = *n - i__;
00263 template_blas_gemv("No transpose", &i__2, &i__3, &c_b5, &w_ref(1, iw
00264 + 1), ldw, &w_ref(i__ + 1, iw), &c__1, &c_b6, &
00265 w_ref(1, iw), &c__1);
00266 }
00267 i__2 = i__ - 1;
00268 template_blas_scal(&i__2, &tau[i__ - 1], &w_ref(1, iw), &c__1);
00269 i__2 = i__ - 1;
00270 alpha = tau[i__ - 1] * -.5 * template_blas_dot(&i__2, &w_ref(1, iw), &
00271 c__1, &a_ref(1, i__), &c__1);
00272 i__2 = i__ - 1;
00273 template_blas_axpy(&i__2, &alpha, &a_ref(1, i__), &c__1, &w_ref(1, iw), &
00274 c__1);
00275 }
00276
00277
00278 }
00279 } else {
00280
00281
00282
00283 i__1 = *nb;
00284 for (i__ = 1; i__ <= i__1; ++i__) {
00285
00286
00287
00288 i__2 = *n - i__ + 1;
00289 i__3 = i__ - 1;
00290 template_blas_gemv("No transpose", &i__2, &i__3, &c_b5, &a_ref(i__, 1), lda, &
00291 w_ref(i__, 1), ldw, &c_b6, &a_ref(i__, i__), &c__1);
00292 i__2 = *n - i__ + 1;
00293 i__3 = i__ - 1;
00294 template_blas_gemv("No transpose", &i__2, &i__3, &c_b5, &w_ref(i__, 1), ldw, &
00295 a_ref(i__, 1), lda, &c_b6, &a_ref(i__, i__), &c__1);
00296 if (i__ < *n) {
00297
00298
00299
00300
00301
00302 i__2 = i__ + 2;
00303 i__3 = *n - i__;
00304 template_lapack_larfg(&i__3, &a_ref(i__ + 1, i__), &a_ref(minMACRO(i__2,*n), i__)
00305 , &c__1, &tau[i__]);
00306 e[i__] = a_ref(i__ + 1, i__);
00307 a_ref(i__ + 1, i__) = 1.;
00308
00309
00310
00311 i__2 = *n - i__;
00312 template_blas_symv("Lower", &i__2, &c_b6, &a_ref(i__ + 1, i__ + 1), lda, &
00313 a_ref(i__ + 1, i__), &c__1, &c_b16, &w_ref(i__ + 1,
00314 i__), &c__1);
00315 i__2 = *n - i__;
00316 i__3 = i__ - 1;
00317 template_blas_gemv("Transpose", &i__2, &i__3, &c_b6, &w_ref(i__ + 1, 1),
00318 ldw, &a_ref(i__ + 1, i__), &c__1, &c_b16, &w_ref(1,
00319 i__), &c__1);
00320 i__2 = *n - i__;
00321 i__3 = i__ - 1;
00322 template_blas_gemv("No transpose", &i__2, &i__3, &c_b5, &a_ref(i__ + 1, 1)
00323 , lda, &w_ref(1, i__), &c__1, &c_b6, &w_ref(i__ + 1,
00324 i__), &c__1);
00325 i__2 = *n - i__;
00326 i__3 = i__ - 1;
00327 template_blas_gemv("Transpose", &i__2, &i__3, &c_b6, &a_ref(i__ + 1, 1),
00328 lda, &a_ref(i__ + 1, i__), &c__1, &c_b16, &w_ref(1,
00329 i__), &c__1);
00330 i__2 = *n - i__;
00331 i__3 = i__ - 1;
00332 template_blas_gemv("No transpose", &i__2, &i__3, &c_b5, &w_ref(i__ + 1, 1)
00333 , ldw, &w_ref(1, i__), &c__1, &c_b6, &w_ref(i__ + 1,
00334 i__), &c__1);
00335 i__2 = *n - i__;
00336 template_blas_scal(&i__2, &tau[i__], &w_ref(i__ + 1, i__), &c__1);
00337 i__2 = *n - i__;
00338 alpha = tau[i__] * -.5 * template_blas_dot(&i__2, &w_ref(i__ + 1, i__), &
00339 c__1, &a_ref(i__ + 1, i__), &c__1);
00340 i__2 = *n - i__;
00341 template_blas_axpy(&i__2, &alpha, &a_ref(i__ + 1, i__), &c__1, &w_ref(i__
00342 + 1, i__), &c__1);
00343 }
00344
00345
00346 }
00347 }
00348
00349 return 0;
00350
00351
00352
00353 }
00354
00355 #undef w_ref
00356 #undef a_ref
00357
00358
00359 #endif