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__4 = 4; 18 static integer c__1 = 1; 19 static integer c__16 = 16; 20 static integer c__0 = 0; 21 22 /* > \brief \b DLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2. 23 24 =========== DOCUMENTATION =========== 25 26 Online html documentation available at 27 http://www.netlib.org/lapack/explore-html/ 28 29 > \htmlonly 30 > Download DLASY2 + dependencies 31 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasy2. 32 f"> 33 > [TGZ]</a> 34 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasy2. 35 f"> 36 > [ZIP]</a> 37 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasy2. 38 f"> 39 > [TXT]</a> 40 > \endhtmlonly 41 42 Definition: 43 =========== 44 45 SUBROUTINE DLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, 46 LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO ) 47 48 LOGICAL LTRANL, LTRANR 49 INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2 50 DOUBLE PRECISION SCALE, XNORM 51 DOUBLE PRECISION B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ), 52 $ X( LDX, * ) 53 54 55 > \par Purpose: 56 ============= 57 > 58 > \verbatim 59 > 60 > DLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in 61 > 62 > op(TL)*X + ISGN*X*op(TR) = SCALE*B, 63 > 64 > where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or 65 > -1. op(T) = T or T**T, where T**T denotes the transpose of T. 66 > \endverbatim 67 68 Arguments: 69 ========== 70 71 > \param[in] LTRANL 72 > \verbatim 73 > LTRANL is LOGICAL 74 > On entry, LTRANL specifies the op(TL): 75 > = .FALSE., op(TL) = TL, 76 > = .TRUE., op(TL) = TL**T. 77 > \endverbatim 78 > 79 > \param[in] LTRANR 80 > \verbatim 81 > LTRANR is LOGICAL 82 > On entry, LTRANR specifies the op(TR): 83 > = .FALSE., op(TR) = TR, 84 > = .TRUE., op(TR) = TR**T. 85 > \endverbatim 86 > 87 > \param[in] ISGN 88 > \verbatim 89 > ISGN is INTEGER 90 > On entry, ISGN specifies the sign of the equation 91 > as described before. ISGN may only be 1 or -1. 92 > \endverbatim 93 > 94 > \param[in] N1 95 > \verbatim 96 > N1 is INTEGER 97 > On entry, N1 specifies the order of matrix TL. 98 > N1 may only be 0, 1 or 2. 99 > \endverbatim 100 > 101 > \param[in] N2 102 > \verbatim 103 > N2 is INTEGER 104 > On entry, N2 specifies the order of matrix TR. 105 > N2 may only be 0, 1 or 2. 106 > \endverbatim 107 > 108 > \param[in] TL 109 > \verbatim 110 > TL is DOUBLE PRECISION array, dimension (LDTL,2) 111 > On entry, TL contains an N1 by N1 matrix. 112 > \endverbatim 113 > 114 > \param[in] LDTL 115 > \verbatim 116 > LDTL is INTEGER 117 > The leading dimension of the matrix TL. LDTL >= max(1,N1). 118 > \endverbatim 119 > 120 > \param[in] TR 121 > \verbatim 122 > TR is DOUBLE PRECISION array, dimension (LDTR,2) 123 > On entry, TR contains an N2 by N2 matrix. 124 > \endverbatim 125 > 126 > \param[in] LDTR 127 > \verbatim 128 > LDTR is INTEGER 129 > The leading dimension of the matrix TR. LDTR >= max(1,N2). 130 > \endverbatim 131 > 132 > \param[in] B 133 > \verbatim 134 > B is DOUBLE PRECISION array, dimension (LDB,2) 135 > On entry, the N1 by N2 matrix B contains the right-hand 136 > side of the equation. 137 > \endverbatim 138 > 139 > \param[in] LDB 140 > \verbatim 141 > LDB is INTEGER 142 > The leading dimension of the matrix B. LDB >= max(1,N1). 143 > \endverbatim 144 > 145 > \param[out] SCALE 146 > \verbatim 147 > SCALE is DOUBLE PRECISION 148 > On exit, SCALE contains the scale factor. SCALE is chosen 149 > less than or equal to 1 to prevent the solution overflowing. 150 > \endverbatim 151 > 152 > \param[out] X 153 > \verbatim 154 > X is DOUBLE PRECISION array, dimension (LDX,2) 155 > On exit, X contains the N1 by N2 solution. 156 > \endverbatim 157 > 158 > \param[in] LDX 159 > \verbatim 160 > LDX is INTEGER 161 > The leading dimension of the matrix X. LDX >= max(1,N1). 162 > \endverbatim 163 > 164 > \param[out] XNORM 165 > \verbatim 166 > XNORM is DOUBLE PRECISION 167 > On exit, XNORM is the infinity-norm of the solution. 168 > \endverbatim 169 > 170 > \param[out] INFO 171 > \verbatim 172 > INFO is INTEGER 173 > On exit, INFO is set to 174 > 0: successful exit. 175 > 1: TL and TR have too close eigenvalues, so TL or 176 > TR is perturbed to get a nonsingular equation. 177 > NOTE: In the interests of speed, this routine does not 178 > check the inputs for errors. 179 > \endverbatim 180 181 Authors: 182 ======== 183 184 > \author Univ. of Tennessee 185 > \author Univ. of California Berkeley 186 > \author Univ. of Colorado Denver 187 > \author NAG Ltd. 188 189 > \date September 2012 190 191 > \ingroup doubleSYauxiliary 192 193 ===================================================================== igraphdlasy2_(logical * ltranl,logical * ltranr,integer * isgn,integer * n1,integer * n2,doublereal * tl,integer * ldtl,doublereal * tr,integer * ldtr,doublereal * b,integer * ldb,doublereal * scale,doublereal * x,integer * ldx,doublereal * xnorm,integer * info)194 Subroutine */ int igraphdlasy2_(logical *ltranl, logical *ltranr, integer *isgn, 195 integer *n1, integer *n2, doublereal *tl, integer *ldtl, doublereal * 196 tr, integer *ldtr, doublereal *b, integer *ldb, doublereal *scale, 197 doublereal *x, integer *ldx, doublereal *xnorm, integer *info) 198 { 199 /* Initialized data */ 200 201 static integer locu12[4] = { 3,4,1,2 }; 202 static integer locl21[4] = { 2,1,4,3 }; 203 static integer locu22[4] = { 4,3,2,1 }; 204 static logical xswpiv[4] = { FALSE_,FALSE_,TRUE_,TRUE_ }; 205 static logical bswpiv[4] = { FALSE_,TRUE_,FALSE_,TRUE_ }; 206 207 /* System generated locals */ 208 integer b_dim1, b_offset, tl_dim1, tl_offset, tr_dim1, tr_offset, x_dim1, 209 x_offset; 210 doublereal d__1, d__2, d__3, d__4, d__5, d__6, d__7, d__8; 211 212 /* Local variables */ 213 integer i__, j, k; 214 doublereal x2[2], l21, u11, u12; 215 integer ip, jp; 216 doublereal u22, t16[16] /* was [4][4] */, gam, bet, eps, sgn, tmp[4], 217 tau1, btmp[4], smin; 218 integer ipiv; 219 doublereal temp; 220 integer jpiv[4]; 221 doublereal xmax; 222 integer ipsv, jpsv; 223 logical bswap; 224 extern /* Subroutine */ int igraphdcopy_(integer *, doublereal *, integer *, 225 doublereal *, integer *), igraphdswap_(integer *, doublereal *, integer 226 *, doublereal *, integer *); 227 logical xswap; 228 extern doublereal igraphdlamch_(char *); 229 extern integer igraphidamax_(integer *, doublereal *, integer *); 230 doublereal smlnum; 231 232 233 /* -- LAPACK auxiliary routine (version 3.4.2) -- 234 -- LAPACK is a software package provided by Univ. of Tennessee, -- 235 -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 236 September 2012 237 238 239 ===================================================================== 240 241 Parameter adjustments */ 242 tl_dim1 = *ldtl; 243 tl_offset = 1 + tl_dim1; 244 tl -= tl_offset; 245 tr_dim1 = *ldtr; 246 tr_offset = 1 + tr_dim1; 247 tr -= tr_offset; 248 b_dim1 = *ldb; 249 b_offset = 1 + b_dim1; 250 b -= b_offset; 251 x_dim1 = *ldx; 252 x_offset = 1 + x_dim1; 253 x -= x_offset; 254 255 /* Function Body 256 257 Do not check the input parameters for errors */ 258 259 *info = 0; 260 261 /* Quick return if possible */ 262 263 if (*n1 == 0 || *n2 == 0) { 264 return 0; 265 } 266 267 /* Set constants to control overflow */ 268 269 eps = igraphdlamch_("P"); 270 smlnum = igraphdlamch_("S") / eps; 271 sgn = (doublereal) (*isgn); 272 273 k = *n1 + *n1 + *n2 - 2; 274 switch (k) { 275 case 1: goto L10; 276 case 2: goto L20; 277 case 3: goto L30; 278 case 4: goto L50; 279 } 280 281 /* 1 by 1: TL11*X + SGN*X*TR11 = B11 */ 282 283 L10: 284 tau1 = tl[tl_dim1 + 1] + sgn * tr[tr_dim1 + 1]; 285 bet = abs(tau1); 286 if (bet <= smlnum) { 287 tau1 = smlnum; 288 bet = smlnum; 289 *info = 1; 290 } 291 292 *scale = 1.; 293 gam = (d__1 = b[b_dim1 + 1], abs(d__1)); 294 if (smlnum * gam > bet) { 295 *scale = 1. / gam; 296 } 297 298 x[x_dim1 + 1] = b[b_dim1 + 1] * *scale / tau1; 299 *xnorm = (d__1 = x[x_dim1 + 1], abs(d__1)); 300 return 0; 301 302 /* 1 by 2: 303 TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12] = [B11 B12] 304 [TR21 TR22] */ 305 306 L20: 307 308 /* Computing MAX 309 Computing MAX */ 310 d__7 = (d__1 = tl[tl_dim1 + 1], abs(d__1)), d__8 = (d__2 = tr[tr_dim1 + 1] 311 , abs(d__2)), d__7 = max(d__7,d__8), d__8 = (d__3 = tr[(tr_dim1 << 312 1) + 1], abs(d__3)), d__7 = max(d__7,d__8), d__8 = (d__4 = tr[ 313 tr_dim1 + 2], abs(d__4)), d__7 = max(d__7,d__8), d__8 = (d__5 = 314 tr[(tr_dim1 << 1) + 2], abs(d__5)); 315 d__6 = eps * max(d__7,d__8); 316 smin = max(d__6,smlnum); 317 tmp[0] = tl[tl_dim1 + 1] + sgn * tr[tr_dim1 + 1]; 318 tmp[3] = tl[tl_dim1 + 1] + sgn * tr[(tr_dim1 << 1) + 2]; 319 if (*ltranr) { 320 tmp[1] = sgn * tr[tr_dim1 + 2]; 321 tmp[2] = sgn * tr[(tr_dim1 << 1) + 1]; 322 } else { 323 tmp[1] = sgn * tr[(tr_dim1 << 1) + 1]; 324 tmp[2] = sgn * tr[tr_dim1 + 2]; 325 } 326 btmp[0] = b[b_dim1 + 1]; 327 btmp[1] = b[(b_dim1 << 1) + 1]; 328 goto L40; 329 330 /* 2 by 1: 331 op[TL11 TL12]*[X11] + ISGN* [X11]*TR11 = [B11] 332 [TL21 TL22] [X21] [X21] [B21] */ 333 334 L30: 335 /* Computing MAX 336 Computing MAX */ 337 d__7 = (d__1 = tr[tr_dim1 + 1], abs(d__1)), d__8 = (d__2 = tl[tl_dim1 + 1] 338 , abs(d__2)), d__7 = max(d__7,d__8), d__8 = (d__3 = tl[(tl_dim1 << 339 1) + 1], abs(d__3)), d__7 = max(d__7,d__8), d__8 = (d__4 = tl[ 340 tl_dim1 + 2], abs(d__4)), d__7 = max(d__7,d__8), d__8 = (d__5 = 341 tl[(tl_dim1 << 1) + 2], abs(d__5)); 342 d__6 = eps * max(d__7,d__8); 343 smin = max(d__6,smlnum); 344 tmp[0] = tl[tl_dim1 + 1] + sgn * tr[tr_dim1 + 1]; 345 tmp[3] = tl[(tl_dim1 << 1) + 2] + sgn * tr[tr_dim1 + 1]; 346 if (*ltranl) { 347 tmp[1] = tl[(tl_dim1 << 1) + 1]; 348 tmp[2] = tl[tl_dim1 + 2]; 349 } else { 350 tmp[1] = tl[tl_dim1 + 2]; 351 tmp[2] = tl[(tl_dim1 << 1) + 1]; 352 } 353 btmp[0] = b[b_dim1 + 1]; 354 btmp[1] = b[b_dim1 + 2]; 355 L40: 356 357 /* Solve 2 by 2 system using complete pivoting. 358 Set pivots less than SMIN to SMIN. */ 359 360 ipiv = igraphidamax_(&c__4, tmp, &c__1); 361 u11 = tmp[ipiv - 1]; 362 if (abs(u11) <= smin) { 363 *info = 1; 364 u11 = smin; 365 } 366 u12 = tmp[locu12[ipiv - 1] - 1]; 367 l21 = tmp[locl21[ipiv - 1] - 1] / u11; 368 u22 = tmp[locu22[ipiv - 1] - 1] - u12 * l21; 369 xswap = xswpiv[ipiv - 1]; 370 bswap = bswpiv[ipiv - 1]; 371 if (abs(u22) <= smin) { 372 *info = 1; 373 u22 = smin; 374 } 375 if (bswap) { 376 temp = btmp[1]; 377 btmp[1] = btmp[0] - l21 * temp; 378 btmp[0] = temp; 379 } else { 380 btmp[1] -= l21 * btmp[0]; 381 } 382 *scale = 1.; 383 if (smlnum * 2. * abs(btmp[1]) > abs(u22) || smlnum * 2. * abs(btmp[0]) > 384 abs(u11)) { 385 /* Computing MAX */ 386 d__1 = abs(btmp[0]), d__2 = abs(btmp[1]); 387 *scale = .5 / max(d__1,d__2); 388 btmp[0] *= *scale; 389 btmp[1] *= *scale; 390 } 391 x2[1] = btmp[1] / u22; 392 x2[0] = btmp[0] / u11 - u12 / u11 * x2[1]; 393 if (xswap) { 394 temp = x2[1]; 395 x2[1] = x2[0]; 396 x2[0] = temp; 397 } 398 x[x_dim1 + 1] = x2[0]; 399 if (*n1 == 1) { 400 x[(x_dim1 << 1) + 1] = x2[1]; 401 *xnorm = (d__1 = x[x_dim1 + 1], abs(d__1)) + (d__2 = x[(x_dim1 << 1) 402 + 1], abs(d__2)); 403 } else { 404 x[x_dim1 + 2] = x2[1]; 405 /* Computing MAX */ 406 d__3 = (d__1 = x[x_dim1 + 1], abs(d__1)), d__4 = (d__2 = x[x_dim1 + 2] 407 , abs(d__2)); 408 *xnorm = max(d__3,d__4); 409 } 410 return 0; 411 412 /* 2 by 2: 413 op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12] 414 [TL21 TL22] [X21 X22] [X21 X22] [TR21 TR22] [B21 B22] 415 416 Solve equivalent 4 by 4 system using complete pivoting. 417 Set pivots less than SMIN to SMIN. */ 418 419 L50: 420 /* Computing MAX */ 421 d__5 = (d__1 = tr[tr_dim1 + 1], abs(d__1)), d__6 = (d__2 = tr[(tr_dim1 << 422 1) + 1], abs(d__2)), d__5 = max(d__5,d__6), d__6 = (d__3 = tr[ 423 tr_dim1 + 2], abs(d__3)), d__5 = max(d__5,d__6), d__6 = (d__4 = 424 tr[(tr_dim1 << 1) + 2], abs(d__4)); 425 smin = max(d__5,d__6); 426 /* Computing MAX */ 427 d__5 = smin, d__6 = (d__1 = tl[tl_dim1 + 1], abs(d__1)), d__5 = max(d__5, 428 d__6), d__6 = (d__2 = tl[(tl_dim1 << 1) + 1], abs(d__2)), d__5 = 429 max(d__5,d__6), d__6 = (d__3 = tl[tl_dim1 + 2], abs(d__3)), d__5 = 430 max(d__5,d__6), d__6 = (d__4 = tl[(tl_dim1 << 1) + 2], abs(d__4)) 431 ; 432 smin = max(d__5,d__6); 433 /* Computing MAX */ 434 d__1 = eps * smin; 435 smin = max(d__1,smlnum); 436 btmp[0] = 0.; 437 igraphdcopy_(&c__16, btmp, &c__0, t16, &c__1); 438 t16[0] = tl[tl_dim1 + 1] + sgn * tr[tr_dim1 + 1]; 439 t16[5] = tl[(tl_dim1 << 1) + 2] + sgn * tr[tr_dim1 + 1]; 440 t16[10] = tl[tl_dim1 + 1] + sgn * tr[(tr_dim1 << 1) + 2]; 441 t16[15] = tl[(tl_dim1 << 1) + 2] + sgn * tr[(tr_dim1 << 1) + 2]; 442 if (*ltranl) { 443 t16[4] = tl[tl_dim1 + 2]; 444 t16[1] = tl[(tl_dim1 << 1) + 1]; 445 t16[14] = tl[tl_dim1 + 2]; 446 t16[11] = tl[(tl_dim1 << 1) + 1]; 447 } else { 448 t16[4] = tl[(tl_dim1 << 1) + 1]; 449 t16[1] = tl[tl_dim1 + 2]; 450 t16[14] = tl[(tl_dim1 << 1) + 1]; 451 t16[11] = tl[tl_dim1 + 2]; 452 } 453 if (*ltranr) { 454 t16[8] = sgn * tr[(tr_dim1 << 1) + 1]; 455 t16[13] = sgn * tr[(tr_dim1 << 1) + 1]; 456 t16[2] = sgn * tr[tr_dim1 + 2]; 457 t16[7] = sgn * tr[tr_dim1 + 2]; 458 } else { 459 t16[8] = sgn * tr[tr_dim1 + 2]; 460 t16[13] = sgn * tr[tr_dim1 + 2]; 461 t16[2] = sgn * tr[(tr_dim1 << 1) + 1]; 462 t16[7] = sgn * tr[(tr_dim1 << 1) + 1]; 463 } 464 btmp[0] = b[b_dim1 + 1]; 465 btmp[1] = b[b_dim1 + 2]; 466 btmp[2] = b[(b_dim1 << 1) + 1]; 467 btmp[3] = b[(b_dim1 << 1) + 2]; 468 469 /* Perform elimination */ 470 471 for (i__ = 1; i__ <= 3; ++i__) { 472 xmax = 0.; 473 for (ip = i__; ip <= 4; ++ip) { 474 for (jp = i__; jp <= 4; ++jp) { 475 if ((d__1 = t16[ip + (jp << 2) - 5], abs(d__1)) >= xmax) { 476 xmax = (d__1 = t16[ip + (jp << 2) - 5], abs(d__1)); 477 ipsv = ip; 478 jpsv = jp; 479 } 480 /* L60: */ 481 } 482 /* L70: */ 483 } 484 if (ipsv != i__) { 485 igraphdswap_(&c__4, &t16[ipsv - 1], &c__4, &t16[i__ - 1], &c__4); 486 temp = btmp[i__ - 1]; 487 btmp[i__ - 1] = btmp[ipsv - 1]; 488 btmp[ipsv - 1] = temp; 489 } 490 if (jpsv != i__) { 491 igraphdswap_(&c__4, &t16[(jpsv << 2) - 4], &c__1, &t16[(i__ << 2) - 4], 492 &c__1); 493 } 494 jpiv[i__ - 1] = jpsv; 495 if ((d__1 = t16[i__ + (i__ << 2) - 5], abs(d__1)) < smin) { 496 *info = 1; 497 t16[i__ + (i__ << 2) - 5] = smin; 498 } 499 for (j = i__ + 1; j <= 4; ++j) { 500 t16[j + (i__ << 2) - 5] /= t16[i__ + (i__ << 2) - 5]; 501 btmp[j - 1] -= t16[j + (i__ << 2) - 5] * btmp[i__ - 1]; 502 for (k = i__ + 1; k <= 4; ++k) { 503 t16[j + (k << 2) - 5] -= t16[j + (i__ << 2) - 5] * t16[i__ + ( 504 k << 2) - 5]; 505 /* L80: */ 506 } 507 /* L90: */ 508 } 509 /* L100: */ 510 } 511 if (abs(t16[15]) < smin) { 512 t16[15] = smin; 513 } 514 *scale = 1.; 515 if (smlnum * 8. * abs(btmp[0]) > abs(t16[0]) || smlnum * 8. * abs(btmp[1]) 516 > abs(t16[5]) || smlnum * 8. * abs(btmp[2]) > abs(t16[10]) || 517 smlnum * 8. * abs(btmp[3]) > abs(t16[15])) { 518 /* Computing MAX */ 519 d__1 = abs(btmp[0]), d__2 = abs(btmp[1]), d__1 = max(d__1,d__2), d__2 520 = abs(btmp[2]), d__1 = max(d__1,d__2), d__2 = abs(btmp[3]); 521 *scale = .125 / max(d__1,d__2); 522 btmp[0] *= *scale; 523 btmp[1] *= *scale; 524 btmp[2] *= *scale; 525 btmp[3] *= *scale; 526 } 527 for (i__ = 1; i__ <= 4; ++i__) { 528 k = 5 - i__; 529 temp = 1. / t16[k + (k << 2) - 5]; 530 tmp[k - 1] = btmp[k - 1] * temp; 531 for (j = k + 1; j <= 4; ++j) { 532 tmp[k - 1] -= temp * t16[k + (j << 2) - 5] * tmp[j - 1]; 533 /* L110: */ 534 } 535 /* L120: */ 536 } 537 for (i__ = 1; i__ <= 3; ++i__) { 538 if (jpiv[4 - i__ - 1] != 4 - i__) { 539 temp = tmp[4 - i__ - 1]; 540 tmp[4 - i__ - 1] = tmp[jpiv[4 - i__ - 1] - 1]; 541 tmp[jpiv[4 - i__ - 1] - 1] = temp; 542 } 543 /* L130: */ 544 } 545 x[x_dim1 + 1] = tmp[0]; 546 x[x_dim1 + 2] = tmp[1]; 547 x[(x_dim1 << 1) + 1] = tmp[2]; 548 x[(x_dim1 << 1) + 2] = tmp[3]; 549 /* Computing MAX */ 550 d__1 = abs(tmp[0]) + abs(tmp[2]), d__2 = abs(tmp[1]) + abs(tmp[3]); 551 *xnorm = max(d__1,d__2); 552 return 0; 553 554 /* End of DLASY2 */ 555 556 } /* igraphdlasy2_ */ 557 558