Télécharger dlasy2.eso

Retour à la liste

Numérotation des lignes :

  1. C DLASY2 SOURCE BP208322 15/12/17 21:15:07 8750
  2. *> \brief \b DLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DLASY2 + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasy2.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasy2.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasy2.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
  23. * LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
  24. *
  25. * .. Scalar Arguments ..
  26. * LOGICAL LTRANL, LTRANR
  27. * INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
  28. * REAL*8 SCALE, XNORM
  29. * ..
  30. * .. Array Arguments ..
  31. * REAL*8 B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
  32. * $ X( LDX, * )
  33. * ..
  34. *
  35. *
  36. *> \par Purpose:
  37. * =============
  38. *>
  39. *> \verbatim
  40. *>
  41. *> DLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
  42. *>
  43. *> op(TL)*X + ISGN*X*op(TR) = SCALE*B,
  44. *>
  45. *> where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
  46. *> -1. op(T) = T or T**T, where T**T denotes the transpose of T.
  47. *> \endverbatim
  48. *
  49. * Arguments:
  50. * ==========
  51. *
  52. *> \param[in] LTRANL
  53. *> \verbatim
  54. *> LTRANL is LOGICAL
  55. *> On entry, LTRANL specifies the op(TL):
  56. *> = .FALSE., op(TL) = TL,
  57. *> = .TRUE., op(TL) = TL**T.
  58. *> \endverbatim
  59. *>
  60. *> \param[in] LTRANR
  61. *> \verbatim
  62. *> LTRANR is LOGICAL
  63. *> On entry, LTRANR specifies the op(TR):
  64. *> = .FALSE., op(TR) = TR,
  65. *> = .TRUE., op(TR) = TR**T.
  66. *> \endverbatim
  67. *>
  68. *> \param[in] ISGN
  69. *> \verbatim
  70. *> ISGN is INTEGER
  71. *> On entry, ISGN specifies the sign of the equation
  72. *> as described before. ISGN may only be 1 or -1.
  73. *> \endverbatim
  74. *>
  75. *> \param[in] N1
  76. *> \verbatim
  77. *> N1 is INTEGER
  78. *> On entry, N1 specifies the order of matrix TL.
  79. *> N1 may only be 0, 1 or 2.
  80. *> \endverbatim
  81. *>
  82. *> \param[in] N2
  83. *> \verbatim
  84. *> N2 is INTEGER
  85. *> On entry, N2 specifies the order of matrix TR.
  86. *> N2 may only be 0, 1 or 2.
  87. *> \endverbatim
  88. *>
  89. *> \param[in] TL
  90. *> \verbatim
  91. *> TL is DOUBLE PRECISION array, dimension (LDTL,2)
  92. *> On entry, TL contains an N1 by N1 matrix.
  93. *> \endverbatim
  94. *>
  95. *> \param[in] LDTL
  96. *> \verbatim
  97. *> LDTL is INTEGER
  98. *> The leading dimension of the matrix TL. LDTL >= max(1,N1).
  99. *> \endverbatim
  100. *>
  101. *> \param[in] TR
  102. *> \verbatim
  103. *> TR is DOUBLE PRECISION array, dimension (LDTR,2)
  104. *> On entry, TR contains an N2 by N2 matrix.
  105. *> \endverbatim
  106. *>
  107. *> \param[in] LDTR
  108. *> \verbatim
  109. *> LDTR is INTEGER
  110. *> The leading dimension of the matrix TR. LDTR >= max(1,N2).
  111. *> \endverbatim
  112. *>
  113. *> \param[in] B
  114. *> \verbatim
  115. *> B is DOUBLE PRECISION array, dimension (LDB,2)
  116. *> On entry, the N1 by N2 matrix B contains the right-hand
  117. *> side of the equation.
  118. *> \endverbatim
  119. *>
  120. *> \param[in] LDB
  121. *> \verbatim
  122. *> LDB is INTEGER
  123. *> The leading dimension of the matrix B. LDB >= max(1,N1).
  124. *> \endverbatim
  125. *>
  126. *> \param[out] SCALE
  127. *> \verbatim
  128. *> SCALE is DOUBLE PRECISION
  129. *> On exit, SCALE contains the scale factor. SCALE is chosen
  130. *> less than or equal to 1 to prevent the solution overflowing.
  131. *> \endverbatim
  132. *>
  133. *> \param[out] X
  134. *> \verbatim
  135. *> X is DOUBLE PRECISION array, dimension (LDX,2)
  136. *> On exit, X contains the N1 by N2 solution.
  137. *> \endverbatim
  138. *>
  139. *> \param[in] LDX
  140. *> \verbatim
  141. *> LDX is INTEGER
  142. *> The leading dimension of the matrix X. LDX >= max(1,N1).
  143. *> \endverbatim
  144. *>
  145. *> \param[out] XNORM
  146. *> \verbatim
  147. *> XNORM is DOUBLE PRECISION
  148. *> On exit, XNORM is the infinity-norm of the solution.
  149. *> \endverbatim
  150. *>
  151. *> \param[out] INFO
  152. *> \verbatim
  153. *> INFO is INTEGER
  154. *> On exit, INFO is set to
  155. *> 0: successful exit.
  156. *> 1: TL and TR have too close eigenvalues, so TL or
  157. *> TR is perturbed to get a nonsingular equation.
  158. *> NOTE: In the interests of speed, this routine does not
  159. *> check the inputs for errors.
  160. *> \endverbatim
  161. *
  162. * Authors:
  163. * ========
  164. *
  165. *> \author Univ. of Tennessee
  166. *> \author Univ. of California Berkeley
  167. *> \author Univ. of Colorado Denver
  168. *> \author NAG Ltd.
  169. *
  170. *> \date September 2012
  171. *
  172. *> \ingroup doubleSYauxiliary
  173. *
  174. * =====================================================================
  175. SUBROUTINE DLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
  176. $ LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
  177. *
  178. * -- LAPACK auxiliary routine (version 3.4.2) --
  179. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  180. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  181. * September 2012
  182. *
  183. * .. Scalar Arguments ..
  184. LOGICAL LTRANL, LTRANR
  185. INTEGER INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
  186. REAL*8 SCALE, XNORM
  187. * ..
  188. * .. Array Arguments ..
  189. REAL*8 B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
  190. $ X( LDX, * )
  191. * ..
  192. *
  193. * =====================================================================
  194. *
  195. * .. Parameters ..
  196. REAL*8 ZERO, ONE
  197. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  198. REAL*8 TWO, HALF, EIGHT
  199. PARAMETER ( TWO = 2.0D+0, HALF = 0.5D+0, EIGHT = 8.0D+0 )
  200. * ..
  201. * .. Local Scalars ..
  202. LOGICAL BSWAP, XSWAP
  203. INTEGER I, IP, IPIV, IPSV, J, JP, JPSV, K
  204. REAL*8 BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1,
  205. $ TEMP, U11, U12, U22, XMAX
  206. * ..
  207. * .. Local Arrays ..
  208. LOGICAL BSWPIV( 4 ), XSWPIV( 4 )
  209. INTEGER JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ),
  210. $ LOCU22( 4 )
  211. REAL*8 BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
  212. * ..
  213. * .. External Functions ..
  214. INTEGER IDAMAX
  215. REAL*8 DLAMCH
  216. EXTERNAL IDAMAX, DLAMCH
  217. * ..
  218. * .. External Subroutines ..
  219. EXTERNAL DCOPY, DSWAP
  220. * ..
  221. ** .. Intrinsic Functions ..
  222. * INTRINSIC ABS, MAX
  223. ** ..
  224. ** .. Data statements ..
  225. DATA LOCU12 / 3, 4, 1, 2 /
  226. DATA LOCL21 / 2, 1, 4, 3 /
  227. DATA LOCU22 / 4, 3, 2, 1 /
  228. DATA XSWPIV / .FALSE., .FALSE., .TRUE., .TRUE. /
  229. DATA BSWPIV / .FALSE., .TRUE., .FALSE., .TRUE. /
  230. ** ..
  231. ** .. Executable Statements ..
  232. *
  233. * Do not check the input parameters for errors
  234. *
  235. INFO = 0
  236. *
  237. * Quick return if possible
  238. *
  239. IF( N1.EQ.0 .OR. N2.EQ.0 )
  240. $ RETURN
  241. *
  242. * Set constants to control overflow
  243. *
  244. EPS = DLAMCH( 'P' )
  245. SMLNUM = DLAMCH( 'S' ) / EPS
  246. SGN = ISGN
  247. *
  248. K = N1 + N1 + N2 - 2
  249. GO TO ( 10, 20, 30, 50 )K
  250. *
  251. * 1 by 1: TL11*X + SGN*X*TR11 = B11
  252. *
  253. 10 CONTINUE
  254. TAU1 = TL( 1, 1 ) + SGN*TR( 1, 1 )
  255. BET = ABS( TAU1 )
  256. IF( BET.LE.SMLNUM ) THEN
  257. TAU1 = SMLNUM
  258. BET = SMLNUM
  259. INFO = 1
  260. END IF
  261. *
  262. SCALE = ONE
  263. GAM = ABS( B( 1, 1 ) )
  264. IF( SMLNUM*GAM.GT.BET )
  265. $ SCALE = ONE / GAM
  266. *
  267. X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / TAU1
  268. XNORM = ABS( X( 1, 1 ) )
  269. RETURN
  270. *
  271. * 1 by 2:
  272. * TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12] = [B11 B12]
  273. * [TR21 TR22]
  274. *
  275. 20 CONTINUE
  276. *
  277. SMIN = MAX( EPS*MAX( ABS( TL( 1, 1 ) ), ABS( TR( 1, 1 ) ),
  278. $ ABS( TR( 1, 2 ) ), ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) ),
  279. $ SMLNUM )
  280. TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
  281. TMP( 4 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
  282. IF( LTRANR ) THEN
  283. TMP( 2 ) = SGN*TR( 2, 1 )
  284. TMP( 3 ) = SGN*TR( 1, 2 )
  285. ELSE
  286. TMP( 2 ) = SGN*TR( 1, 2 )
  287. TMP( 3 ) = SGN*TR( 2, 1 )
  288. END IF
  289. BTMP( 1 ) = B( 1, 1 )
  290. BTMP( 2 ) = B( 1, 2 )
  291. GO TO 40
  292. *
  293. * 2 by 1:
  294. * op[TL11 TL12]*[X11] + ISGN* [X11]*TR11 = [B11]
  295. * [TL21 TL22] [X21] [X21] [B21]
  296. *
  297. 30 CONTINUE
  298. SMIN = MAX( EPS*MAX( ABS( TR( 1, 1 ) ), ABS( TL( 1, 1 ) ),
  299. $ ABS( TL( 1, 2 ) ), ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) ),
  300. $ SMLNUM )
  301. TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
  302. TMP( 4 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
  303. IF( LTRANL ) THEN
  304. TMP( 2 ) = TL( 1, 2 )
  305. TMP( 3 ) = TL( 2, 1 )
  306. ELSE
  307. TMP( 2 ) = TL( 2, 1 )
  308. TMP( 3 ) = TL( 1, 2 )
  309. END IF
  310. BTMP( 1 ) = B( 1, 1 )
  311. BTMP( 2 ) = B( 2, 1 )
  312. 40 CONTINUE
  313. *
  314. * Solve 2 by 2 system using complete pivoting.
  315. * Set pivots less than SMIN to SMIN.
  316. *
  317. IPIV = IDAMAX( 4, TMP, 1 )
  318. U11 = TMP( IPIV )
  319. IF( ABS( U11 ).LE.SMIN ) THEN
  320. INFO = 1
  321. U11 = SMIN
  322. END IF
  323. U12 = TMP( LOCU12( IPIV ) )
  324. L21 = TMP( LOCL21( IPIV ) ) / U11
  325. U22 = TMP( LOCU22( IPIV ) ) - U12*L21
  326. XSWAP = XSWPIV( IPIV )
  327. BSWAP = BSWPIV( IPIV )
  328. IF( ABS( U22 ).LE.SMIN ) THEN
  329. INFO = 1
  330. U22 = SMIN
  331. END IF
  332. IF( BSWAP ) THEN
  333. TEMP = BTMP( 2 )
  334. BTMP( 2 ) = BTMP( 1 ) - L21*TEMP
  335. BTMP( 1 ) = TEMP
  336. ELSE
  337. BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 )
  338. END IF
  339. SCALE = ONE
  340. IF( ( TWO*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( U22 ) .OR.
  341. $ ( TWO*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( U11 ) ) THEN
  342. SCALE = HALF / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ) )
  343. BTMP( 1 ) = BTMP( 1 )*SCALE
  344. BTMP( 2 ) = BTMP( 2 )*SCALE
  345. END IF
  346. X2( 2 ) = BTMP( 2 ) / U22
  347. X2( 1 ) = BTMP( 1 ) / U11 - ( U12 / U11 )*X2( 2 )
  348. IF( XSWAP ) THEN
  349. TEMP = X2( 2 )
  350. X2( 2 ) = X2( 1 )
  351. X2( 1 ) = TEMP
  352. END IF
  353. X( 1, 1 ) = X2( 1 )
  354. IF( N1.EQ.1 ) THEN
  355. X( 1, 2 ) = X2( 2 )
  356. XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
  357. ELSE
  358. X( 2, 1 ) = X2( 2 )
  359. XNORM = MAX( ABS( X( 1, 1 ) ), ABS( X( 2, 1 ) ) )
  360. END IF
  361. RETURN
  362. *
  363. * 2 by 2:
  364. * op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
  365. * [TL21 TL22] [X21 X22] [X21 X22] [TR21 TR22] [B21 B22]
  366. *
  367. * Solve equivalent 4 by 4 system using complete pivoting.
  368. * Set pivots less than SMIN to SMIN.
  369. *
  370. 50 CONTINUE
  371. SMIN = MAX( ABS( TR( 1, 1 ) ), ABS( TR( 1, 2 ) ),
  372. $ ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) )
  373. SMIN = MAX( SMIN, ABS( TL( 1, 1 ) ), ABS( TL( 1, 2 ) ),
  374. $ ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) )
  375. SMIN = MAX( EPS*SMIN, SMLNUM )
  376. BTMP( 1 ) = ZERO
  377. CALL DCOPY( 16, BTMP, 0, T16, 1 )
  378. T16( 1, 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
  379. T16( 2, 2 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
  380. T16( 3, 3 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
  381. T16( 4, 4 ) = TL( 2, 2 ) + SGN*TR( 2, 2 )
  382. IF( LTRANL ) THEN
  383. T16( 1, 2 ) = TL( 2, 1 )
  384. T16( 2, 1 ) = TL( 1, 2 )
  385. T16( 3, 4 ) = TL( 2, 1 )
  386. T16( 4, 3 ) = TL( 1, 2 )
  387. ELSE
  388. T16( 1, 2 ) = TL( 1, 2 )
  389. T16( 2, 1 ) = TL( 2, 1 )
  390. T16( 3, 4 ) = TL( 1, 2 )
  391. T16( 4, 3 ) = TL( 2, 1 )
  392. END IF
  393. IF( LTRANR ) THEN
  394. T16( 1, 3 ) = SGN*TR( 1, 2 )
  395. T16( 2, 4 ) = SGN*TR( 1, 2 )
  396. T16( 3, 1 ) = SGN*TR( 2, 1 )
  397. T16( 4, 2 ) = SGN*TR( 2, 1 )
  398. ELSE
  399. T16( 1, 3 ) = SGN*TR( 2, 1 )
  400. T16( 2, 4 ) = SGN*TR( 2, 1 )
  401. T16( 3, 1 ) = SGN*TR( 1, 2 )
  402. T16( 4, 2 ) = SGN*TR( 1, 2 )
  403. END IF
  404. BTMP( 1 ) = B( 1, 1 )
  405. BTMP( 2 ) = B( 2, 1 )
  406. BTMP( 3 ) = B( 1, 2 )
  407. BTMP( 4 ) = B( 2, 2 )
  408. *
  409. * Perform elimination
  410. *
  411. DO 100 I = 1, 3
  412. XMAX = ZERO
  413. DO 70 IP = I, 4
  414. DO 60 JP = I, 4
  415. IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN
  416. XMAX = ABS( T16( IP, JP ) )
  417. IPSV = IP
  418. JPSV = JP
  419. END IF
  420. 60 CONTINUE
  421. 70 CONTINUE
  422. IF( IPSV.NE.I ) THEN
  423. CALL DSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 )
  424. TEMP = BTMP( I )
  425. BTMP( I ) = BTMP( IPSV )
  426. BTMP( IPSV ) = TEMP
  427. END IF
  428. IF( JPSV.NE.I )
  429. $ CALL DSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 )
  430. JPIV( I ) = JPSV
  431. IF( ABS( T16( I, I ) ).LT.SMIN ) THEN
  432. INFO = 1
  433. T16( I, I ) = SMIN
  434. END IF
  435. DO 90 J = I + 1, 4
  436. T16( J, I ) = T16( J, I ) / T16( I, I )
  437. BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I )
  438. DO 80 K = I + 1, 4
  439. T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K )
  440. 80 CONTINUE
  441. 90 CONTINUE
  442. 100 CONTINUE
  443. IF( ABS( T16( 4, 4 ) ).LT.SMIN )
  444. $ T16( 4, 4 ) = SMIN
  445. SCALE = ONE
  446. IF( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
  447. $ ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
  448. $ ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
  449. $ ( EIGHT*SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN
  450. SCALE = ( ONE / EIGHT ) / MAX( ABS( BTMP( 1 ) ),
  451. $ ABS( BTMP( 2 ) ), ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) )
  452. BTMP( 1 ) = BTMP( 1 )*SCALE
  453. BTMP( 2 ) = BTMP( 2 )*SCALE
  454. BTMP( 3 ) = BTMP( 3 )*SCALE
  455. BTMP( 4 ) = BTMP( 4 )*SCALE
  456. END IF
  457. DO 120 I = 1, 4
  458. K = 5 - I
  459. TEMP = ONE / T16( K, K )
  460. TMP( K ) = BTMP( K )*TEMP
  461. DO 110 J = K + 1, 4
  462. TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J )
  463. 110 CONTINUE
  464. 120 CONTINUE
  465. DO 130 I = 1, 3
  466. IF( JPIV( 4-I ).NE.4-I ) THEN
  467. TEMP = TMP( 4-I )
  468. TMP( 4-I ) = TMP( JPIV( 4-I ) )
  469. TMP( JPIV( 4-I ) ) = TEMP
  470. END IF
  471. 130 CONTINUE
  472. X( 1, 1 ) = TMP( 1 )
  473. X( 2, 1 ) = TMP( 2 )
  474. X( 1, 2 ) = TMP( 3 )
  475. X( 2, 2 ) = TMP( 4 )
  476. XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ),
  477. $ ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) )
  478. RETURN
  479. *
  480. * End of DLASY2
  481. *
  482. END
  483.  
  484.  
  485.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales