Télécharger dgebal.eso

Retour à la liste

Numérotation des lignes :

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

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