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