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_false = FALSE_; 19 static integer c__2 = 2; 20 static doublereal c_b21 = 1.; 21 static doublereal c_b25 = 0.; 22 static logical c_true = TRUE_; 23 24 /* > \brief \b DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system 25 of special form, in real arithmetic. 26 27 =========== DOCUMENTATION =========== 28 29 Online html documentation available at 30 http://www.netlib.org/lapack/explore-html/ 31 32 > \htmlonly 33 > Download DLAQTR + dependencies 34 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqtr. 35 f"> 36 > [TGZ]</a> 37 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqtr. 38 f"> 39 > [ZIP]</a> 40 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqtr. 41 f"> 42 > [TXT]</a> 43 > \endhtmlonly 44 45 Definition: 46 =========== 47 48 SUBROUTINE DLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, 49 INFO ) 50 51 LOGICAL LREAL, LTRAN 52 INTEGER INFO, LDT, N 53 DOUBLE PRECISION SCALE, W 54 DOUBLE PRECISION B( * ), T( LDT, * ), WORK( * ), X( * ) 55 56 57 > \par Purpose: 58 ============= 59 > 60 > \verbatim 61 > 62 > DLAQTR solves the real quasi-triangular system 63 > 64 > op(T)*p = scale*c, if LREAL = .TRUE. 65 > 66 > or the complex quasi-triangular systems 67 > 68 > op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE. 69 > 70 > in real arithmetic, where T is upper quasi-triangular. 71 > If LREAL = .FALSE., then the first diagonal block of T must be 72 > 1 by 1, B is the specially structured matrix 73 > 74 > B = [ b(1) b(2) ... b(n) ] 75 > [ w ] 76 > [ w ] 77 > [ . ] 78 > [ w ] 79 > 80 > op(A) = A or A**T, A**T denotes the transpose of 81 > matrix A. 82 > 83 > On input, X = [ c ]. On output, X = [ p ]. 84 > [ d ] [ q ] 85 > 86 > This subroutine is designed for the condition number estimation 87 > in routine DTRSNA. 88 > \endverbatim 89 90 Arguments: 91 ========== 92 93 > \param[in] LTRAN 94 > \verbatim 95 > LTRAN is LOGICAL 96 > On entry, LTRAN specifies the option of conjugate transpose: 97 > = .FALSE., op(T+i*B) = T+i*B, 98 > = .TRUE., op(T+i*B) = (T+i*B)**T. 99 > \endverbatim 100 > 101 > \param[in] LREAL 102 > \verbatim 103 > LREAL is LOGICAL 104 > On entry, LREAL specifies the input matrix structure: 105 > = .FALSE., the input is complex 106 > = .TRUE., the input is real 107 > \endverbatim 108 > 109 > \param[in] N 110 > \verbatim 111 > N is INTEGER 112 > On entry, N specifies the order of T+i*B. N >= 0. 113 > \endverbatim 114 > 115 > \param[in] T 116 > \verbatim 117 > T is DOUBLE PRECISION array, dimension (LDT,N) 118 > On entry, T contains a matrix in Schur canonical form. 119 > If LREAL = .FALSE., then the first diagonal block of T mu 120 > be 1 by 1. 121 > \endverbatim 122 > 123 > \param[in] LDT 124 > \verbatim 125 > LDT is INTEGER 126 > The leading dimension of the matrix T. LDT >= max(1,N). 127 > \endverbatim 128 > 129 > \param[in] B 130 > \verbatim 131 > B is DOUBLE PRECISION array, dimension (N) 132 > On entry, B contains the elements to form the matrix 133 > B as described above. 134 > If LREAL = .TRUE., B is not referenced. 135 > \endverbatim 136 > 137 > \param[in] W 138 > \verbatim 139 > W is DOUBLE PRECISION 140 > On entry, W is the diagonal element of the matrix B. 141 > If LREAL = .TRUE., W is not referenced. 142 > \endverbatim 143 > 144 > \param[out] SCALE 145 > \verbatim 146 > SCALE is DOUBLE PRECISION 147 > On exit, SCALE is the scale factor. 148 > \endverbatim 149 > 150 > \param[in,out] X 151 > \verbatim 152 > X is DOUBLE PRECISION array, dimension (2*N) 153 > On entry, X contains the right hand side of the system. 154 > On exit, X is overwritten by the solution. 155 > \endverbatim 156 > 157 > \param[out] WORK 158 > \verbatim 159 > WORK is DOUBLE PRECISION array, dimension (N) 160 > \endverbatim 161 > 162 > \param[out] INFO 163 > \verbatim 164 > INFO is INTEGER 165 > On exit, INFO is set to 166 > 0: successful exit. 167 > 1: the some diagonal 1 by 1 block has been perturbed by 168 > a small number SMIN to keep nonsingularity. 169 > 2: the some diagonal 2 by 2 block has been perturbed by 170 > a small number in DLALN2 to keep nonsingularity. 171 > NOTE: In the interests of speed, this routine does not 172 > check the inputs for errors. 173 > \endverbatim 174 175 Authors: 176 ======== 177 178 > \author Univ. of Tennessee 179 > \author Univ. of California Berkeley 180 > \author Univ. of Colorado Denver 181 > \author NAG Ltd. 182 183 > \date September 2012 184 185 > \ingroup doubleOTHERauxiliary 186 187 ===================================================================== igraphdlaqtr_(logical * ltran,logical * lreal,integer * n,doublereal * t,integer * ldt,doublereal * b,doublereal * w,doublereal * scale,doublereal * x,doublereal * work,integer * info)188 Subroutine */ int igraphdlaqtr_(logical *ltran, logical *lreal, integer *n, 189 doublereal *t, integer *ldt, doublereal *b, doublereal *w, doublereal 190 *scale, doublereal *x, doublereal *work, integer *info) 191 { 192 /* System generated locals */ 193 integer t_dim1, t_offset, i__1, i__2; 194 doublereal d__1, d__2, d__3, d__4, d__5, d__6; 195 196 /* Local variables */ 197 doublereal d__[4] /* was [2][2] */; 198 integer i__, j, k; 199 doublereal v[4] /* was [2][2] */, z__; 200 integer j1, j2, n1, n2; 201 doublereal si, xj, sr, rec, eps, tjj, tmp; 202 extern doublereal igraphddot_(integer *, doublereal *, integer *, doublereal *, 203 integer *); 204 integer ierr; 205 doublereal smin, xmax; 206 extern /* Subroutine */ int igraphdscal_(integer *, doublereal *, doublereal *, 207 integer *); 208 extern doublereal igraphdasum_(integer *, doublereal *, integer *); 209 extern /* Subroutine */ int igraphdaxpy_(integer *, doublereal *, doublereal *, 210 integer *, doublereal *, integer *); 211 integer jnext; 212 doublereal sminw, xnorm; 213 extern /* Subroutine */ int igraphdlaln2_(logical *, integer *, integer *, 214 doublereal *, doublereal *, doublereal *, integer *, doublereal *, 215 doublereal *, doublereal *, integer *, doublereal *, doublereal * 216 , doublereal *, integer *, doublereal *, doublereal *, integer *); 217 extern doublereal igraphdlamch_(char *), igraphdlange_(char *, integer *, 218 integer *, doublereal *, integer *, doublereal *); 219 extern integer igraphidamax_(integer *, doublereal *, integer *); 220 doublereal scaloc; 221 extern /* Subroutine */ int igraphdladiv_(doublereal *, doublereal *, 222 doublereal *, doublereal *, doublereal *, doublereal *); 223 doublereal bignum; 224 logical notran; 225 doublereal smlnum; 226 227 228 /* -- LAPACK auxiliary routine (version 3.4.2) -- 229 -- LAPACK is a software package provided by Univ. of Tennessee, -- 230 -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 231 September 2012 232 233 234 ===================================================================== 235 236 237 Do not test the input parameters for errors 238 239 Parameter adjustments */ 240 t_dim1 = *ldt; 241 t_offset = 1 + t_dim1; 242 t -= t_offset; 243 --b; 244 --x; 245 --work; 246 247 /* Function Body */ 248 notran = ! (*ltran); 249 *info = 0; 250 251 /* Quick return if possible */ 252 253 if (*n == 0) { 254 return 0; 255 } 256 257 /* Set constants to control overflow */ 258 259 eps = igraphdlamch_("P"); 260 smlnum = igraphdlamch_("S") / eps; 261 bignum = 1. / smlnum; 262 263 xnorm = igraphdlange_("M", n, n, &t[t_offset], ldt, d__); 264 if (! (*lreal)) { 265 /* Computing MAX */ 266 d__1 = xnorm, d__2 = abs(*w), d__1 = max(d__1,d__2), d__2 = igraphdlange_( 267 "M", n, &c__1, &b[1], n, d__); 268 xnorm = max(d__1,d__2); 269 } 270 /* Computing MAX */ 271 d__1 = smlnum, d__2 = eps * xnorm; 272 smin = max(d__1,d__2); 273 274 /* Compute 1-norm of each column of strictly upper triangular 275 part of T to control overflow in triangular solver. */ 276 277 work[1] = 0.; 278 i__1 = *n; 279 for (j = 2; j <= i__1; ++j) { 280 i__2 = j - 1; 281 work[j] = igraphdasum_(&i__2, &t[j * t_dim1 + 1], &c__1); 282 /* L10: */ 283 } 284 285 if (! (*lreal)) { 286 i__1 = *n; 287 for (i__ = 2; i__ <= i__1; ++i__) { 288 work[i__] += (d__1 = b[i__], abs(d__1)); 289 /* L20: */ 290 } 291 } 292 293 n2 = *n << 1; 294 n1 = *n; 295 if (! (*lreal)) { 296 n1 = n2; 297 } 298 k = igraphidamax_(&n1, &x[1], &c__1); 299 xmax = (d__1 = x[k], abs(d__1)); 300 *scale = 1.; 301 302 if (xmax > bignum) { 303 *scale = bignum / xmax; 304 igraphdscal_(&n1, scale, &x[1], &c__1); 305 xmax = bignum; 306 } 307 308 if (*lreal) { 309 310 if (notran) { 311 312 /* Solve T*p = scale*c */ 313 314 jnext = *n; 315 for (j = *n; j >= 1; --j) { 316 if (j > jnext) { 317 goto L30; 318 } 319 j1 = j; 320 j2 = j; 321 jnext = j - 1; 322 if (j > 1) { 323 if (t[j + (j - 1) * t_dim1] != 0.) { 324 j1 = j - 1; 325 jnext = j - 2; 326 } 327 } 328 329 if (j1 == j2) { 330 331 /* Meet 1 by 1 diagonal block 332 333 Scale to avoid overflow when computing 334 x(j) = b(j)/T(j,j) */ 335 336 xj = (d__1 = x[j1], abs(d__1)); 337 tjj = (d__1 = t[j1 + j1 * t_dim1], abs(d__1)); 338 tmp = t[j1 + j1 * t_dim1]; 339 if (tjj < smin) { 340 tmp = smin; 341 tjj = smin; 342 *info = 1; 343 } 344 345 if (xj == 0.) { 346 goto L30; 347 } 348 349 if (tjj < 1.) { 350 if (xj > bignum * tjj) { 351 rec = 1. / xj; 352 igraphdscal_(n, &rec, &x[1], &c__1); 353 *scale *= rec; 354 xmax *= rec; 355 } 356 } 357 x[j1] /= tmp; 358 xj = (d__1 = x[j1], abs(d__1)); 359 360 /* Scale x if necessary to avoid overflow when adding a 361 multiple of column j1 of T. */ 362 363 if (xj > 1.) { 364 rec = 1. / xj; 365 if (work[j1] > (bignum - xmax) * rec) { 366 igraphdscal_(n, &rec, &x[1], &c__1); 367 *scale *= rec; 368 } 369 } 370 if (j1 > 1) { 371 i__1 = j1 - 1; 372 d__1 = -x[j1]; 373 igraphdaxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[1] 374 , &c__1); 375 i__1 = j1 - 1; 376 k = igraphidamax_(&i__1, &x[1], &c__1); 377 xmax = (d__1 = x[k], abs(d__1)); 378 } 379 380 } else { 381 382 /* Meet 2 by 2 diagonal block 383 384 Call 2 by 2 linear system solve, to take 385 care of possible overflow by scaling factor. */ 386 387 d__[0] = x[j1]; 388 d__[1] = x[j2]; 389 igraphdlaln2_(&c_false, &c__2, &c__1, &smin, &c_b21, &t[j1 + j1 390 * t_dim1], ldt, &c_b21, &c_b21, d__, &c__2, & 391 c_b25, &c_b25, v, &c__2, &scaloc, &xnorm, &ierr); 392 if (ierr != 0) { 393 *info = 2; 394 } 395 396 if (scaloc != 1.) { 397 igraphdscal_(n, &scaloc, &x[1], &c__1); 398 *scale *= scaloc; 399 } 400 x[j1] = v[0]; 401 x[j2] = v[1]; 402 403 /* Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2)) 404 to avoid overflow in updating right-hand side. 405 406 Computing MAX */ 407 d__1 = abs(v[0]), d__2 = abs(v[1]); 408 xj = max(d__1,d__2); 409 if (xj > 1.) { 410 rec = 1. / xj; 411 /* Computing MAX */ 412 d__1 = work[j1], d__2 = work[j2]; 413 if (max(d__1,d__2) > (bignum - xmax) * rec) { 414 igraphdscal_(n, &rec, &x[1], &c__1); 415 *scale *= rec; 416 } 417 } 418 419 /* Update right-hand side */ 420 421 if (j1 > 1) { 422 i__1 = j1 - 1; 423 d__1 = -x[j1]; 424 igraphdaxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[1] 425 , &c__1); 426 i__1 = j1 - 1; 427 d__1 = -x[j2]; 428 igraphdaxpy_(&i__1, &d__1, &t[j2 * t_dim1 + 1], &c__1, &x[1] 429 , &c__1); 430 i__1 = j1 - 1; 431 k = igraphidamax_(&i__1, &x[1], &c__1); 432 xmax = (d__1 = x[k], abs(d__1)); 433 } 434 435 } 436 437 L30: 438 ; 439 } 440 441 } else { 442 443 /* Solve T**T*p = scale*c */ 444 445 jnext = 1; 446 i__1 = *n; 447 for (j = 1; j <= i__1; ++j) { 448 if (j < jnext) { 449 goto L40; 450 } 451 j1 = j; 452 j2 = j; 453 jnext = j + 1; 454 if (j < *n) { 455 if (t[j + 1 + j * t_dim1] != 0.) { 456 j2 = j + 1; 457 jnext = j + 2; 458 } 459 } 460 461 if (j1 == j2) { 462 463 /* 1 by 1 diagonal block 464 465 Scale if necessary to avoid overflow in forming the 466 right-hand side element by inner product. */ 467 468 xj = (d__1 = x[j1], abs(d__1)); 469 if (xmax > 1.) { 470 rec = 1. / xmax; 471 if (work[j1] > (bignum - xj) * rec) { 472 igraphdscal_(n, &rec, &x[1], &c__1); 473 *scale *= rec; 474 xmax *= rec; 475 } 476 } 477 478 i__2 = j1 - 1; 479 x[j1] -= igraphddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1, &x[1], & 480 c__1); 481 482 xj = (d__1 = x[j1], abs(d__1)); 483 tjj = (d__1 = t[j1 + j1 * t_dim1], abs(d__1)); 484 tmp = t[j1 + j1 * t_dim1]; 485 if (tjj < smin) { 486 tmp = smin; 487 tjj = smin; 488 *info = 1; 489 } 490 491 if (tjj < 1.) { 492 if (xj > bignum * tjj) { 493 rec = 1. / xj; 494 igraphdscal_(n, &rec, &x[1], &c__1); 495 *scale *= rec; 496 xmax *= rec; 497 } 498 } 499 x[j1] /= tmp; 500 /* Computing MAX */ 501 d__2 = xmax, d__3 = (d__1 = x[j1], abs(d__1)); 502 xmax = max(d__2,d__3); 503 504 } else { 505 506 /* 2 by 2 diagonal block 507 508 Scale if necessary to avoid overflow in forming the 509 right-hand side elements by inner product. 510 511 Computing MAX */ 512 d__3 = (d__1 = x[j1], abs(d__1)), d__4 = (d__2 = x[j2], 513 abs(d__2)); 514 xj = max(d__3,d__4); 515 if (xmax > 1.) { 516 rec = 1. / xmax; 517 /* Computing MAX */ 518 d__1 = work[j2], d__2 = work[j1]; 519 if (max(d__1,d__2) > (bignum - xj) * rec) { 520 igraphdscal_(n, &rec, &x[1], &c__1); 521 *scale *= rec; 522 xmax *= rec; 523 } 524 } 525 526 i__2 = j1 - 1; 527 d__[0] = x[j1] - igraphddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1, 528 &x[1], &c__1); 529 i__2 = j1 - 1; 530 d__[1] = x[j2] - igraphddot_(&i__2, &t[j2 * t_dim1 + 1], &c__1, 531 &x[1], &c__1); 532 533 igraphdlaln2_(&c_true, &c__2, &c__1, &smin, &c_b21, &t[j1 + j1 * 534 t_dim1], ldt, &c_b21, &c_b21, d__, &c__2, &c_b25, 535 &c_b25, v, &c__2, &scaloc, &xnorm, &ierr); 536 if (ierr != 0) { 537 *info = 2; 538 } 539 540 if (scaloc != 1.) { 541 igraphdscal_(n, &scaloc, &x[1], &c__1); 542 *scale *= scaloc; 543 } 544 x[j1] = v[0]; 545 x[j2] = v[1]; 546 /* Computing MAX */ 547 d__3 = (d__1 = x[j1], abs(d__1)), d__4 = (d__2 = x[j2], 548 abs(d__2)), d__3 = max(d__3,d__4); 549 xmax = max(d__3,xmax); 550 551 } 552 L40: 553 ; 554 } 555 } 556 557 } else { 558 559 /* Computing MAX */ 560 d__1 = eps * abs(*w); 561 sminw = max(d__1,smin); 562 if (notran) { 563 564 /* Solve (T + iB)*(p+iq) = c+id */ 565 566 jnext = *n; 567 for (j = *n; j >= 1; --j) { 568 if (j > jnext) { 569 goto L70; 570 } 571 j1 = j; 572 j2 = j; 573 jnext = j - 1; 574 if (j > 1) { 575 if (t[j + (j - 1) * t_dim1] != 0.) { 576 j1 = j - 1; 577 jnext = j - 2; 578 } 579 } 580 581 if (j1 == j2) { 582 583 /* 1 by 1 diagonal block 584 585 Scale if necessary to avoid overflow in division */ 586 587 z__ = *w; 588 if (j1 == 1) { 589 z__ = b[1]; 590 } 591 xj = (d__1 = x[j1], abs(d__1)) + (d__2 = x[*n + j1], abs( 592 d__2)); 593 tjj = (d__1 = t[j1 + j1 * t_dim1], abs(d__1)) + abs(z__); 594 tmp = t[j1 + j1 * t_dim1]; 595 if (tjj < sminw) { 596 tmp = sminw; 597 tjj = sminw; 598 *info = 1; 599 } 600 601 if (xj == 0.) { 602 goto L70; 603 } 604 605 if (tjj < 1.) { 606 if (xj > bignum * tjj) { 607 rec = 1. / xj; 608 igraphdscal_(&n2, &rec, &x[1], &c__1); 609 *scale *= rec; 610 xmax *= rec; 611 } 612 } 613 igraphdladiv_(&x[j1], &x[*n + j1], &tmp, &z__, &sr, &si); 614 x[j1] = sr; 615 x[*n + j1] = si; 616 xj = (d__1 = x[j1], abs(d__1)) + (d__2 = x[*n + j1], abs( 617 d__2)); 618 619 /* Scale x if necessary to avoid overflow when adding a 620 multiple of column j1 of T. */ 621 622 if (xj > 1.) { 623 rec = 1. / xj; 624 if (work[j1] > (bignum - xmax) * rec) { 625 igraphdscal_(&n2, &rec, &x[1], &c__1); 626 *scale *= rec; 627 } 628 } 629 630 if (j1 > 1) { 631 i__1 = j1 - 1; 632 d__1 = -x[j1]; 633 igraphdaxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[1] 634 , &c__1); 635 i__1 = j1 - 1; 636 d__1 = -x[*n + j1]; 637 igraphdaxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[* 638 n + 1], &c__1); 639 640 x[1] += b[j1] * x[*n + j1]; 641 x[*n + 1] -= b[j1] * x[j1]; 642 643 xmax = 0.; 644 i__1 = j1 - 1; 645 for (k = 1; k <= i__1; ++k) { 646 /* Computing MAX */ 647 d__3 = xmax, d__4 = (d__1 = x[k], abs(d__1)) + ( 648 d__2 = x[k + *n], abs(d__2)); 649 xmax = max(d__3,d__4); 650 /* L50: */ 651 } 652 } 653 654 } else { 655 656 /* Meet 2 by 2 diagonal block */ 657 658 d__[0] = x[j1]; 659 d__[1] = x[j2]; 660 d__[2] = x[*n + j1]; 661 d__[3] = x[*n + j2]; 662 d__1 = -(*w); 663 igraphdlaln2_(&c_false, &c__2, &c__2, &sminw, &c_b21, &t[j1 + 664 j1 * t_dim1], ldt, &c_b21, &c_b21, d__, &c__2, & 665 c_b25, &d__1, v, &c__2, &scaloc, &xnorm, &ierr); 666 if (ierr != 0) { 667 *info = 2; 668 } 669 670 if (scaloc != 1.) { 671 i__1 = *n << 1; 672 igraphdscal_(&i__1, &scaloc, &x[1], &c__1); 673 *scale = scaloc * *scale; 674 } 675 x[j1] = v[0]; 676 x[j2] = v[1]; 677 x[*n + j1] = v[2]; 678 x[*n + j2] = v[3]; 679 680 /* Scale X(J1), .... to avoid overflow in 681 updating right hand side. 682 683 Computing MAX */ 684 d__1 = abs(v[0]) + abs(v[2]), d__2 = abs(v[1]) + abs(v[3]) 685 ; 686 xj = max(d__1,d__2); 687 if (xj > 1.) { 688 rec = 1. / xj; 689 /* Computing MAX */ 690 d__1 = work[j1], d__2 = work[j2]; 691 if (max(d__1,d__2) > (bignum - xmax) * rec) { 692 igraphdscal_(&n2, &rec, &x[1], &c__1); 693 *scale *= rec; 694 } 695 } 696 697 /* Update the right-hand side. */ 698 699 if (j1 > 1) { 700 i__1 = j1 - 1; 701 d__1 = -x[j1]; 702 igraphdaxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[1] 703 , &c__1); 704 i__1 = j1 - 1; 705 d__1 = -x[j2]; 706 igraphdaxpy_(&i__1, &d__1, &t[j2 * t_dim1 + 1], &c__1, &x[1] 707 , &c__1); 708 709 i__1 = j1 - 1; 710 d__1 = -x[*n + j1]; 711 igraphdaxpy_(&i__1, &d__1, &t[j1 * t_dim1 + 1], &c__1, &x[* 712 n + 1], &c__1); 713 i__1 = j1 - 1; 714 d__1 = -x[*n + j2]; 715 igraphdaxpy_(&i__1, &d__1, &t[j2 * t_dim1 + 1], &c__1, &x[* 716 n + 1], &c__1); 717 718 x[1] = x[1] + b[j1] * x[*n + j1] + b[j2] * x[*n + j2]; 719 x[*n + 1] = x[*n + 1] - b[j1] * x[j1] - b[j2] * x[j2]; 720 721 xmax = 0.; 722 i__1 = j1 - 1; 723 for (k = 1; k <= i__1; ++k) { 724 /* Computing MAX */ 725 d__3 = (d__1 = x[k], abs(d__1)) + (d__2 = x[k + * 726 n], abs(d__2)); 727 xmax = max(d__3,xmax); 728 /* L60: */ 729 } 730 } 731 732 } 733 L70: 734 ; 735 } 736 737 } else { 738 739 /* Solve (T + iB)**T*(p+iq) = c+id */ 740 741 jnext = 1; 742 i__1 = *n; 743 for (j = 1; j <= i__1; ++j) { 744 if (j < jnext) { 745 goto L80; 746 } 747 j1 = j; 748 j2 = j; 749 jnext = j + 1; 750 if (j < *n) { 751 if (t[j + 1 + j * t_dim1] != 0.) { 752 j2 = j + 1; 753 jnext = j + 2; 754 } 755 } 756 757 if (j1 == j2) { 758 759 /* 1 by 1 diagonal block 760 761 Scale if necessary to avoid overflow in forming the 762 right-hand side element by inner product. */ 763 764 xj = (d__1 = x[j1], abs(d__1)) + (d__2 = x[j1 + *n], abs( 765 d__2)); 766 if (xmax > 1.) { 767 rec = 1. / xmax; 768 if (work[j1] > (bignum - xj) * rec) { 769 igraphdscal_(&n2, &rec, &x[1], &c__1); 770 *scale *= rec; 771 xmax *= rec; 772 } 773 } 774 775 i__2 = j1 - 1; 776 x[j1] -= igraphddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1, &x[1], & 777 c__1); 778 i__2 = j1 - 1; 779 x[*n + j1] -= igraphddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1, &x[ 780 *n + 1], &c__1); 781 if (j1 > 1) { 782 x[j1] -= b[j1] * x[*n + 1]; 783 x[*n + j1] += b[j1] * x[1]; 784 } 785 xj = (d__1 = x[j1], abs(d__1)) + (d__2 = x[j1 + *n], abs( 786 d__2)); 787 788 z__ = *w; 789 if (j1 == 1) { 790 z__ = b[1]; 791 } 792 793 /* Scale if necessary to avoid overflow in 794 complex division */ 795 796 tjj = (d__1 = t[j1 + j1 * t_dim1], abs(d__1)) + abs(z__); 797 tmp = t[j1 + j1 * t_dim1]; 798 if (tjj < sminw) { 799 tmp = sminw; 800 tjj = sminw; 801 *info = 1; 802 } 803 804 if (tjj < 1.) { 805 if (xj > bignum * tjj) { 806 rec = 1. / xj; 807 igraphdscal_(&n2, &rec, &x[1], &c__1); 808 *scale *= rec; 809 xmax *= rec; 810 } 811 } 812 d__1 = -z__; 813 igraphdladiv_(&x[j1], &x[*n + j1], &tmp, &d__1, &sr, &si); 814 x[j1] = sr; 815 x[j1 + *n] = si; 816 /* Computing MAX */ 817 d__3 = (d__1 = x[j1], abs(d__1)) + (d__2 = x[j1 + *n], 818 abs(d__2)); 819 xmax = max(d__3,xmax); 820 821 } else { 822 823 /* 2 by 2 diagonal block 824 825 Scale if necessary to avoid overflow in forming the 826 right-hand side element by inner product. 827 828 Computing MAX */ 829 d__5 = (d__1 = x[j1], abs(d__1)) + (d__2 = x[*n + j1], 830 abs(d__2)), d__6 = (d__3 = x[j2], abs(d__3)) + ( 831 d__4 = x[*n + j2], abs(d__4)); 832 xj = max(d__5,d__6); 833 if (xmax > 1.) { 834 rec = 1. / xmax; 835 /* Computing MAX */ 836 d__1 = work[j1], d__2 = work[j2]; 837 if (max(d__1,d__2) > (bignum - xj) / xmax) { 838 igraphdscal_(&n2, &rec, &x[1], &c__1); 839 *scale *= rec; 840 xmax *= rec; 841 } 842 } 843 844 i__2 = j1 - 1; 845 d__[0] = x[j1] - igraphddot_(&i__2, &t[j1 * t_dim1 + 1], &c__1, 846 &x[1], &c__1); 847 i__2 = j1 - 1; 848 d__[1] = x[j2] - igraphddot_(&i__2, &t[j2 * t_dim1 + 1], &c__1, 849 &x[1], &c__1); 850 i__2 = j1 - 1; 851 d__[2] = x[*n + j1] - igraphddot_(&i__2, &t[j1 * t_dim1 + 1], & 852 c__1, &x[*n + 1], &c__1); 853 i__2 = j1 - 1; 854 d__[3] = x[*n + j2] - igraphddot_(&i__2, &t[j2 * t_dim1 + 1], & 855 c__1, &x[*n + 1], &c__1); 856 d__[0] -= b[j1] * x[*n + 1]; 857 d__[1] -= b[j2] * x[*n + 1]; 858 d__[2] += b[j1] * x[1]; 859 d__[3] += b[j2] * x[1]; 860 861 igraphdlaln2_(&c_true, &c__2, &c__2, &sminw, &c_b21, &t[j1 + j1 862 * t_dim1], ldt, &c_b21, &c_b21, d__, &c__2, & 863 c_b25, w, v, &c__2, &scaloc, &xnorm, &ierr); 864 if (ierr != 0) { 865 *info = 2; 866 } 867 868 if (scaloc != 1.) { 869 igraphdscal_(&n2, &scaloc, &x[1], &c__1); 870 *scale = scaloc * *scale; 871 } 872 x[j1] = v[0]; 873 x[j2] = v[1]; 874 x[*n + j1] = v[2]; 875 x[*n + j2] = v[3]; 876 /* Computing MAX */ 877 d__5 = (d__1 = x[j1], abs(d__1)) + (d__2 = x[*n + j1], 878 abs(d__2)), d__6 = (d__3 = x[j2], abs(d__3)) + ( 879 d__4 = x[*n + j2], abs(d__4)), d__5 = max(d__5, 880 d__6); 881 xmax = max(d__5,xmax); 882 883 } 884 885 L80: 886 ; 887 } 888 889 } 890 891 } 892 893 return 0; 894 895 /* End of DLAQTR */ 896 897 } /* igraphdlaqtr_ */ 898 899