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_SYGS2_HEADER
00038 #define TEMPLATE_LAPACK_SYGS2_HEADER
00039
00040 #include "template_lapack_common.h"
00041
00042 template<class Treal>
00043 int template_lapack_sygs2(const integer *itype, const char *uplo, const integer *n,
00044 Treal *a, const integer *lda, Treal *b, const integer *ldb, integer *
00045 info)
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 Treal c_b6 = -1.;
00117 integer c__1 = 1;
00118 Treal c_b27 = 1.;
00119
00120
00121 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
00122 Treal d__1;
00123
00124 integer k;
00125 logical upper;
00126 Treal ct;
00127 Treal akk, bkk;
00128 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00129 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00130
00131
00132 a_dim1 = *lda;
00133 a_offset = 1 + a_dim1 * 1;
00134 a -= a_offset;
00135 b_dim1 = *ldb;
00136 b_offset = 1 + b_dim1 * 1;
00137 b -= b_offset;
00138
00139
00140 *info = 0;
00141 upper = template_blas_lsame(uplo, "U");
00142 if (*itype < 1 || *itype > 3) {
00143 *info = -1;
00144 } else if (! upper && ! template_blas_lsame(uplo, "L")) {
00145 *info = -2;
00146 } else if (*n < 0) {
00147 *info = -3;
00148 } else if (*lda < maxMACRO(1,*n)) {
00149 *info = -5;
00150 } else if (*ldb < maxMACRO(1,*n)) {
00151 *info = -7;
00152 }
00153 if (*info != 0) {
00154 i__1 = -(*info);
00155 template_blas_erbla("SYGS2 ", &i__1);
00156 return 0;
00157 }
00158
00159 if (*itype == 1) {
00160 if (upper) {
00161
00162
00163
00164 i__1 = *n;
00165 for (k = 1; k <= i__1; ++k) {
00166
00167
00168
00169 akk = a_ref(k, k);
00170 bkk = b_ref(k, k);
00171
00172 d__1 = bkk;
00173 akk /= d__1 * d__1;
00174 a_ref(k, k) = akk;
00175 if (k < *n) {
00176 i__2 = *n - k;
00177 d__1 = 1. / bkk;
00178 template_blas_scal(&i__2, &d__1, &a_ref(k, k + 1), lda);
00179 ct = akk * -.5;
00180 i__2 = *n - k;
00181 template_blas_axpy(&i__2, &ct, &b_ref(k, k + 1), ldb, &a_ref(k, k + 1)
00182 , lda);
00183 i__2 = *n - k;
00184 template_blas_syr2(uplo, &i__2, &c_b6, &a_ref(k, k + 1),
00185 lda, &b_ref(k, k + 1), ldb,
00186 &a_ref(k + 1, k + 1), lda);
00187 i__2 = *n - k;
00188 template_blas_axpy(&i__2, &ct, &b_ref(k, k + 1), ldb, &a_ref(k, k + 1)
00189 , lda);
00190 i__2 = *n - k;
00191 template_blas_trsv(uplo, "Transpose", "Non-unit", &i__2, &b_ref(k + 1,
00192 k + 1), ldb, &a_ref(k, k + 1), lda);
00193 }
00194
00195 }
00196 } else {
00197
00198
00199
00200 i__1 = *n;
00201 for (k = 1; k <= i__1; ++k) {
00202
00203
00204
00205 akk = a_ref(k, k);
00206 bkk = b_ref(k, k);
00207
00208 d__1 = bkk;
00209 akk /= d__1 * d__1;
00210 a_ref(k, k) = akk;
00211 if (k < *n) {
00212 i__2 = *n - k;
00213 d__1 = 1. / bkk;
00214 template_blas_scal(&i__2, &d__1, &a_ref(k + 1, k), &c__1);
00215 ct = akk * -.5;
00216 i__2 = *n - k;
00217 template_blas_axpy(&i__2, &ct, &b_ref(k + 1, k), &c__1, &a_ref(k + 1,
00218 k), &c__1);
00219 i__2 = *n - k;
00220 template_blas_syr2(uplo, &i__2, &c_b6, &a_ref(k + 1, k),
00221 &c__1, &b_ref(k + 1, k),
00222 &c__1, &a_ref(k + 1, k + 1),
00223 lda);
00224
00225 i__2 = *n - k;
00226 template_blas_axpy(&i__2, &ct, &b_ref(k + 1, k), &c__1, &a_ref(k + 1,
00227 k), &c__1);
00228 i__2 = *n - k;
00229 template_blas_trsv(uplo, "No transpose", "Non-unit", &i__2, &b_ref(k
00230 + 1, k + 1), ldb, &a_ref(k + 1, k), &c__1);
00231 }
00232
00233 }
00234 }
00235 } else {
00236 if (upper) {
00237
00238
00239
00240 i__1 = *n;
00241 for (k = 1; k <= i__1; ++k) {
00242
00243
00244
00245 akk = a_ref(k, k);
00246 bkk = b_ref(k, k);
00247 i__2 = k - 1;
00248 template_blas_trmv(uplo, "No transpose", "Non-unit", &i__2, &b[b_offset],
00249 ldb, &a_ref(1, k), &c__1);
00250 ct = akk * .5;
00251 i__2 = k - 1;
00252 template_blas_axpy(&i__2, &ct, &b_ref(1, k), &c__1, &a_ref(1, k), &c__1);
00253 i__2 = k - 1;
00254 template_blas_syr2(uplo, &i__2, &c_b27, &a_ref(1, k), &c__1, &b_ref(1, k),
00255 &c__1, &a[a_offset], lda);
00256 i__2 = k - 1;
00257 template_blas_axpy(&i__2, &ct, &b_ref(1, k), &c__1, &a_ref(1, k), &c__1);
00258 i__2 = k - 1;
00259 template_blas_scal(&i__2, &bkk, &a_ref(1, k), &c__1);
00260
00261 d__1 = bkk;
00262 a_ref(k, k) = akk * (d__1 * d__1);
00263
00264 }
00265 } else {
00266
00267
00268
00269 i__1 = *n;
00270 for (k = 1; k <= i__1; ++k) {
00271
00272
00273
00274 akk = a_ref(k, k);
00275 bkk = b_ref(k, k);
00276 i__2 = k - 1;
00277 template_blas_trmv(uplo, "Transpose", "Non-unit", &i__2, &b[b_offset],
00278 ldb, &a_ref(k, 1), lda);
00279 ct = akk * .5;
00280 i__2 = k - 1;
00281 template_blas_axpy(&i__2, &ct, &b_ref(k, 1), ldb, &a_ref(k, 1), lda);
00282 i__2 = k - 1;
00283 template_blas_syr2(uplo, &i__2, &c_b27, &a_ref(k, 1), lda, &b_ref(k, 1),
00284 ldb, &a[a_offset], lda);
00285 i__2 = k - 1;
00286 template_blas_axpy(&i__2, &ct, &b_ref(k, 1), ldb, &a_ref(k, 1), lda);
00287 i__2 = k - 1;
00288 template_blas_scal(&i__2, &bkk, &a_ref(k, 1), lda);
00289
00290 d__1 = bkk;
00291 a_ref(k, k) = akk * (d__1 * d__1);
00292
00293 }
00294 }
00295 }
00296 return 0;
00297
00298
00299
00300 }
00301
00302 #undef b_ref
00303 #undef a_ref
00304
00305
00306 #endif