00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
00017
00018 static integer c__1 = 1;
00019
00020 int ztrsna_(char *job, char *howmny, logical *select,
00021 integer *n, doublecomplex *t, integer *ldt, doublecomplex *vl,
00022 integer *ldvl, doublecomplex *vr, integer *ldvr, doublereal *s,
00023 doublereal *sep, integer *mm, integer *m, doublecomplex *work,
00024 integer *ldwork, doublereal *rwork, integer *info)
00025 {
00026
00027 integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset,
00028 work_dim1, work_offset, i__1, i__2, i__3, i__4, i__5;
00029 doublereal d__1, d__2;
00030 doublecomplex z__1;
00031
00032
00033 double z_abs(doublecomplex *), d_imag(doublecomplex *);
00034
00035
00036 integer i__, j, k, ks, ix;
00037 doublereal eps, est;
00038 integer kase, ierr;
00039 doublecomplex prod;
00040 doublereal lnrm, rnrm, scale;
00041 extern logical lsame_(char *, char *);
00042 integer isave[3];
00043 extern VOID zdotc_(doublecomplex *, integer *,
00044 doublecomplex *, integer *, doublecomplex *, integer *);
00045 doublecomplex dummy[1];
00046 logical wants;
00047 doublereal xnorm;
00048 extern int zlacn2_(integer *, doublecomplex *,
00049 doublecomplex *, doublereal *, integer *, integer *), dlabad_(
00050 doublereal *, doublereal *);
00051 extern doublereal dznrm2_(integer *, doublecomplex *, integer *), dlamch_(
00052 char *);
00053 extern int xerbla_(char *, integer *);
00054 doublereal bignum;
00055 logical wantbh;
00056 extern integer izamax_(integer *, doublecomplex *, integer *);
00057 logical somcon;
00058 extern int zdrscl_(integer *, doublereal *,
00059 doublecomplex *, integer *);
00060 char normin[1];
00061 extern int zlacpy_(char *, integer *, integer *,
00062 doublecomplex *, integer *, doublecomplex *, integer *);
00063 doublereal smlnum;
00064 logical wantsp;
00065 extern int zlatrs_(char *, char *, char *, char *,
00066 integer *, doublecomplex *, integer *, doublecomplex *,
00067 doublereal *, doublereal *, integer *), ztrexc_(char *, integer *, doublecomplex *, integer *,
00068 doublecomplex *, integer *, integer *, integer *, integer *);
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
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 --select;
00244 t_dim1 = *ldt;
00245 t_offset = 1 + t_dim1;
00246 t -= t_offset;
00247 vl_dim1 = *ldvl;
00248 vl_offset = 1 + vl_dim1;
00249 vl -= vl_offset;
00250 vr_dim1 = *ldvr;
00251 vr_offset = 1 + vr_dim1;
00252 vr -= vr_offset;
00253 --s;
00254 --sep;
00255 work_dim1 = *ldwork;
00256 work_offset = 1 + work_dim1;
00257 work -= work_offset;
00258 --rwork;
00259
00260
00261 wantbh = lsame_(job, "B");
00262 wants = lsame_(job, "E") || wantbh;
00263 wantsp = lsame_(job, "V") || wantbh;
00264
00265 somcon = lsame_(howmny, "S");
00266
00267
00268
00269
00270 if (somcon) {
00271 *m = 0;
00272 i__1 = *n;
00273 for (j = 1; j <= i__1; ++j) {
00274 if (select[j]) {
00275 ++(*m);
00276 }
00277
00278 }
00279 } else {
00280 *m = *n;
00281 }
00282
00283 *info = 0;
00284 if (! wants && ! wantsp) {
00285 *info = -1;
00286 } else if (! lsame_(howmny, "A") && ! somcon) {
00287 *info = -2;
00288 } else if (*n < 0) {
00289 *info = -4;
00290 } else if (*ldt < max(1,*n)) {
00291 *info = -6;
00292 } else if (*ldvl < 1 || wants && *ldvl < *n) {
00293 *info = -8;
00294 } else if (*ldvr < 1 || wants && *ldvr < *n) {
00295 *info = -10;
00296 } else if (*mm < *m) {
00297 *info = -13;
00298 } else if (*ldwork < 1 || wantsp && *ldwork < *n) {
00299 *info = -16;
00300 }
00301 if (*info != 0) {
00302 i__1 = -(*info);
00303 xerbla_("ZTRSNA", &i__1);
00304 return 0;
00305 }
00306
00307
00308
00309 if (*n == 0) {
00310 return 0;
00311 }
00312
00313 if (*n == 1) {
00314 if (somcon) {
00315 if (! select[1]) {
00316 return 0;
00317 }
00318 }
00319 if (wants) {
00320 s[1] = 1.;
00321 }
00322 if (wantsp) {
00323 sep[1] = z_abs(&t[t_dim1 + 1]);
00324 }
00325 return 0;
00326 }
00327
00328
00329
00330 eps = dlamch_("P");
00331 smlnum = dlamch_("S") / eps;
00332 bignum = 1. / smlnum;
00333 dlabad_(&smlnum, &bignum);
00334
00335 ks = 1;
00336 i__1 = *n;
00337 for (k = 1; k <= i__1; ++k) {
00338
00339 if (somcon) {
00340 if (! select[k]) {
00341 goto L50;
00342 }
00343 }
00344
00345 if (wants) {
00346
00347
00348
00349
00350 zdotc_(&z__1, n, &vr[ks * vr_dim1 + 1], &c__1, &vl[ks * vl_dim1 +
00351 1], &c__1);
00352 prod.r = z__1.r, prod.i = z__1.i;
00353 rnrm = dznrm2_(n, &vr[ks * vr_dim1 + 1], &c__1);
00354 lnrm = dznrm2_(n, &vl[ks * vl_dim1 + 1], &c__1);
00355 s[ks] = z_abs(&prod) / (rnrm * lnrm);
00356
00357 }
00358
00359 if (wantsp) {
00360
00361
00362
00363
00364
00365
00366
00367 zlacpy_("Full", n, n, &t[t_offset], ldt, &work[work_offset],
00368 ldwork);
00369 ztrexc_("No Q", n, &work[work_offset], ldwork, dummy, &c__1, &k, &
00370 c__1, &ierr);
00371
00372
00373
00374 i__2 = *n;
00375 for (i__ = 2; i__ <= i__2; ++i__) {
00376 i__3 = i__ + i__ * work_dim1;
00377 i__4 = i__ + i__ * work_dim1;
00378 i__5 = work_dim1 + 1;
00379 z__1.r = work[i__4].r - work[i__5].r, z__1.i = work[i__4].i -
00380 work[i__5].i;
00381 work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00382
00383 }
00384
00385
00386
00387
00388 sep[ks] = 0.;
00389 est = 0.;
00390 kase = 0;
00391 *(unsigned char *)normin = 'N';
00392 L30:
00393 i__2 = *n - 1;
00394 zlacn2_(&i__2, &work[(*n + 1) * work_dim1 + 1], &work[work_offset]
00395 , &est, &kase, isave);
00396
00397 if (kase != 0) {
00398 if (kase == 1) {
00399
00400
00401
00402 i__2 = *n - 1;
00403 zlatrs_("Upper", "Conjugate transpose", "Nonunit", normin,
00404 &i__2, &work[(work_dim1 << 1) + 2], ldwork, &
00405 work[work_offset], &scale, &rwork[1], &ierr);
00406 } else {
00407
00408
00409
00410 i__2 = *n - 1;
00411 zlatrs_("Upper", "No transpose", "Nonunit", normin, &i__2,
00412 &work[(work_dim1 << 1) + 2], ldwork, &work[
00413 work_offset], &scale, &rwork[1], &ierr);
00414 }
00415 *(unsigned char *)normin = 'Y';
00416 if (scale != 1.) {
00417
00418
00419
00420
00421 i__2 = *n - 1;
00422 ix = izamax_(&i__2, &work[work_offset], &c__1);
00423 i__2 = ix + work_dim1;
00424 xnorm = (d__1 = work[i__2].r, abs(d__1)) + (d__2 = d_imag(
00425 &work[ix + work_dim1]), abs(d__2));
00426 if (scale < xnorm * smlnum || scale == 0.) {
00427 goto L40;
00428 }
00429 zdrscl_(n, &scale, &work[work_offset], &c__1);
00430 }
00431 goto L30;
00432 }
00433
00434 sep[ks] = 1. / max(est,smlnum);
00435 }
00436
00437 L40:
00438 ++ks;
00439 L50:
00440 ;
00441 }
00442 return 0;
00443
00444
00445
00446 }