1 /* -- translated by f2c (version 20191129). 2 You must link the resulting object file with libf2c: 3 on Microsoft Windows system, link with libf2c.lib; 4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 5 or, if you install libf2c.a in a standard place, with -lf2c -lm 6 -- in that order, at the end of the command line, as in 7 cc *.o -lf2c -lm 8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 9 10 http://www.netlib.org/f2c/libf2c.zip 11 */ 12 13 #include "f2c.h" 14 15 /* Table of constant values */ 16 17 static integer c__1 = 1; 18 static logical c_true = TRUE_; 19 static logical c_false = FALSE_; 20 21 /* > \brief \b DTRSNA 22 23 =========== DOCUMENTATION =========== 24 25 Online html documentation available at 26 http://www.netlib.org/lapack/explore-html/ 27 28 > \htmlonly 29 > Download DTRSNA + dependencies 30 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsna. 31 f"> 32 > [TGZ]</a> 33 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsna. 34 f"> 35 > [ZIP]</a> 36 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsna. 37 f"> 38 > [TXT]</a> 39 > \endhtmlonly 40 41 Definition: 42 =========== 43 44 SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 45 LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, 46 INFO ) 47 48 CHARACTER HOWMNY, JOB 49 INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N 50 LOGICAL SELECT( * ) 51 INTEGER IWORK( * ) 52 DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ), 53 $ VR( LDVR, * ), WORK( LDWORK, * ) 54 55 56 > \par Purpose: 57 ============= 58 > 59 > \verbatim 60 > 61 > DTRSNA estimates reciprocal condition numbers for specified 62 > eigenvalues and/or right eigenvectors of a real upper 63 > quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q 64 > orthogonal). 65 > 66 > T must be in Schur canonical form (as returned by DHSEQR), that is, 67 > block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 68 > 2-by-2 diagonal block has its diagonal elements equal and its 69 > off-diagonal elements of opposite sign. 70 > \endverbatim 71 72 Arguments: 73 ========== 74 75 > \param[in] JOB 76 > \verbatim 77 > JOB is CHARACTER*1 78 > Specifies whether condition numbers are required for 79 > eigenvalues (S) or eigenvectors (SEP): 80 > = 'E': for eigenvalues only (S); 81 > = 'V': for eigenvectors only (SEP); 82 > = 'B': for both eigenvalues and eigenvectors (S and SEP). 83 > \endverbatim 84 > 85 > \param[in] HOWMNY 86 > \verbatim 87 > HOWMNY is CHARACTER*1 88 > = 'A': compute condition numbers for all eigenpairs; 89 > = 'S': compute condition numbers for selected eigenpairs 90 > specified by the array SELECT. 91 > \endverbatim 92 > 93 > \param[in] SELECT 94 > \verbatim 95 > SELECT is LOGICAL array, dimension (N) 96 > If HOWMNY = 'S', SELECT specifies the eigenpairs for which 97 > condition numbers are required. To select condition numbers 98 > for the eigenpair corresponding to a real eigenvalue w(j), 99 > SELECT(j) must be set to .TRUE.. To select condition numbers 100 > corresponding to a complex conjugate pair of eigenvalues w(j) 101 > and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be 102 > set to .TRUE.. 103 > If HOWMNY = 'A', SELECT is not referenced. 104 > \endverbatim 105 > 106 > \param[in] N 107 > \verbatim 108 > N is INTEGER 109 > The order of the matrix T. N >= 0. 110 > \endverbatim 111 > 112 > \param[in] T 113 > \verbatim 114 > T is DOUBLE PRECISION array, dimension (LDT,N) 115 > The upper quasi-triangular matrix T, in Schur canonical form. 116 > \endverbatim 117 > 118 > \param[in] LDT 119 > \verbatim 120 > LDT is INTEGER 121 > The leading dimension of the array T. LDT >= max(1,N). 122 > \endverbatim 123 > 124 > \param[in] VL 125 > \verbatim 126 > VL is DOUBLE PRECISION array, dimension (LDVL,M) 127 > If JOB = 'E' or 'B', VL must contain left eigenvectors of T 128 > (or of any Q*T*Q**T with Q orthogonal), corresponding to the 129 > eigenpairs specified by HOWMNY and SELECT. The eigenvectors 130 > must be stored in consecutive columns of VL, as returned by 131 > DHSEIN or DTREVC. 132 > If JOB = 'V', VL is not referenced. 133 > \endverbatim 134 > 135 > \param[in] LDVL 136 > \verbatim 137 > LDVL is INTEGER 138 > The leading dimension of the array VL. 139 > LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N. 140 > \endverbatim 141 > 142 > \param[in] VR 143 > \verbatim 144 > VR is DOUBLE PRECISION array, dimension (LDVR,M) 145 > If JOB = 'E' or 'B', VR must contain right eigenvectors of T 146 > (or of any Q*T*Q**T with Q orthogonal), corresponding to the 147 > eigenpairs specified by HOWMNY and SELECT. The eigenvectors 148 > must be stored in consecutive columns of VR, as returned by 149 > DHSEIN or DTREVC. 150 > If JOB = 'V', VR is not referenced. 151 > \endverbatim 152 > 153 > \param[in] LDVR 154 > \verbatim 155 > LDVR is INTEGER 156 > The leading dimension of the array VR. 157 > LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N. 158 > \endverbatim 159 > 160 > \param[out] S 161 > \verbatim 162 > S is DOUBLE PRECISION array, dimension (MM) 163 > If JOB = 'E' or 'B', the reciprocal condition numbers of the 164 > selected eigenvalues, stored in consecutive elements of the 165 > array. For a complex conjugate pair of eigenvalues two 166 > consecutive elements of S are set to the same value. Thus 167 > S(j), SEP(j), and the j-th columns of VL and VR all 168 > correspond to the same eigenpair (but not in general the 169 > j-th eigenpair, unless all eigenpairs are selected). 170 > If JOB = 'V', S is not referenced. 171 > \endverbatim 172 > 173 > \param[out] SEP 174 > \verbatim 175 > SEP is DOUBLE PRECISION array, dimension (MM) 176 > If JOB = 'V' or 'B', the estimated reciprocal condition 177 > numbers of the selected eigenvectors, stored in consecutive 178 > elements of the array. For a complex eigenvector two 179 > consecutive elements of SEP are set to the same value. If 180 > the eigenvalues cannot be reordered to compute SEP(j), SEP(j) 181 > is set to 0; this can only occur when the true value would be 182 > very small anyway. 183 > If JOB = 'E', SEP is not referenced. 184 > \endverbatim 185 > 186 > \param[in] MM 187 > \verbatim 188 > MM is INTEGER 189 > The number of elements in the arrays S (if JOB = 'E' or 'B') 190 > and/or SEP (if JOB = 'V' or 'B'). MM >= M. 191 > \endverbatim 192 > 193 > \param[out] M 194 > \verbatim 195 > M is INTEGER 196 > The number of elements of the arrays S and/or SEP actually 197 > used to store the estimated condition numbers. 198 > If HOWMNY = 'A', M is set to N. 199 > \endverbatim 200 > 201 > \param[out] WORK 202 > \verbatim 203 > WORK is DOUBLE PRECISION array, dimension (LDWORK,N+6) 204 > If JOB = 'E', WORK is not referenced. 205 > \endverbatim 206 > 207 > \param[in] LDWORK 208 > \verbatim 209 > LDWORK is INTEGER 210 > The leading dimension of the array WORK. 211 > LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N. 212 > \endverbatim 213 > 214 > \param[out] IWORK 215 > \verbatim 216 > IWORK is INTEGER array, dimension (2*(N-1)) 217 > If JOB = 'E', IWORK is not referenced. 218 > \endverbatim 219 > 220 > \param[out] INFO 221 > \verbatim 222 > INFO is INTEGER 223 > = 0: successful exit 224 > < 0: if INFO = -i, the i-th argument had an illegal value 225 > \endverbatim 226 227 Authors: 228 ======== 229 230 > \author Univ. of Tennessee 231 > \author Univ. of California Berkeley 232 > \author Univ. of Colorado Denver 233 > \author NAG Ltd. 234 235 > \date November 2011 236 237 > \ingroup doubleOTHERcomputational 238 239 > \par Further Details: 240 ===================== 241 > 242 > \verbatim 243 > 244 > The reciprocal of the condition number of an eigenvalue lambda is 245 > defined as 246 > 247 > S(lambda) = |v**T*u| / (norm(u)*norm(v)) 248 > 249 > where u and v are the right and left eigenvectors of T corresponding 250 > to lambda; v**T denotes the transpose of v, and norm(u) 251 > denotes the Euclidean norm. These reciprocal condition numbers always 252 > lie between zero (very badly conditioned) and one (very well 253 > conditioned). If n = 1, S(lambda) is defined to be 1. 254 > 255 > An approximate error bound for a computed eigenvalue W(i) is given by 256 > 257 > EPS * norm(T) / S(i) 258 > 259 > where EPS is the machine precision. 260 > 261 > The reciprocal of the condition number of the right eigenvector u 262 > corresponding to lambda is defined as follows. Suppose 263 > 264 > T = ( lambda c ) 265 > ( 0 T22 ) 266 > 267 > Then the reciprocal condition number is 268 > 269 > SEP( lambda, T22 ) = sigma-min( T22 - lambda*I ) 270 > 271 > where sigma-min denotes the smallest singular value. We approximate 272 > the smallest singular value by the reciprocal of an estimate of the 273 > one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is 274 > defined to be abs(T(1,1)). 275 > 276 > An approximate error bound for a computed right eigenvector VR(i) 277 > is given by 278 > 279 > EPS * norm(T) / SEP(i) 280 > \endverbatim 281 > 282 ===================================================================== igraphdtrsna_(char * job,char * howmny,logical * select,integer * n,doublereal * t,integer * ldt,doublereal * vl,integer * ldvl,doublereal * vr,integer * ldvr,doublereal * s,doublereal * sep,integer * mm,integer * m,doublereal * work,integer * ldwork,integer * iwork,integer * info)283 Subroutine */ int igraphdtrsna_(char *job, char *howmny, logical *select, 284 integer *n, doublereal *t, integer *ldt, doublereal *vl, integer * 285 ldvl, doublereal *vr, integer *ldvr, doublereal *s, doublereal *sep, 286 integer *mm, integer *m, doublereal *work, integer *ldwork, integer * 287 iwork, integer *info) 288 { 289 /* System generated locals */ 290 integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, 291 work_dim1, work_offset, i__1, i__2; 292 doublereal d__1, d__2; 293 294 /* Builtin functions */ 295 double sqrt(doublereal); 296 297 /* Local variables */ 298 integer i__, j, k, n2; 299 doublereal cs; 300 integer nn, ks; 301 doublereal sn, mu, eps, est; 302 integer kase; 303 doublereal cond; 304 extern doublereal igraphddot_(integer *, doublereal *, integer *, doublereal *, 305 integer *); 306 logical pair; 307 integer ierr; 308 doublereal dumm, prod; 309 integer ifst; 310 doublereal lnrm; 311 integer ilst; 312 doublereal rnrm; 313 extern doublereal igraphdnrm2_(integer *, doublereal *, integer *); 314 doublereal prod1, prod2, scale, delta; 315 extern logical igraphlsame_(char *, char *); 316 integer isave[3]; 317 logical wants; 318 doublereal dummy[1]; 319 extern /* Subroutine */ int igraphdlacn2_(integer *, doublereal *, doublereal *, 320 integer *, doublereal *, integer *, integer *); 321 extern doublereal igraphdlapy2_(doublereal *, doublereal *); 322 extern /* Subroutine */ int igraphdlabad_(doublereal *, doublereal *); 323 extern doublereal igraphdlamch_(char *); 324 extern /* Subroutine */ int igraphdlacpy_(char *, integer *, integer *, 325 doublereal *, integer *, doublereal *, integer *), 326 igraphxerbla_(char *, integer *, ftnlen); 327 doublereal bignum; 328 logical wantbh; 329 extern /* Subroutine */ int igraphdlaqtr_(logical *, logical *, integer *, 330 doublereal *, integer *, doublereal *, doublereal *, doublereal *, 331 doublereal *, doublereal *, integer *), igraphdtrexc_(char *, integer * 332 , doublereal *, integer *, doublereal *, integer *, integer *, 333 integer *, doublereal *, integer *); 334 logical somcon; 335 doublereal smlnum; 336 logical wantsp; 337 338 339 /* -- LAPACK computational routine (version 3.4.0) -- 340 -- LAPACK is a software package provided by Univ. of Tennessee, -- 341 -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 342 November 2011 343 344 345 ===================================================================== 346 347 348 Decode and test the input parameters 349 350 Parameter adjustments */ 351 --select; 352 t_dim1 = *ldt; 353 t_offset = 1 + t_dim1; 354 t -= t_offset; 355 vl_dim1 = *ldvl; 356 vl_offset = 1 + vl_dim1; 357 vl -= vl_offset; 358 vr_dim1 = *ldvr; 359 vr_offset = 1 + vr_dim1; 360 vr -= vr_offset; 361 --s; 362 --sep; 363 work_dim1 = *ldwork; 364 work_offset = 1 + work_dim1; 365 work -= work_offset; 366 --iwork; 367 368 /* Function Body */ 369 wantbh = igraphlsame_(job, "B"); 370 wants = igraphlsame_(job, "E") || wantbh; 371 wantsp = igraphlsame_(job, "V") || wantbh; 372 373 somcon = igraphlsame_(howmny, "S"); 374 375 *info = 0; 376 if (! wants && ! wantsp) { 377 *info = -1; 378 } else if (! igraphlsame_(howmny, "A") && ! somcon) { 379 *info = -2; 380 } else if (*n < 0) { 381 *info = -4; 382 } else if (*ldt < max(1,*n)) { 383 *info = -6; 384 } else if (*ldvl < 1 || wants && *ldvl < *n) { 385 *info = -8; 386 } else if (*ldvr < 1 || wants && *ldvr < *n) { 387 *info = -10; 388 } else { 389 390 /* Set M to the number of eigenpairs for which condition numbers 391 are required, and test MM. */ 392 393 if (somcon) { 394 *m = 0; 395 pair = FALSE_; 396 i__1 = *n; 397 for (k = 1; k <= i__1; ++k) { 398 if (pair) { 399 pair = FALSE_; 400 } else { 401 if (k < *n) { 402 if (t[k + 1 + k * t_dim1] == 0.) { 403 if (select[k]) { 404 ++(*m); 405 } 406 } else { 407 pair = TRUE_; 408 if (select[k] || select[k + 1]) { 409 *m += 2; 410 } 411 } 412 } else { 413 if (select[*n]) { 414 ++(*m); 415 } 416 } 417 } 418 /* L10: */ 419 } 420 } else { 421 *m = *n; 422 } 423 424 if (*mm < *m) { 425 *info = -13; 426 } else if (*ldwork < 1 || wantsp && *ldwork < *n) { 427 *info = -16; 428 } 429 } 430 if (*info != 0) { 431 i__1 = -(*info); 432 igraphxerbla_("DTRSNA", &i__1, (ftnlen)6); 433 return 0; 434 } 435 436 /* Quick return if possible */ 437 438 if (*n == 0) { 439 return 0; 440 } 441 442 if (*n == 1) { 443 if (somcon) { 444 if (! select[1]) { 445 return 0; 446 } 447 } 448 if (wants) { 449 s[1] = 1.; 450 } 451 if (wantsp) { 452 sep[1] = (d__1 = t[t_dim1 + 1], abs(d__1)); 453 } 454 return 0; 455 } 456 457 /* Get machine constants */ 458 459 eps = igraphdlamch_("P"); 460 smlnum = igraphdlamch_("S") / eps; 461 bignum = 1. / smlnum; 462 igraphdlabad_(&smlnum, &bignum); 463 464 ks = 0; 465 pair = FALSE_; 466 i__1 = *n; 467 for (k = 1; k <= i__1; ++k) { 468 469 /* Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block. */ 470 471 if (pair) { 472 pair = FALSE_; 473 goto L60; 474 } else { 475 if (k < *n) { 476 pair = t[k + 1 + k * t_dim1] != 0.; 477 } 478 } 479 480 /* Determine whether condition numbers are required for the k-th 481 eigenpair. */ 482 483 if (somcon) { 484 if (pair) { 485 if (! select[k] && ! select[k + 1]) { 486 goto L60; 487 } 488 } else { 489 if (! select[k]) { 490 goto L60; 491 } 492 } 493 } 494 495 ++ks; 496 497 if (wants) { 498 499 /* Compute the reciprocal condition number of the k-th 500 eigenvalue. */ 501 502 if (! pair) { 503 504 /* Real eigenvalue. */ 505 506 prod = igraphddot_(n, &vr[ks * vr_dim1 + 1], &c__1, &vl[ks * 507 vl_dim1 + 1], &c__1); 508 rnrm = igraphdnrm2_(n, &vr[ks * vr_dim1 + 1], &c__1); 509 lnrm = igraphdnrm2_(n, &vl[ks * vl_dim1 + 1], &c__1); 510 s[ks] = abs(prod) / (rnrm * lnrm); 511 } else { 512 513 /* Complex eigenvalue. */ 514 515 prod1 = igraphddot_(n, &vr[ks * vr_dim1 + 1], &c__1, &vl[ks * 516 vl_dim1 + 1], &c__1); 517 prod1 += igraphddot_(n, &vr[(ks + 1) * vr_dim1 + 1], &c__1, &vl[(ks 518 + 1) * vl_dim1 + 1], &c__1); 519 prod2 = igraphddot_(n, &vl[ks * vl_dim1 + 1], &c__1, &vr[(ks + 1) * 520 vr_dim1 + 1], &c__1); 521 prod2 -= igraphddot_(n, &vl[(ks + 1) * vl_dim1 + 1], &c__1, &vr[ks * 522 vr_dim1 + 1], &c__1); 523 d__1 = igraphdnrm2_(n, &vr[ks * vr_dim1 + 1], &c__1); 524 d__2 = igraphdnrm2_(n, &vr[(ks + 1) * vr_dim1 + 1], &c__1); 525 rnrm = igraphdlapy2_(&d__1, &d__2); 526 d__1 = igraphdnrm2_(n, &vl[ks * vl_dim1 + 1], &c__1); 527 d__2 = igraphdnrm2_(n, &vl[(ks + 1) * vl_dim1 + 1], &c__1); 528 lnrm = igraphdlapy2_(&d__1, &d__2); 529 cond = igraphdlapy2_(&prod1, &prod2) / (rnrm * lnrm); 530 s[ks] = cond; 531 s[ks + 1] = cond; 532 } 533 } 534 535 if (wantsp) { 536 537 /* Estimate the reciprocal condition number of the k-th 538 eigenvector. 539 540 Copy the matrix T to the array WORK and swap the diagonal 541 block beginning at T(k,k) to the (1,1) position. */ 542 543 igraphdlacpy_("Full", n, n, &t[t_offset], ldt, &work[work_offset], 544 ldwork); 545 ifst = k; 546 ilst = 1; 547 igraphdtrexc_("No Q", n, &work[work_offset], ldwork, dummy, &c__1, & 548 ifst, &ilst, &work[(*n + 1) * work_dim1 + 1], &ierr); 549 550 if (ierr == 1 || ierr == 2) { 551 552 /* Could not swap because blocks not well separated */ 553 554 scale = 1.; 555 est = bignum; 556 } else { 557 558 /* Reordering successful */ 559 560 if (work[work_dim1 + 2] == 0.) { 561 562 /* Form C = T22 - lambda*I in WORK(2:N,2:N). */ 563 564 i__2 = *n; 565 for (i__ = 2; i__ <= i__2; ++i__) { 566 work[i__ + i__ * work_dim1] -= work[work_dim1 + 1]; 567 /* L20: */ 568 } 569 n2 = 1; 570 nn = *n - 1; 571 } else { 572 573 /* Triangularize the 2 by 2 block by unitary 574 transformation U = [ cs i*ss ] 575 [ i*ss cs ]. 576 such that the (1,1) position of WORK is complex 577 eigenvalue lambda with positive imaginary part. (2,2) 578 position of WORK is the complex eigenvalue lambda 579 with negative imaginary part. */ 580 581 mu = sqrt((d__1 = work[(work_dim1 << 1) + 1], abs(d__1))) 582 * sqrt((d__2 = work[work_dim1 + 2], abs(d__2))); 583 delta = igraphdlapy2_(&mu, &work[work_dim1 + 2]); 584 cs = mu / delta; 585 sn = -work[work_dim1 + 2] / delta; 586 587 /* Form 588 589 C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ] 590 [ mu ] 591 [ .. ] 592 [ .. ] 593 [ mu ] 594 where C**T is transpose of matrix C, 595 and RWORK is stored starting in the N+1-st column of 596 WORK. */ 597 598 i__2 = *n; 599 for (j = 3; j <= i__2; ++j) { 600 work[j * work_dim1 + 2] = cs * work[j * work_dim1 + 2] 601 ; 602 work[j + j * work_dim1] -= work[work_dim1 + 1]; 603 /* L30: */ 604 } 605 work[(work_dim1 << 1) + 2] = 0.; 606 607 work[(*n + 1) * work_dim1 + 1] = mu * 2.; 608 i__2 = *n - 1; 609 for (i__ = 2; i__ <= i__2; ++i__) { 610 work[i__ + (*n + 1) * work_dim1] = sn * work[(i__ + 1) 611 * work_dim1 + 1]; 612 /* L40: */ 613 } 614 n2 = 2; 615 nn = *n - 1 << 1; 616 } 617 618 /* Estimate norm(inv(C**T)) */ 619 620 est = 0.; 621 kase = 0; 622 L50: 623 igraphdlacn2_(&nn, &work[(*n + 2) * work_dim1 + 1], &work[(*n + 4) * 624 work_dim1 + 1], &iwork[1], &est, &kase, isave); 625 if (kase != 0) { 626 if (kase == 1) { 627 if (n2 == 1) { 628 629 /* Real eigenvalue: solve C**T*x = scale*c. */ 630 631 i__2 = *n - 1; 632 igraphdlaqtr_(&c_true, &c_true, &i__2, &work[(work_dim1 633 << 1) + 2], ldwork, dummy, &dumm, &scale, 634 &work[(*n + 4) * work_dim1 + 1], &work[(* 635 n + 6) * work_dim1 + 1], &ierr); 636 } else { 637 638 /* Complex eigenvalue: solve 639 C**T*(p+iq) = scale*(c+id) in real arithmetic. */ 640 641 i__2 = *n - 1; 642 igraphdlaqtr_(&c_true, &c_false, &i__2, &work[( 643 work_dim1 << 1) + 2], ldwork, &work[(*n + 644 1) * work_dim1 + 1], &mu, &scale, &work[(* 645 n + 4) * work_dim1 + 1], &work[(*n + 6) * 646 work_dim1 + 1], &ierr); 647 } 648 } else { 649 if (n2 == 1) { 650 651 /* Real eigenvalue: solve C*x = scale*c. */ 652 653 i__2 = *n - 1; 654 igraphdlaqtr_(&c_false, &c_true, &i__2, &work[( 655 work_dim1 << 1) + 2], ldwork, dummy, & 656 dumm, &scale, &work[(*n + 4) * work_dim1 657 + 1], &work[(*n + 6) * work_dim1 + 1], & 658 ierr); 659 } else { 660 661 /* Complex eigenvalue: solve 662 C*(p+iq) = scale*(c+id) in real arithmetic. */ 663 664 i__2 = *n - 1; 665 igraphdlaqtr_(&c_false, &c_false, &i__2, &work[( 666 work_dim1 << 1) + 2], ldwork, &work[(*n + 667 1) * work_dim1 + 1], &mu, &scale, &work[(* 668 n + 4) * work_dim1 + 1], &work[(*n + 6) * 669 work_dim1 + 1], &ierr); 670 671 } 672 } 673 674 goto L50; 675 } 676 } 677 678 sep[ks] = scale / max(est,smlnum); 679 if (pair) { 680 sep[ks + 1] = sep[ks]; 681 } 682 } 683 684 if (pair) { 685 ++ks; 686 } 687 688 L60: 689 ; 690 } 691 return 0; 692 693 /* End of DTRSNA */ 694 695 } /* igraphdtrsna_ */ 696 697