You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

zgebal.f 11 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399
  1. *> \brief \b ZGEBAL
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download ZGEBAL + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgebal.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgebal.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgebal.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE ZGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
  22. *
  23. * .. Scalar Arguments ..
  24. * CHARACTER JOB
  25. * INTEGER IHI, ILO, INFO, LDA, N
  26. * ..
  27. * .. Array Arguments ..
  28. * DOUBLE PRECISION SCALE( * )
  29. * COMPLEX*16 A( LDA, * )
  30. * ..
  31. *
  32. *
  33. *> \par Purpose:
  34. * =============
  35. *>
  36. *> \verbatim
  37. *>
  38. *> ZGEBAL balances a general complex matrix A. This involves, first,
  39. *> permuting A by a similarity transformation to isolate eigenvalues
  40. *> in the first 1 to ILO-1 and last IHI+1 to N elements on the
  41. *> diagonal; and second, applying a diagonal similarity transformation
  42. *> to rows and columns ILO to IHI to make the rows and columns as
  43. *> close in norm as possible. Both steps are optional.
  44. *>
  45. *> Balancing may reduce the 1-norm of the matrix, and improve the
  46. *> accuracy of the computed eigenvalues and/or eigenvectors.
  47. *> \endverbatim
  48. *
  49. * Arguments:
  50. * ==========
  51. *
  52. *> \param[in] JOB
  53. *> \verbatim
  54. *> JOB is CHARACTER*1
  55. *> Specifies the operations to be performed on A:
  56. *> = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0
  57. *> for i = 1,...,N;
  58. *> = 'P': permute only;
  59. *> = 'S': scale only;
  60. *> = 'B': both permute and scale.
  61. *> \endverbatim
  62. *>
  63. *> \param[in] N
  64. *> \verbatim
  65. *> N is INTEGER
  66. *> The order of the matrix A. N >= 0.
  67. *> \endverbatim
  68. *>
  69. *> \param[in,out] A
  70. *> \verbatim
  71. *> A is COMPLEX*16 array, dimension (LDA,N)
  72. *> On entry, the input matrix A.
  73. *> On exit, A is overwritten by the balanced matrix.
  74. *> If JOB = 'N', A is not referenced.
  75. *> See Further Details.
  76. *> \endverbatim
  77. *>
  78. *> \param[in] LDA
  79. *> \verbatim
  80. *> LDA is INTEGER
  81. *> The leading dimension of the array A. LDA >= max(1,N).
  82. *> \endverbatim
  83. *>
  84. *> \param[out] ILO
  85. *> \verbatim
  86. *> \endverbatim
  87. *>
  88. *> \param[out] IHI
  89. *> \verbatim
  90. *> ILO and IHI are set to INTEGER such that on exit
  91. *> A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
  92. *> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
  93. *> \endverbatim
  94. *>
  95. *> \param[out] SCALE
  96. *> \verbatim
  97. *> SCALE is DOUBLE PRECISION array, dimension (N)
  98. *> Details of the permutations and scaling factors applied to
  99. *> A. If P(j) is the index of the row and column interchanged
  100. *> with row and column j and D(j) is the scaling factor
  101. *> applied to row and column j, then
  102. *> SCALE(j) = P(j) for j = 1,...,ILO-1
  103. *> = D(j) for j = ILO,...,IHI
  104. *> = P(j) for j = IHI+1,...,N.
  105. *> The order in which the interchanges are made is N to IHI+1,
  106. *> then 1 to ILO-1.
  107. *> \endverbatim
  108. *>
  109. *> \param[out] INFO
  110. *> \verbatim
  111. *> INFO is INTEGER
  112. *> = 0: successful exit.
  113. *> < 0: if INFO = -i, the i-th argument had an illegal value.
  114. *> \endverbatim
  115. *
  116. * Authors:
  117. * ========
  118. *
  119. *> \author Univ. of Tennessee
  120. *> \author Univ. of California Berkeley
  121. *> \author Univ. of Colorado Denver
  122. *> \author NAG Ltd.
  123. *
  124. *> \date November 2015
  125. *
  126. *> \ingroup complex16GEcomputational
  127. *
  128. *> \par Further Details:
  129. * =====================
  130. *>
  131. *> \verbatim
  132. *>
  133. *> The permutations consist of row and column interchanges which put
  134. *> the matrix in the form
  135. *>
  136. *> ( T1 X Y )
  137. *> P A P = ( 0 B Z )
  138. *> ( 0 0 T2 )
  139. *>
  140. *> where T1 and T2 are upper triangular matrices whose eigenvalues lie
  141. *> along the diagonal. The column indices ILO and IHI mark the starting
  142. *> and ending columns of the submatrix B. Balancing consists of applying
  143. *> a diagonal similarity transformation inv(D) * B * D to make the
  144. *> 1-norms of each row of B and its corresponding column nearly equal.
  145. *> The output matrix is
  146. *>
  147. *> ( T1 X*D Y )
  148. *> ( 0 inv(D)*B*D inv(D)*Z ).
  149. *> ( 0 0 T2 )
  150. *>
  151. *> Information about the permutations P and the diagonal matrix D is
  152. *> returned in the vector SCALE.
  153. *>
  154. *> This subroutine is based on the EISPACK routine CBAL.
  155. *>
  156. *> Modified by Tzu-Yi Chen, Computer Science Division, University of
  157. *> California at Berkeley, USA
  158. *> \endverbatim
  159. *>
  160. * =====================================================================
  161. SUBROUTINE ZGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
  162. *
  163. * -- LAPACK computational routine (version 3.6.0) --
  164. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  165. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  166. * November 2015
  167. *
  168. * .. Scalar Arguments ..
  169. CHARACTER JOB
  170. INTEGER IHI, ILO, INFO, LDA, N
  171. * ..
  172. * .. Array Arguments ..
  173. DOUBLE PRECISION SCALE( * )
  174. COMPLEX*16 A( LDA, * )
  175. * ..
  176. *
  177. * =====================================================================
  178. *
  179. * .. Parameters ..
  180. DOUBLE PRECISION ZERO, ONE
  181. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  182. DOUBLE PRECISION SCLFAC
  183. PARAMETER ( SCLFAC = 2.0D+0 )
  184. DOUBLE PRECISION FACTOR
  185. PARAMETER ( FACTOR = 0.95D+0 )
  186. * ..
  187. * .. Local Scalars ..
  188. LOGICAL NOCONV
  189. INTEGER I, ICA, IEXC, IRA, J, K, L, M
  190. DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
  191. $ SFMIN2
  192. COMPLEX*16 CDUM
  193. * ..
  194. * .. External Functions ..
  195. LOGICAL DISNAN, LSAME
  196. INTEGER IZAMAX
  197. DOUBLE PRECISION DLAMCH, DZNRM2
  198. EXTERNAL DISNAN, LSAME, IZAMAX, DLAMCH, DZNRM2
  199. * ..
  200. * .. External Subroutines ..
  201. EXTERNAL XERBLA, ZDSCAL, ZSWAP
  202. * ..
  203. * .. Intrinsic Functions ..
  204. INTRINSIC ABS, DBLE, DIMAG, MAX, MIN
  205. *
  206. * Test the input parameters
  207. *
  208. INFO = 0
  209. IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
  210. $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
  211. INFO = -1
  212. ELSE IF( N.LT.0 ) THEN
  213. INFO = -2
  214. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  215. INFO = -4
  216. END IF
  217. IF( INFO.NE.0 ) THEN
  218. CALL XERBLA( 'ZGEBAL', -INFO )
  219. RETURN
  220. END IF
  221. *
  222. K = 1
  223. L = N
  224. *
  225. IF( N.EQ.0 )
  226. $ GO TO 210
  227. *
  228. IF( LSAME( JOB, 'N' ) ) THEN
  229. DO 10 I = 1, N
  230. SCALE( I ) = ONE
  231. 10 CONTINUE
  232. GO TO 210
  233. END IF
  234. *
  235. IF( LSAME( JOB, 'S' ) )
  236. $ GO TO 120
  237. *
  238. * Permutation to isolate eigenvalues if possible
  239. *
  240. GO TO 50
  241. *
  242. * Row and column exchange.
  243. *
  244. 20 CONTINUE
  245. SCALE( M ) = J
  246. IF( J.EQ.M )
  247. $ GO TO 30
  248. *
  249. CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
  250. CALL ZSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
  251. *
  252. 30 CONTINUE
  253. GO TO ( 40, 80 )IEXC
  254. *
  255. * Search for rows isolating an eigenvalue and push them down.
  256. *
  257. 40 CONTINUE
  258. IF( L.EQ.1 )
  259. $ GO TO 210
  260. L = L - 1
  261. *
  262. 50 CONTINUE
  263. DO 70 J = L, 1, -1
  264. *
  265. DO 60 I = 1, L
  266. IF( I.EQ.J )
  267. $ GO TO 60
  268. IF( DBLE( A( J, I ) ).NE.ZERO .OR. DIMAG( A( J, I ) ).NE.
  269. $ ZERO )GO TO 70
  270. 60 CONTINUE
  271. *
  272. M = L
  273. IEXC = 1
  274. GO TO 20
  275. 70 CONTINUE
  276. *
  277. GO TO 90
  278. *
  279. * Search for columns isolating an eigenvalue and push them left.
  280. *
  281. 80 CONTINUE
  282. K = K + 1
  283. *
  284. 90 CONTINUE
  285. DO 110 J = K, L
  286. *
  287. DO 100 I = K, L
  288. IF( I.EQ.J )
  289. $ GO TO 100
  290. IF( DBLE( A( I, J ) ).NE.ZERO .OR. DIMAG( A( I, J ) ).NE.
  291. $ ZERO )GO TO 110
  292. 100 CONTINUE
  293. *
  294. M = K
  295. IEXC = 2
  296. GO TO 20
  297. 110 CONTINUE
  298. *
  299. 120 CONTINUE
  300. DO 130 I = K, L
  301. SCALE( I ) = ONE
  302. 130 CONTINUE
  303. *
  304. IF( LSAME( JOB, 'P' ) )
  305. $ GO TO 210
  306. *
  307. * Balance the submatrix in rows K to L.
  308. *
  309. * Iterative loop for norm reduction
  310. *
  311. SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
  312. SFMAX1 = ONE / SFMIN1
  313. SFMIN2 = SFMIN1*SCLFAC
  314. SFMAX2 = ONE / SFMIN2
  315. 140 CONTINUE
  316. NOCONV = .FALSE.
  317. *
  318. DO 200 I = K, L
  319. *
  320. C = DZNRM2( L-K+1, A( K, I ), 1 )
  321. R = DZNRM2( L-K+1, A( I, K ), LDA )
  322. ICA = IZAMAX( L, A( 1, I ), 1 )
  323. CA = ABS( A( ICA, I ) )
  324. IRA = IZAMAX( N-K+1, A( I, K ), LDA )
  325. RA = ABS( A( I, IRA+K-1 ) )
  326. *
  327. * Guard against zero C or R due to underflow.
  328. *
  329. IF( C.EQ.ZERO .OR. R.EQ.ZERO )
  330. $ GO TO 200
  331. G = R / SCLFAC
  332. F = ONE
  333. S = C + R
  334. 160 CONTINUE
  335. IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
  336. $ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
  337. IF( DISNAN( C+F+CA+R+G+RA ) ) THEN
  338. *
  339. * Exit if NaN to avoid infinite loop
  340. *
  341. INFO = -3
  342. CALL XERBLA( 'ZGEBAL', -INFO )
  343. RETURN
  344. END IF
  345. F = F*SCLFAC
  346. C = C*SCLFAC
  347. CA = CA*SCLFAC
  348. R = R / SCLFAC
  349. G = G / SCLFAC
  350. RA = RA / SCLFAC
  351. GO TO 160
  352. *
  353. 170 CONTINUE
  354. G = C / SCLFAC
  355. 180 CONTINUE
  356. IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
  357. $ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
  358. F = F / SCLFAC
  359. C = C / SCLFAC
  360. G = G / SCLFAC
  361. CA = CA / SCLFAC
  362. R = R*SCLFAC
  363. RA = RA*SCLFAC
  364. GO TO 180
  365. *
  366. * Now balance.
  367. *
  368. 190 CONTINUE
  369. IF( ( C+R ).GE.FACTOR*S )
  370. $ GO TO 200
  371. IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
  372. IF( F*SCALE( I ).LE.SFMIN1 )
  373. $ GO TO 200
  374. END IF
  375. IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
  376. IF( SCALE( I ).GE.SFMAX1 / F )
  377. $ GO TO 200
  378. END IF
  379. G = ONE / F
  380. SCALE( I ) = SCALE( I )*F
  381. NOCONV = .TRUE.
  382. *
  383. CALL ZDSCAL( N-K+1, G, A( I, K ), LDA )
  384. CALL ZDSCAL( L, F, A( 1, I ), 1 )
  385. *
  386. 200 CONTINUE
  387. *
  388. IF( NOCONV )
  389. $ GO TO 140
  390. *
  391. 210 CONTINUE
  392. ILO = K
  393. IHI = L
  394. *
  395. RETURN
  396. *
  397. * End of ZGEBAL
  398. *
  399. END