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_SPGST_HEADER
00038 #define TEMPLATE_LAPACK_SPGST_HEADER
00039
00040 #include "template_lapack_common.h"
00041
00042 template<class Treal>
00043 int template_lapack_spgst(const integer *itype, const char *uplo, const integer *n,
00044 Treal *ap, const Treal *bp, integer *info)
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 integer c__1 = 1;
00108 Treal c_b9 = -1.;
00109 Treal c_b11 = 1.;
00110
00111
00112 integer i__1, i__2;
00113 Treal d__1;
00114
00115 integer j, k;
00116 logical upper;
00117 integer j1, k1;
00118 integer jj, kk;
00119 Treal ct;
00120 Treal ajj;
00121 integer j1j1;
00122 Treal akk;
00123 integer k1k1;
00124 Treal bjj, bkk;
00125
00126
00127 --bp;
00128 --ap;
00129
00130
00131 *info = 0;
00132 upper = template_blas_lsame(uplo, "U");
00133 if (*itype < 1 || *itype > 3) {
00134 *info = -1;
00135 } else if (! upper && ! template_blas_lsame(uplo, "L")) {
00136 *info = -2;
00137 } else if (*n < 0) {
00138 *info = -3;
00139 }
00140 if (*info != 0) {
00141 i__1 = -(*info);
00142 template_blas_erbla("SPGST ", &i__1);
00143 return 0;
00144 }
00145
00146 if (*itype == 1) {
00147 if (upper) {
00148
00149
00150
00151
00152
00153 jj = 0;
00154 i__1 = *n;
00155 for (j = 1; j <= i__1; ++j) {
00156 j1 = jj + 1;
00157 jj += j;
00158
00159
00160
00161 bjj = bp[jj];
00162 template_blas_tpsv(uplo, "Transpose", "Nonunit", &j, &bp[1], &ap[j1], &
00163 c__1);
00164 i__2 = j - 1;
00165 template_blas_spmv(uplo, &i__2, &c_b9, &ap[1], &bp[j1], &c__1, &c_b11, &
00166 ap[j1], &c__1);
00167 i__2 = j - 1;
00168 d__1 = 1. / bjj;
00169 template_blas_scal(&i__2, &d__1, &ap[j1], &c__1);
00170 i__2 = j - 1;
00171 ap[jj] = (ap[jj] - template_blas_dot(&i__2, &ap[j1], &c__1, &bp[j1], &
00172 c__1)) / bjj;
00173
00174 }
00175 } else {
00176
00177
00178
00179
00180
00181 kk = 1;
00182 i__1 = *n;
00183 for (k = 1; k <= i__1; ++k) {
00184 k1k1 = kk + *n - k + 1;
00185
00186
00187
00188 akk = ap[kk];
00189 bkk = bp[kk];
00190
00191 d__1 = bkk;
00192 akk /= d__1 * d__1;
00193 ap[kk] = akk;
00194 if (k < *n) {
00195 i__2 = *n - k;
00196 d__1 = 1. / bkk;
00197 template_blas_scal(&i__2, &d__1, &ap[kk + 1], &c__1);
00198 ct = akk * -.5;
00199 i__2 = *n - k;
00200 template_blas_axpy(&i__2, &ct, &bp[kk + 1], &c__1, &ap[kk + 1], &c__1)
00201 ;
00202 i__2 = *n - k;
00203 template_blas_spr2(uplo, &i__2, &c_b9, &ap[kk + 1], &c__1, &bp[kk + 1]
00204 , &c__1, &ap[k1k1]);
00205 i__2 = *n - k;
00206 template_blas_axpy(&i__2, &ct, &bp[kk + 1], &c__1, &ap[kk + 1], &c__1)
00207 ;
00208 i__2 = *n - k;
00209 template_blas_tpsv(uplo, "No transpose", "Non-unit", &i__2, &bp[k1k1],
00210 &ap[kk + 1], &c__1);
00211 }
00212 kk = k1k1;
00213
00214 }
00215 }
00216 } else {
00217 if (upper) {
00218
00219
00220
00221
00222
00223 kk = 0;
00224 i__1 = *n;
00225 for (k = 1; k <= i__1; ++k) {
00226 k1 = kk + 1;
00227 kk += k;
00228
00229
00230
00231 akk = ap[kk];
00232 bkk = bp[kk];
00233 i__2 = k - 1;
00234 template_blas_tpmv(uplo, "No transpose", "Non-unit", &i__2, &bp[1], &ap[
00235 k1], &c__1);
00236 ct = akk * .5;
00237 i__2 = k - 1;
00238 template_blas_axpy(&i__2, &ct, &bp[k1], &c__1, &ap[k1], &c__1);
00239 i__2 = k - 1;
00240 template_blas_spr2(uplo, &i__2, &c_b11, &ap[k1], &c__1, &bp[k1], &c__1, &
00241 ap[1]);
00242 i__2 = k - 1;
00243 template_blas_axpy(&i__2, &ct, &bp[k1], &c__1, &ap[k1], &c__1);
00244 i__2 = k - 1;
00245 template_blas_scal(&i__2, &bkk, &ap[k1], &c__1);
00246
00247 d__1 = bkk;
00248 ap[kk] = akk * (d__1 * d__1);
00249
00250 }
00251 } else {
00252
00253
00254
00255
00256
00257 jj = 1;
00258 i__1 = *n;
00259 for (j = 1; j <= i__1; ++j) {
00260 j1j1 = jj + *n - j + 1;
00261
00262
00263
00264 ajj = ap[jj];
00265 bjj = bp[jj];
00266 i__2 = *n - j;
00267 ap[jj] = ajj * bjj + template_blas_dot(&i__2, &ap[jj + 1], &c__1, &bp[jj
00268 + 1], &c__1);
00269 i__2 = *n - j;
00270 template_blas_scal(&i__2, &bjj, &ap[jj + 1], &c__1);
00271 i__2 = *n - j;
00272 template_blas_spmv(uplo, &i__2, &c_b11, &ap[j1j1], &bp[jj + 1], &c__1, &
00273 c_b11, &ap[jj + 1], &c__1);
00274 i__2 = *n - j + 1;
00275 template_blas_tpmv(uplo, "Transpose", "Non-unit", &i__2, &bp[jj], &ap[jj],
00276 &c__1);
00277 jj = j1j1;
00278
00279 }
00280 }
00281 }
00282 return 0;
00283
00284
00285
00286 }
00287
00288 #endif