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_STEIN_HEADER
00038 #define TEMPLATE_LAPACK_STEIN_HEADER
00039
00040
00041 template<class Treal>
00042 int template_lapack_stein(const integer *n, const Treal *d__, const Treal *e,
00043 const integer *m, const Treal *w, const integer *iblock, const integer *isplit,
00044 Treal *z__, const integer *ldz, Treal *work, integer *iwork,
00045 integer *ifail, integer *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
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 integer c__2 = 2;
00144 integer c__1 = 1;
00145 integer c_n1 = -1;
00146
00147
00148 integer z_dim1, z_offset, i__1, i__2, i__3;
00149 Treal d__1, d__2, d__3, d__4, d__5;
00150
00151 integer jblk, nblk;
00152 integer jmax;
00153 integer i__, j;
00154 integer iseed[4], gpind, iinfo;
00155 integer b1;
00156 integer j1;
00157 Treal ortol;
00158 integer indrv1, indrv2, indrv3, indrv4, indrv5, bn;
00159 Treal xj;
00160 integer nrmchk;
00161 integer blksiz;
00162 Treal onenrm, dtpcrt, pertol, scl, eps, sep, nrm, tol;
00163 integer its;
00164 Treal xjm, ztr, eps1;
00165 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
00166
00167
00168 --d__;
00169 --e;
00170 --w;
00171 --iblock;
00172 --isplit;
00173 z_dim1 = *ldz;
00174 z_offset = 1 + z_dim1 * 1;
00175 z__ -= z_offset;
00176 --work;
00177 --iwork;
00178 --ifail;
00179
00180
00181 ortol = dtpcrt = xjm = onenrm = gpind = 0;
00182
00183 *info = 0;
00184 i__1 = *m;
00185 for (i__ = 1; i__ <= i__1; ++i__) {
00186 ifail[i__] = 0;
00187
00188 }
00189
00190 if (*n < 0) {
00191 *info = -1;
00192 } else if (*m < 0 || *m > *n) {
00193 *info = -4;
00194 } else if (*ldz < maxMACRO(1,*n)) {
00195 *info = -9;
00196 } else {
00197 i__1 = *m;
00198 for (j = 2; j <= i__1; ++j) {
00199 if (iblock[j] < iblock[j - 1]) {
00200 *info = -6;
00201 goto L30;
00202 }
00203 if (iblock[j] == iblock[j - 1] && w[j] < w[j - 1]) {
00204 *info = -5;
00205 goto L30;
00206 }
00207
00208 }
00209 L30:
00210 ;
00211 }
00212
00213 if (*info != 0) {
00214 i__1 = -(*info);
00215 template_blas_erbla("STEIN ", &i__1);
00216 return 0;
00217 }
00218
00219
00220
00221 if (*n == 0 || *m == 0) {
00222 return 0;
00223 } else if (*n == 1) {
00224 z___ref(1, 1) = 1.;
00225 return 0;
00226 }
00227
00228
00229
00230 eps = template_lapack_lamch("Precision", (Treal)0);
00231
00232
00233
00234 for (i__ = 1; i__ <= 4; ++i__) {
00235 iseed[i__ - 1] = 1;
00236
00237 }
00238
00239
00240
00241 indrv1 = 0;
00242 indrv2 = indrv1 + *n;
00243 indrv3 = indrv2 + *n;
00244 indrv4 = indrv3 + *n;
00245 indrv5 = indrv4 + *n;
00246
00247
00248
00249 j1 = 1;
00250 i__1 = iblock[*m];
00251 for (nblk = 1; nblk <= i__1; ++nblk) {
00252
00253
00254
00255 if (nblk == 1) {
00256 b1 = 1;
00257 } else {
00258 b1 = isplit[nblk - 1] + 1;
00259 }
00260 bn = isplit[nblk];
00261 blksiz = bn - b1 + 1;
00262 if (blksiz == 1) {
00263 goto L60;
00264 }
00265 gpind = b1;
00266
00267
00268
00269 onenrm = (d__1 = d__[b1], absMACRO(d__1)) + (d__2 = e[b1], absMACRO(d__2));
00270
00271 d__3 = onenrm, d__4 = (d__1 = d__[bn], absMACRO(d__1)) + (d__2 = e[bn - 1],
00272 absMACRO(d__2));
00273 onenrm = maxMACRO(d__3,d__4);
00274 i__2 = bn - 1;
00275 for (i__ = b1 + 1; i__ <= i__2; ++i__) {
00276
00277 d__4 = onenrm, d__5 = (d__1 = d__[i__], absMACRO(d__1)) + (d__2 = e[
00278 i__ - 1], absMACRO(d__2)) + (d__3 = e[i__], absMACRO(d__3));
00279 onenrm = maxMACRO(d__4,d__5);
00280
00281 }
00282 ortol = onenrm * .001;
00283
00284 dtpcrt = template_blas_sqrt(.1 / blksiz);
00285
00286
00287
00288 L60:
00289 jblk = 0;
00290 i__2 = *m;
00291 for (j = j1; j <= i__2; ++j) {
00292 if (iblock[j] != nblk) {
00293 j1 = j;
00294 goto L160;
00295 }
00296 ++jblk;
00297 xj = w[j];
00298
00299
00300
00301 if (blksiz == 1) {
00302 work[indrv1 + 1] = 1.;
00303 goto L120;
00304 }
00305
00306
00307
00308
00309 if (jblk > 1) {
00310 eps1 = (d__1 = eps * xj, absMACRO(d__1));
00311 pertol = eps1 * 10.;
00312 sep = xj - xjm;
00313 if (sep < pertol) {
00314 xj = xjm + pertol;
00315 }
00316 }
00317
00318 its = 0;
00319 nrmchk = 0;
00320
00321
00322
00323 template_lapack_larnv(&c__2, iseed, &blksiz, &work[indrv1 + 1]);
00324
00325
00326
00327 template_blas_copy(&blksiz, &d__[b1], &c__1, &work[indrv4 + 1], &c__1);
00328 i__3 = blksiz - 1;
00329 template_blas_copy(&i__3, &e[b1], &c__1, &work[indrv2 + 2], &c__1);
00330 i__3 = blksiz - 1;
00331 template_blas_copy(&i__3, &e[b1], &c__1, &work[indrv3 + 1], &c__1);
00332
00333
00334
00335 tol = 0.;
00336 template_lapack_lagtf(&blksiz, &work[indrv4 + 1], &xj, &work[indrv2 + 2], &work[
00337 indrv3 + 1], &tol, &work[indrv5 + 1], &iwork[1], &iinfo);
00338
00339
00340
00341 L70:
00342 ++its;
00343 if (its > 5) {
00344 goto L100;
00345 }
00346
00347
00348
00349
00350 d__2 = eps, d__3 = (d__1 = work[indrv4 + blksiz], absMACRO(d__1));
00351 scl = blksiz * onenrm * maxMACRO(d__2,d__3) / template_blas_asum(&blksiz, &work[
00352 indrv1 + 1], &c__1);
00353 template_blas_scal(&blksiz, &scl, &work[indrv1 + 1], &c__1);
00354
00355
00356
00357 template_lapack_lagts(&c_n1, &blksiz, &work[indrv4 + 1], &work[indrv2 + 2], &
00358 work[indrv3 + 1], &work[indrv5 + 1], &iwork[1], &work[
00359 indrv1 + 1], &tol, &iinfo);
00360
00361
00362
00363
00364 if (jblk == 1) {
00365 goto L90;
00366 }
00367 if ((d__1 = xj - xjm, absMACRO(d__1)) > ortol) {
00368 gpind = j;
00369 }
00370 if (gpind != j) {
00371 i__3 = j - 1;
00372 for (i__ = gpind; i__ <= i__3; ++i__) {
00373 ztr = -template_blas_dot(&blksiz, &work[indrv1 + 1], &c__1, &z___ref(
00374 b1, i__), &c__1);
00375 template_blas_axpy(&blksiz, &ztr, &z___ref(b1, i__), &c__1, &work[
00376 indrv1 + 1], &c__1);
00377
00378 }
00379 }
00380
00381
00382
00383 L90:
00384 jmax = template_blas_idamax(&blksiz, &work[indrv1 + 1], &c__1);
00385 nrm = (d__1 = work[indrv1 + jmax], absMACRO(d__1));
00386
00387
00388
00389
00390 if (nrm < dtpcrt) {
00391 goto L70;
00392 }
00393 ++nrmchk;
00394 if (nrmchk < 3) {
00395 goto L70;
00396 }
00397
00398 goto L110;
00399
00400
00401
00402
00403 L100:
00404 ++(*info);
00405 ifail[*info] = j;
00406
00407
00408
00409 L110:
00410 scl = 1. / template_blas_nrm2(&blksiz, &work[indrv1 + 1], &c__1);
00411 jmax = template_blas_idamax(&blksiz, &work[indrv1 + 1], &c__1);
00412 if (work[indrv1 + jmax] < 0.) {
00413 scl = -scl;
00414 }
00415 template_blas_scal(&blksiz, &scl, &work[indrv1 + 1], &c__1);
00416 L120:
00417 i__3 = *n;
00418 for (i__ = 1; i__ <= i__3; ++i__) {
00419 z___ref(i__, j) = 0.;
00420
00421 }
00422 i__3 = blksiz;
00423 for (i__ = 1; i__ <= i__3; ++i__) {
00424 z___ref(b1 + i__ - 1, j) = work[indrv1 + i__];
00425
00426 }
00427
00428
00429
00430
00431 xjm = xj;
00432
00433
00434 }
00435 L160:
00436 ;
00437 }
00438
00439 return 0;
00440
00441
00442
00443 }
00444
00445 #undef z___ref
00446
00447
00448 #endif