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