Télécharger dtrsyl.eso

Retour à la liste

Numérotation des lignes :

  1. C DTRSYL SOURCE BP208322 15/10/13 21:16:03 8670
  2. *> \brief \b DTRSYL
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. *> \htmlonly
  10. *> Download DTRSYL + dependencies
  11. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsyl.f">
  12. *> [TGZ]</a>
  13. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsyl.f">
  14. *> [ZIP]</a>
  15. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsyl.f">
  16. *> [TXT]</a>
  17. *> \endhtmlonly
  18. *
  19. * Definition:
  20. * ===========
  21. *
  22. * SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
  23. * LDC, SCALE, INFO )
  24. *
  25. * .. Scalar Arguments ..
  26. * CHARACTER TRANA, TRANB
  27. * INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
  28. * REAL*8 SCALE
  29. * ..
  30. * .. Array Arguments ..
  31. * REAL*8 A( LDA, * ), B( LDB, * ), C( LDC, * )
  32. * ..
  33. *
  34. *
  35. *> \par Purpose:
  36. * =============
  37. *>
  38. *> \verbatim
  39. *>
  40. *> DTRSYL solves the real Sylvester matrix equation:
  41. *>
  42. *> op(A)*X + X*op(B) = scale*C or
  43. *> op(A)*X - X*op(B) = scale*C,
  44. *>
  45. *> where op(A) = A or A**T, and A and B are both upper quasi-
  46. *> triangular. A is M-by-M and B is N-by-N; the right hand side C and
  47. *> the solution X are M-by-N; and scale is an output scale factor, set
  48. *> <= 1 to avoid overflow in X.
  49. *>
  50. *> A and B must be in Schur canonical form (as returned by DHSEQR), that
  51. *> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
  52. *> each 2-by-2 diagonal block has its diagonal elements equal and its
  53. *> off-diagonal elements of opposite sign.
  54. *> \endverbatim
  55. *
  56. * Arguments:
  57. * ==========
  58. *
  59. *> \param[in] TRANA
  60. *> \verbatim
  61. *> TRANA is CHARACTER*1
  62. *> Specifies the option op(A):
  63. *> = 'N': op(A) = A (No transpose)
  64. *> = 'T': op(A) = A**T (Transpose)
  65. *> = 'C': op(A) = A**H (Conjugate transpose = Transpose)
  66. *> \endverbatim
  67. *>
  68. *> \param[in] TRANB
  69. *> \verbatim
  70. *> TRANB is CHARACTER*1
  71. *> Specifies the option op(B):
  72. *> = 'N': op(B) = B (No transpose)
  73. *> = 'T': op(B) = B**T (Transpose)
  74. *> = 'C': op(B) = B**H (Conjugate transpose = Transpose)
  75. *> \endverbatim
  76. *>
  77. *> \param[in] ISGN
  78. *> \verbatim
  79. *> ISGN is INTEGER
  80. *> Specifies the sign in the equation:
  81. *> = +1: solve op(A)*X + X*op(B) = scale*C
  82. *> = -1: solve op(A)*X - X*op(B) = scale*C
  83. *> \endverbatim
  84. *>
  85. *> \param[in] M
  86. *> \verbatim
  87. *> M is INTEGER
  88. *> The order of the matrix A, and the number of rows in the
  89. *> matrices X and C. M >= 0.
  90. *> \endverbatim
  91. *>
  92. *> \param[in] N
  93. *> \verbatim
  94. *> N is INTEGER
  95. *> The order of the matrix B, and the number of columns in the
  96. *> matrices X and C. N >= 0.
  97. *> \endverbatim
  98. *>
  99. *> \param[in] A
  100. *> \verbatim
  101. *> A is DOUBLE PRECISION array, dimension (LDA,M)
  102. *> The upper quasi-triangular matrix A, in Schur canonical form.
  103. *> \endverbatim
  104. *>
  105. *> \param[in] LDA
  106. *> \verbatim
  107. *> LDA is INTEGER
  108. *> The leading dimension of the array A. LDA >= max(1,M).
  109. *> \endverbatim
  110. *>
  111. *> \param[in] B
  112. *> \verbatim
  113. *> B is DOUBLE PRECISION array, dimension (LDB,N)
  114. *> The upper quasi-triangular matrix B, in Schur canonical form.
  115. *> \endverbatim
  116. *>
  117. *> \param[in] LDB
  118. *> \verbatim
  119. *> LDB is INTEGER
  120. *> The leading dimension of the array B. LDB >= max(1,N).
  121. *> \endverbatim
  122. *>
  123. *> \param[in,out] C
  124. *> \verbatim
  125. *> C is DOUBLE PRECISION array, dimension (LDC,N)
  126. *> On entry, the M-by-N right hand side matrix C.
  127. *> On exit, C is overwritten by the solution matrix X.
  128. *> \endverbatim
  129. *>
  130. *> \param[in] LDC
  131. *> \verbatim
  132. *> LDC is INTEGER
  133. *> The leading dimension of the array C. LDC >= max(1,M)
  134. *> \endverbatim
  135. *>
  136. *> \param[out] SCALE
  137. *> \verbatim
  138. *> SCALE is DOUBLE PRECISION
  139. *> The scale factor, scale, set <= 1 to avoid overflow in X.
  140. *> \endverbatim
  141. *>
  142. *> \param[out] INFO
  143. *> \verbatim
  144. *> INFO is INTEGER
  145. *> = 0: successful exit
  146. *> < 0: if INFO = -i, the i-th argument had an illegal value
  147. *> = 1: A and B have common or very close eigenvalues; perturbed
  148. *> values were used to solve the equation (but the matrices
  149. *> A and B are unchanged).
  150. *> \endverbatim
  151. *
  152. * Authors:
  153. * ========
  154. *
  155. *> \author Univ. of Tennessee
  156. *> \author Univ. of California Berkeley
  157. *> \author Univ. of Colorado Denver
  158. *> \author NAG Ltd.
  159. *
  160. *> \date November 2011
  161. *
  162. *> \ingroup doubleSYcomputational
  163. *
  164. * =====================================================================
  165. SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
  166. $ LDC, SCALE, INFO )
  167. *
  168. * -- LAPACK computational routine (version 3.4.0) --
  169. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  170. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  171. * November 2011
  172. *
  173. * .. Scalar Arguments ..
  174. CHARACTER TRANA, TRANB
  175. INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
  176. REAL*8 SCALE
  177. * ..
  178. * .. Array Arguments ..
  179. REAL*8 A( LDA, * ), B( LDB, * ), C( LDC, * )
  180. * ..
  181. *
  182. * =====================================================================
  183. *
  184. * .. Parameters ..
  185. REAL*8 ZERO, ONE
  186. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  187. * ..
  188. * .. Local Scalars ..
  189. LOGICAL NOTRNA, NOTRNB
  190. INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
  191. REAL*8 A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
  192. $ SMLNUM, SUML, SUMR, XNORM
  193. * ..
  194. * .. Local Arrays ..
  195. REAL*8 DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
  196. * ..
  197. * .. External Functions ..
  198. LOGICAL LSAME
  199. REAL*8 DDOT, DLAMCH, DLANGE
  200. EXTERNAL LSAME, DDOT, DLAMCH, DLANGE
  201. * ..
  202. * .. External Subroutines ..
  203. * ..
  204. ** .. Intrinsic Functions ..
  205. * INTRINSIC ABS, DBLE, MAX, MIN
  206. ** ..
  207. ** .. Executable Statements ..
  208. *
  209. * Decode and Test input parameters
  210. *
  211. NOTRNA = LSAME( TRANA, 'N' )
  212. NOTRNB = LSAME( TRANB, 'N' )
  213. *
  214. INFO = 0
  215. IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT.
  216. $ LSAME( TRANA, 'C' ) ) THEN
  217. INFO = -1
  218. ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT.
  219. $ LSAME( TRANB, 'C' ) ) THEN
  220. INFO = -2
  221. ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
  222. INFO = -3
  223. ELSE IF( M.LT.0 ) THEN
  224. INFO = -4
  225. ELSE IF( N.LT.0 ) THEN
  226. INFO = -5
  227. ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  228. INFO = -7
  229. ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  230. INFO = -9
  231. ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
  232. INFO = -11
  233. END IF
  234. IF( INFO.NE.0 ) THEN
  235. CALL XERBLA( 'DTRSYL', -INFO )
  236. RETURN
  237. END IF
  238. *
  239. * Quick return if possible
  240. *
  241. SCALE = ONE
  242. IF( M.EQ.0 .OR. N.EQ.0 )
  243. $ RETURN
  244. *
  245. * Set constants to control overflow
  246. *
  247. EPS = DLAMCH( 'P' )
  248. SMLNUM = DLAMCH( 'S' )
  249. BIGNUM = ONE / SMLNUM
  250. CALL DLABAD( SMLNUM, BIGNUM )
  251. SMLNUM = SMLNUM*DBLE( M*N ) / EPS
  252. BIGNUM = ONE / SMLNUM
  253. *
  254. SMIN = MAX( SMLNUM, EPS*DLANGE( 'M', M, M, A, LDA, DUM ),
  255. $ EPS*DLANGE( 'M', N, N, B, LDB, DUM ) )
  256. *
  257. SGN = ISGN
  258. *
  259. IF( NOTRNA .AND. NOTRNB ) THEN
  260. *
  261. * Solve A*X + ISGN*X*B = scale*C.
  262. *
  263. * The (K,L)th block of X is determined starting from
  264. * bottom-left corner column by column by
  265. *
  266. * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
  267. *
  268. * Where
  269. * M L-1
  270. * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
  271. * I=K+1 J=1
  272. *
  273. * Start column loop (index = L)
  274. * L1 (L2) : column index of the first (first) row of X(K,L).
  275. *
  276. LNEXT = 1
  277. DO 60 L = 1, N
  278. IF( L.LT.LNEXT )
  279. $ GO TO 60
  280. IF( L.EQ.N ) THEN
  281. L1 = L
  282. L2 = L
  283. ELSE
  284. IF( B( L+1, L ).NE.ZERO ) THEN
  285. L1 = L
  286. L2 = L + 1
  287. LNEXT = L + 2
  288. ELSE
  289. L1 = L
  290. L2 = L
  291. LNEXT = L + 1
  292. END IF
  293. END IF
  294. *
  295. * Start row loop (index = K)
  296. * K1 (K2): row index of the first (last) row of X(K,L).
  297. *
  298. KNEXT = M
  299. DO 50 K = M, 1, -1
  300. IF( K.GT.KNEXT )
  301. $ GO TO 50
  302. IF( K.EQ.1 ) THEN
  303. K1 = K
  304. K2 = K
  305. ELSE
  306. IF( A( K, K-1 ).NE.ZERO ) THEN
  307. K1 = K - 1
  308. K2 = K
  309. KNEXT = K - 2
  310. ELSE
  311. K1 = K
  312. K2 = K
  313. KNEXT = K - 1
  314. END IF
  315. END IF
  316. *
  317. IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
  318. SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  319. $ C( MIN( K1+1, M ), L1 ), 1 )
  320. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  321. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  322. SCALOC = ONE
  323. *
  324. A11 = A( K1, K1 ) + SGN*B( L1, L1 )
  325. DA11 = ABS( A11 )
  326. IF( DA11.LE.SMIN ) THEN
  327. A11 = SMIN
  328. DA11 = SMIN
  329. INFO = 1
  330. END IF
  331. DB = ABS( VEC( 1, 1 ) )
  332. IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  333. IF( DB.GT.BIGNUM*DA11 )
  334. $ SCALOC = ONE / DB
  335. END IF
  336. X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
  337. *
  338. IF( SCALOC.NE.ONE ) THEN
  339. DO 10 J = 1, N
  340. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  341. 10 CONTINUE
  342. SCALE = SCALE*SCALOC
  343. END IF
  344. C( K1, L1 ) = X( 1, 1 )
  345. *
  346. ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
  347. *
  348. SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  349. $ C( MIN( K2+1, M ), L1 ), 1 )
  350. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  351. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  352. *
  353. SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  354. $ C( MIN( K2+1, M ), L1 ), 1 )
  355. SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
  356. VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  357. *
  358. CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
  359. $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
  360. $ ZERO, X, 2, SCALOC, XNORM, IERR )
  361. IF( IERR.NE.0 )
  362. $ INFO = 1
  363. *
  364. IF( SCALOC.NE.ONE ) THEN
  365. DO 20 J = 1, N
  366. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  367. 20 CONTINUE
  368. SCALE = SCALE*SCALOC
  369. END IF
  370. C( K1, L1 ) = X( 1, 1 )
  371. C( K2, L1 ) = X( 2, 1 )
  372. *
  373. ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
  374. *
  375. SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  376. $ C( MIN( K1+1, M ), L1 ), 1 )
  377. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  378. VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
  379. *
  380. SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  381. $ C( MIN( K1+1, M ), L2 ), 1 )
  382. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
  383. VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
  384. *
  385. CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
  386. $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
  387. $ ZERO, X, 2, SCALOC, XNORM, IERR )
  388. IF( IERR.NE.0 )
  389. $ INFO = 1
  390. *
  391. IF( SCALOC.NE.ONE ) THEN
  392. DO 30 J = 1, N
  393. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  394. 30 CONTINUE
  395. SCALE = SCALE*SCALOC
  396. END IF
  397. C( K1, L1 ) = X( 1, 1 )
  398. C( K1, L2 ) = X( 2, 1 )
  399. *
  400. ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
  401. *
  402. SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  403. $ C( MIN( K2+1, M ), L1 ), 1 )
  404. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  405. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  406. *
  407. SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  408. $ C( MIN( K2+1, M ), L2 ), 1 )
  409. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
  410. VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
  411. *
  412. SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  413. $ C( MIN( K2+1, M ), L1 ), 1 )
  414. SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
  415. VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  416. *
  417. SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  418. $ C( MIN( K2+1, M ), L2 ), 1 )
  419. SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
  420. VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
  421. *
  422. CALL DLASY2( .FALSE., .FALSE., ISGN, 2, 2,
  423. $ A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC,
  424. $ 2, SCALOC, X, 2, XNORM, IERR )
  425. IF( IERR.NE.0 )
  426. $ INFO = 1
  427. *
  428. IF( SCALOC.NE.ONE ) THEN
  429. DO 40 J = 1, N
  430. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  431. 40 CONTINUE
  432. SCALE = SCALE*SCALOC
  433. END IF
  434. C( K1, L1 ) = X( 1, 1 )
  435. C( K1, L2 ) = X( 1, 2 )
  436. C( K2, L1 ) = X( 2, 1 )
  437. C( K2, L2 ) = X( 2, 2 )
  438. END IF
  439. *
  440. 50 CONTINUE
  441. *
  442. 60 CONTINUE
  443. *
  444. ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
  445. *
  446. * Solve A**T *X + ISGN*X*B = scale*C.
  447. *
  448. * The (K,L)th block of X is determined starting from
  449. * upper-left corner column by column by
  450. *
  451. * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
  452. *
  453. * Where
  454. * K-1 T L-1
  455. * R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
  456. * I=1 J=1
  457. *
  458. * Start column loop (index = L)
  459. * L1 (L2): column index of the first (last) row of X(K,L)
  460. *
  461. LNEXT = 1
  462. DO 120 L = 1, N
  463. IF( L.LT.LNEXT )
  464. $ GO TO 120
  465. IF( L.EQ.N ) THEN
  466. L1 = L
  467. L2 = L
  468. ELSE
  469. IF( B( L+1, L ).NE.ZERO ) THEN
  470. L1 = L
  471. L2 = L + 1
  472. LNEXT = L + 2
  473. ELSE
  474. L1 = L
  475. L2 = L
  476. LNEXT = L + 1
  477. END IF
  478. END IF
  479. *
  480. * Start row loop (index = K)
  481. * K1 (K2): row index of the first (last) row of X(K,L)
  482. *
  483. KNEXT = 1
  484. DO 110 K = 1, M
  485. IF( K.LT.KNEXT )
  486. $ GO TO 110
  487. IF( K.EQ.M ) THEN
  488. K1 = K
  489. K2 = K
  490. ELSE
  491. IF( A( K+1, K ).NE.ZERO ) THEN
  492. K1 = K
  493. K2 = K + 1
  494. KNEXT = K + 2
  495. ELSE
  496. K1 = K
  497. K2 = K
  498. KNEXT = K + 1
  499. END IF
  500. END IF
  501. *
  502. IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
  503. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  504. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  505. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  506. SCALOC = ONE
  507. *
  508. A11 = A( K1, K1 ) + SGN*B( L1, L1 )
  509. DA11 = ABS( A11 )
  510. IF( DA11.LE.SMIN ) THEN
  511. A11 = SMIN
  512. DA11 = SMIN
  513. INFO = 1
  514. END IF
  515. DB = ABS( VEC( 1, 1 ) )
  516. IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  517. IF( DB.GT.BIGNUM*DA11 )
  518. $ SCALOC = ONE / DB
  519. END IF
  520. X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
  521. *
  522. IF( SCALOC.NE.ONE ) THEN
  523. DO 70 J = 1, N
  524. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  525. 70 CONTINUE
  526. SCALE = SCALE*SCALOC
  527. END IF
  528. C( K1, L1 ) = X( 1, 1 )
  529. *
  530. ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
  531. *
  532. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  533. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  534. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  535. *
  536. SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
  537. SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
  538. VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  539. *
  540. CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
  541. $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
  542. $ ZERO, X, 2, SCALOC, XNORM, IERR )
  543. IF( IERR.NE.0 )
  544. $ INFO = 1
  545. *
  546. IF( SCALOC.NE.ONE ) THEN
  547. DO 80 J = 1, N
  548. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  549. 80 CONTINUE
  550. SCALE = SCALE*SCALOC
  551. END IF
  552. C( K1, L1 ) = X( 1, 1 )
  553. C( K2, L1 ) = X( 2, 1 )
  554. *
  555. ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
  556. *
  557. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  558. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  559. VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
  560. *
  561. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
  562. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
  563. VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
  564. *
  565. CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
  566. $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
  567. $ ZERO, X, 2, SCALOC, XNORM, IERR )
  568. IF( IERR.NE.0 )
  569. $ INFO = 1
  570. *
  571. IF( SCALOC.NE.ONE ) THEN
  572. DO 90 J = 1, N
  573. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  574. 90 CONTINUE
  575. SCALE = SCALE*SCALOC
  576. END IF
  577. C( K1, L1 ) = X( 1, 1 )
  578. C( K1, L2 ) = X( 2, 1 )
  579. *
  580. ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
  581. *
  582. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  583. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
  584. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  585. *
  586. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
  587. SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
  588. VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
  589. *
  590. SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
  591. SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
  592. VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  593. *
  594. SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
  595. SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
  596. VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
  597. *
  598. CALL DLASY2( .TRUE., .FALSE., ISGN, 2, 2, A( K1, K1 ),
  599. $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
  600. $ 2, XNORM, IERR )
  601. IF( IERR.NE.0 )
  602. $ INFO = 1
  603. *
  604. IF( SCALOC.NE.ONE ) THEN
  605. DO 100 J = 1, N
  606. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  607. 100 CONTINUE
  608. SCALE = SCALE*SCALOC
  609. END IF
  610. C( K1, L1 ) = X( 1, 1 )
  611. C( K1, L2 ) = X( 1, 2 )
  612. C( K2, L1 ) = X( 2, 1 )
  613. C( K2, L2 ) = X( 2, 2 )
  614. END IF
  615. *
  616. 110 CONTINUE
  617. 120 CONTINUE
  618. *
  619. ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
  620. *
  621. * Solve A**T*X + ISGN*X*B**T = scale*C.
  622. *
  623. * The (K,L)th block of X is determined starting from
  624. * top-right corner column by column by
  625. *
  626. * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
  627. *
  628. * Where
  629. * K-1 N
  630. * R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
  631. * I=1 J=L+1
  632. *
  633. * Start column loop (index = L)
  634. * L1 (L2): column index of the first (last) row of X(K,L)
  635. *
  636. LNEXT = N
  637. DO 180 L = N, 1, -1
  638. IF( L.GT.LNEXT )
  639. $ GO TO 180
  640. IF( L.EQ.1 ) THEN
  641. L1 = L
  642. L2 = L
  643. ELSE
  644. IF( B( L, L-1 ).NE.ZERO ) THEN
  645. L1 = L - 1
  646. L2 = L
  647. LNEXT = L - 2
  648. ELSE
  649. L1 = L
  650. L2 = L
  651. LNEXT = L - 1
  652. END IF
  653. END IF
  654. *
  655. * Start row loop (index = K)
  656. * K1 (K2): row index of the first (last) row of X(K,L)
  657. *
  658. KNEXT = 1
  659. DO 170 K = 1, M
  660. IF( K.LT.KNEXT )
  661. $ GO TO 170
  662. IF( K.EQ.M ) THEN
  663. K1 = K
  664. K2 = K
  665. ELSE
  666. IF( A( K+1, K ).NE.ZERO ) THEN
  667. K1 = K
  668. K2 = K + 1
  669. KNEXT = K + 2
  670. ELSE
  671. K1 = K
  672. K2 = K
  673. KNEXT = K + 1
  674. END IF
  675. END IF
  676. *
  677. IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
  678. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  679. SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
  680. $ B( L1, MIN( L1+1, N ) ), LDB )
  681. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  682. SCALOC = ONE
  683. *
  684. A11 = A( K1, K1 ) + SGN*B( L1, L1 )
  685. DA11 = ABS( A11 )
  686. IF( DA11.LE.SMIN ) THEN
  687. A11 = SMIN
  688. DA11 = SMIN
  689. INFO = 1
  690. END IF
  691. DB = ABS( VEC( 1, 1 ) )
  692. IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  693. IF( DB.GT.BIGNUM*DA11 )
  694. $ SCALOC = ONE / DB
  695. END IF
  696. X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
  697. *
  698. IF( SCALOC.NE.ONE ) THEN
  699. DO 130 J = 1, N
  700. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  701. 130 CONTINUE
  702. SCALE = SCALE*SCALOC
  703. END IF
  704. C( K1, L1 ) = X( 1, 1 )
  705. *
  706. ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
  707. *
  708. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  709. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  710. $ B( L1, MIN( L2+1, N ) ), LDB )
  711. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  712. *
  713. SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
  714. SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  715. $ B( L1, MIN( L2+1, N ) ), LDB )
  716. VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  717. *
  718. CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
  719. $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
  720. $ ZERO, X, 2, SCALOC, XNORM, IERR )
  721. IF( IERR.NE.0 )
  722. $ INFO = 1
  723. *
  724. IF( SCALOC.NE.ONE ) THEN
  725. DO 140 J = 1, N
  726. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  727. 140 CONTINUE
  728. SCALE = SCALE*SCALOC
  729. END IF
  730. C( K1, L1 ) = X( 1, 1 )
  731. C( K2, L1 ) = X( 2, 1 )
  732. *
  733. ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
  734. *
  735. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  736. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  737. $ B( L1, MIN( L2+1, N ) ), LDB )
  738. VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
  739. *
  740. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
  741. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  742. $ B( L2, MIN( L2+1, N ) ), LDB )
  743. VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
  744. *
  745. CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
  746. $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
  747. $ ZERO, X, 2, SCALOC, XNORM, IERR )
  748. IF( IERR.NE.0 )
  749. $ INFO = 1
  750. *
  751. IF( SCALOC.NE.ONE ) THEN
  752. DO 150 J = 1, N
  753. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  754. 150 CONTINUE
  755. SCALE = SCALE*SCALOC
  756. END IF
  757. C( K1, L1 ) = X( 1, 1 )
  758. C( K1, L2 ) = X( 2, 1 )
  759. *
  760. ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
  761. *
  762. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
  763. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  764. $ B( L1, MIN( L2+1, N ) ), LDB )
  765. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  766. *
  767. SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
  768. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  769. $ B( L2, MIN( L2+1, N ) ), LDB )
  770. VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
  771. *
  772. SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
  773. SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  774. $ B( L1, MIN( L2+1, N ) ), LDB )
  775. VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  776. *
  777. SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
  778. SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  779. $ B( L2, MIN( L2+1, N ) ), LDB )
  780. VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
  781. *
  782. CALL DLASY2( .TRUE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
  783. $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
  784. $ 2, XNORM, IERR )
  785. IF( IERR.NE.0 )
  786. $ INFO = 1
  787. *
  788. IF( SCALOC.NE.ONE ) THEN
  789. DO 160 J = 1, N
  790. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  791. 160 CONTINUE
  792. SCALE = SCALE*SCALOC
  793. END IF
  794. C( K1, L1 ) = X( 1, 1 )
  795. C( K1, L2 ) = X( 1, 2 )
  796. C( K2, L1 ) = X( 2, 1 )
  797. C( K2, L2 ) = X( 2, 2 )
  798. END IF
  799. *
  800. 170 CONTINUE
  801. 180 CONTINUE
  802. *
  803. ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
  804. *
  805. * Solve A*X + ISGN*X*B**T = scale*C.
  806. *
  807. * The (K,L)th block of X is determined starting from
  808. * bottom-right corner column by column by
  809. *
  810. * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
  811. *
  812. * Where
  813. * M N
  814. * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
  815. * I=K+1 J=L+1
  816. *
  817. * Start column loop (index = L)
  818. * L1 (L2): column index of the first (last) row of X(K,L)
  819. *
  820. LNEXT = N
  821. DO 240 L = N, 1, -1
  822. IF( L.GT.LNEXT )
  823. $ GO TO 240
  824. IF( L.EQ.1 ) THEN
  825. L1 = L
  826. L2 = L
  827. ELSE
  828. IF( B( L, L-1 ).NE.ZERO ) THEN
  829. L1 = L - 1
  830. L2 = L
  831. LNEXT = L - 2
  832. ELSE
  833. L1 = L
  834. L2 = L
  835. LNEXT = L - 1
  836. END IF
  837. END IF
  838. *
  839. * Start row loop (index = K)
  840. * K1 (K2): row index of the first (last) row of X(K,L)
  841. *
  842. KNEXT = M
  843. DO 230 K = M, 1, -1
  844. IF( K.GT.KNEXT )
  845. $ GO TO 230
  846. IF( K.EQ.1 ) THEN
  847. K1 = K
  848. K2 = K
  849. ELSE
  850. IF( A( K, K-1 ).NE.ZERO ) THEN
  851. K1 = K - 1
  852. K2 = K
  853. KNEXT = K - 2
  854. ELSE
  855. K1 = K
  856. K2 = K
  857. KNEXT = K - 1
  858. END IF
  859. END IF
  860. *
  861. IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
  862. SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  863. $ C( MIN( K1+1, M ), L1 ), 1 )
  864. SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
  865. $ B( L1, MIN( L1+1, N ) ), LDB )
  866. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  867. SCALOC = ONE
  868. *
  869. A11 = A( K1, K1 ) + SGN*B( L1, L1 )
  870. DA11 = ABS( A11 )
  871. IF( DA11.LE.SMIN ) THEN
  872. A11 = SMIN
  873. DA11 = SMIN
  874. INFO = 1
  875. END IF
  876. DB = ABS( VEC( 1, 1 ) )
  877. IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  878. IF( DB.GT.BIGNUM*DA11 )
  879. $ SCALOC = ONE / DB
  880. END IF
  881. X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
  882. *
  883. IF( SCALOC.NE.ONE ) THEN
  884. DO 190 J = 1, N
  885. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  886. 190 CONTINUE
  887. SCALE = SCALE*SCALOC
  888. END IF
  889. C( K1, L1 ) = X( 1, 1 )
  890. *
  891. ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
  892. *
  893. SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  894. $ C( MIN( K2+1, M ), L1 ), 1 )
  895. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  896. $ B( L1, MIN( L2+1, N ) ), LDB )
  897. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  898. *
  899. SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  900. $ C( MIN( K2+1, M ), L1 ), 1 )
  901. SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  902. $ B( L1, MIN( L2+1, N ) ), LDB )
  903. VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  904. *
  905. CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
  906. $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
  907. $ ZERO, X, 2, SCALOC, XNORM, IERR )
  908. IF( IERR.NE.0 )
  909. $ INFO = 1
  910. *
  911. IF( SCALOC.NE.ONE ) THEN
  912. DO 200 J = 1, N
  913. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  914. 200 CONTINUE
  915. SCALE = SCALE*SCALOC
  916. END IF
  917. C( K1, L1 ) = X( 1, 1 )
  918. C( K2, L1 ) = X( 2, 1 )
  919. *
  920. ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
  921. *
  922. SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  923. $ C( MIN( K1+1, M ), L1 ), 1 )
  924. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  925. $ B( L1, MIN( L2+1, N ) ), LDB )
  926. VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
  927. *
  928. SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
  929. $ C( MIN( K1+1, M ), L2 ), 1 )
  930. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  931. $ B( L2, MIN( L2+1, N ) ), LDB )
  932. VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
  933. *
  934. CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
  935. $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
  936. $ ZERO, X, 2, SCALOC, XNORM, IERR )
  937. IF( IERR.NE.0 )
  938. $ INFO = 1
  939. *
  940. IF( SCALOC.NE.ONE ) THEN
  941. DO 210 J = 1, N
  942. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  943. 210 CONTINUE
  944. SCALE = SCALE*SCALOC
  945. END IF
  946. C( K1, L1 ) = X( 1, 1 )
  947. C( K1, L2 ) = X( 2, 1 )
  948. *
  949. ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
  950. *
  951. SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  952. $ C( MIN( K2+1, M ), L1 ), 1 )
  953. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  954. $ B( L1, MIN( L2+1, N ) ), LDB )
  955. VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
  956. *
  957. SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
  958. $ C( MIN( K2+1, M ), L2 ), 1 )
  959. SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
  960. $ B( L2, MIN( L2+1, N ) ), LDB )
  961. VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
  962. *
  963. SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  964. $ C( MIN( K2+1, M ), L1 ), 1 )
  965. SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  966. $ B( L1, MIN( L2+1, N ) ), LDB )
  967. VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
  968. *
  969. SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
  970. $ C( MIN( K2+1, M ), L2 ), 1 )
  971. SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
  972. $ B( L2, MIN( L2+1, N ) ), LDB )
  973. VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
  974. *
  975. CALL DLASY2( .FALSE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
  976. $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
  977. $ 2, XNORM, IERR )
  978. IF( IERR.NE.0 )
  979. $ INFO = 1
  980. *
  981. IF( SCALOC.NE.ONE ) THEN
  982. DO 220 J = 1, N
  983. CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
  984. 220 CONTINUE
  985. SCALE = SCALE*SCALOC
  986. END IF
  987. C( K1, L1 ) = X( 1, 1 )
  988. C( K1, L2 ) = X( 1, 2 )
  989. C( K2, L1 ) = X( 2, 1 )
  990. C( K2, L2 ) = X( 2, 2 )
  991. END IF
  992. *
  993. 230 CONTINUE
  994. 240 CONTINUE
  995. *
  996. END IF
  997. *
  998. RETURN
  999. *
  1000. * End of DTRSYL
  1001. *
  1002. END
  1003.  
  1004.  

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