6ad934d3d5f570ea2200c3eed846c1292848bf47
[platform/upstream/lapack.git] / SRC / zggev.f
1 *> \brief <b> ZGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at 
6 *            http://www.netlib.org/lapack/explore-html/ 
7 *
8 *> \htmlonly
9 *> Download ZGGEV + dependencies 
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggev.f"> 
11 *> [TGZ]</a> 
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggev.f"> 
13 *> [ZIP]</a> 
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggev.f"> 
15 *> [TXT]</a>
16 *> \endhtmlonly 
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE ZGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
22 *                         VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
23
24 *       .. Scalar Arguments ..
25 *       CHARACTER          JOBVL, JOBVR
26 *       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
27 *       ..
28 *       .. Array Arguments ..
29 *       DOUBLE PRECISION   RWORK( * )
30 *       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
31 *      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
32 *      $                   WORK( * )
33 *       ..
34 *  
35 *
36 *> \par Purpose:
37 *  =============
38 *>
39 *> \verbatim
40 *>
41 *> ZGGEV computes for a pair of N-by-N complex nonsymmetric matrices
42 *> (A,B), the generalized eigenvalues, and optionally, the left and/or
43 *> right generalized eigenvectors.
44 *>
45 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
46 *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
47 *> singular. It is usually represented as the pair (alpha,beta), as
48 *> there is a reasonable interpretation for beta=0, and even for both
49 *> being zero.
50 *>
51 *> The right generalized eigenvector v(j) corresponding to the
52 *> generalized eigenvalue lambda(j) of (A,B) satisfies
53 *>
54 *>              A * v(j) = lambda(j) * B * v(j).
55 *>
56 *> The left generalized eigenvector u(j) corresponding to the
57 *> generalized eigenvalues lambda(j) of (A,B) satisfies
58 *>
59 *>              u(j)**H * A = lambda(j) * u(j)**H * B
60 *>
61 *> where u(j)**H is the conjugate-transpose of u(j).
62 *> \endverbatim
63 *
64 *  Arguments:
65 *  ==========
66 *
67 *> \param[in] JOBVL
68 *> \verbatim
69 *>          JOBVL is CHARACTER*1
70 *>          = 'N':  do not compute the left generalized eigenvectors;
71 *>          = 'V':  compute the left generalized eigenvectors.
72 *> \endverbatim
73 *>
74 *> \param[in] JOBVR
75 *> \verbatim
76 *>          JOBVR is CHARACTER*1
77 *>          = 'N':  do not compute the right generalized eigenvectors;
78 *>          = 'V':  compute the right generalized eigenvectors.
79 *> \endverbatim
80 *>
81 *> \param[in] N
82 *> \verbatim
83 *>          N is INTEGER
84 *>          The order of the matrices A, B, VL, and VR.  N >= 0.
85 *> \endverbatim
86 *>
87 *> \param[in,out] A
88 *> \verbatim
89 *>          A is COMPLEX*16 array, dimension (LDA, N)
90 *>          On entry, the matrix A in the pair (A,B).
91 *>          On exit, A has been overwritten.
92 *> \endverbatim
93 *>
94 *> \param[in] LDA
95 *> \verbatim
96 *>          LDA is INTEGER
97 *>          The leading dimension of A.  LDA >= max(1,N).
98 *> \endverbatim
99 *>
100 *> \param[in,out] B
101 *> \verbatim
102 *>          B is COMPLEX*16 array, dimension (LDB, N)
103 *>          On entry, the matrix B in the pair (A,B).
104 *>          On exit, B has been overwritten.
105 *> \endverbatim
106 *>
107 *> \param[in] LDB
108 *> \verbatim
109 *>          LDB is INTEGER
110 *>          The leading dimension of B.  LDB >= max(1,N).
111 *> \endverbatim
112 *>
113 *> \param[out] ALPHA
114 *> \verbatim
115 *>          ALPHA is COMPLEX*16 array, dimension (N)
116 *> \endverbatim
117 *>
118 *> \param[out] BETA
119 *> \verbatim
120 *>          BETA is COMPLEX*16 array, dimension (N)
121 *>          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
122 *>          generalized eigenvalues.
123 *>
124 *>          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
125 *>          underflow, and BETA(j) may even be zero.  Thus, the user
126 *>          should avoid naively computing the ratio alpha/beta.
127 *>          However, ALPHA will be always less than and usually
128 *>          comparable with norm(A) in magnitude, and BETA always less
129 *>          than and usually comparable with norm(B).
130 *> \endverbatim
131 *>
132 *> \param[out] VL
133 *> \verbatim
134 *>          VL is COMPLEX*16 array, dimension (LDVL,N)
135 *>          If JOBVL = 'V', the left generalized eigenvectors u(j) are
136 *>          stored one after another in the columns of VL, in the same
137 *>          order as their eigenvalues.
138 *>          Each eigenvector is scaled so the largest component has
139 *>          abs(real part) + abs(imag. part) = 1.
140 *>          Not referenced if JOBVL = 'N'.
141 *> \endverbatim
142 *>
143 *> \param[in] LDVL
144 *> \verbatim
145 *>          LDVL is INTEGER
146 *>          The leading dimension of the matrix VL. LDVL >= 1, and
147 *>          if JOBVL = 'V', LDVL >= N.
148 *> \endverbatim
149 *>
150 *> \param[out] VR
151 *> \verbatim
152 *>          VR is COMPLEX*16 array, dimension (LDVR,N)
153 *>          If JOBVR = 'V', the right generalized eigenvectors v(j) are
154 *>          stored one after another in the columns of VR, in the same
155 *>          order as their eigenvalues.
156 *>          Each eigenvector is scaled so the largest component has
157 *>          abs(real part) + abs(imag. part) = 1.
158 *>          Not referenced if JOBVR = 'N'.
159 *> \endverbatim
160 *>
161 *> \param[in] LDVR
162 *> \verbatim
163 *>          LDVR is INTEGER
164 *>          The leading dimension of the matrix VR. LDVR >= 1, and
165 *>          if JOBVR = 'V', LDVR >= N.
166 *> \endverbatim
167 *>
168 *> \param[out] WORK
169 *> \verbatim
170 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
171 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
172 *> \endverbatim
173 *>
174 *> \param[in] LWORK
175 *> \verbatim
176 *>          LWORK is INTEGER
177 *>          The dimension of the array WORK.  LWORK >= max(1,2*N).
178 *>          For good performance, LWORK must generally be larger.
179 *>
180 *>          If LWORK = -1, then a workspace query is assumed; the routine
181 *>          only calculates the optimal size of the WORK array, returns
182 *>          this value as the first entry of the WORK array, and no error
183 *>          message related to LWORK is issued by XERBLA.
184 *> \endverbatim
185 *>
186 *> \param[out] RWORK
187 *> \verbatim
188 *>          RWORK is DOUBLE PRECISION array, dimension (8*N)
189 *> \endverbatim
190 *>
191 *> \param[out] INFO
192 *> \verbatim
193 *>          INFO is INTEGER
194 *>          = 0:  successful exit
195 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
196 *>          =1,...,N:
197 *>                The QZ iteration failed.  No eigenvectors have been
198 *>                calculated, but ALPHA(j) and BETA(j) should be
199 *>                correct for j=INFO+1,...,N.
200 *>          > N:  =N+1: other then QZ iteration failed in DHGEQZ,
201 *>                =N+2: error return from DTGEVC.
202 *> \endverbatim
203 *
204 *  Authors:
205 *  ========
206 *
207 *> \author Univ. of Tennessee 
208 *> \author Univ. of California Berkeley 
209 *> \author Univ. of Colorado Denver 
210 *> \author NAG Ltd. 
211 *
212 *> \date April 2012
213 *
214 *> \ingroup complex16GEeigen
215 *
216 *  =====================================================================
217       SUBROUTINE ZGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
218      $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
219 *
220 *  -- LAPACK driver routine (version 3.4.1) --
221 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
222 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
223 *     April 2012
224 *
225 *     .. Scalar Arguments ..
226       CHARACTER          JOBVL, JOBVR
227       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
228 *     ..
229 *     .. Array Arguments ..
230       DOUBLE PRECISION   RWORK( * )
231       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
232      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
233      $                   WORK( * )
234 *     ..
235 *
236 *  =====================================================================
237 *
238 *     .. Parameters ..
239       DOUBLE PRECISION   ZERO, ONE
240       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
241       COMPLEX*16         CZERO, CONE
242       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
243      $                   CONE = ( 1.0D0, 0.0D0 ) )
244 *     ..
245 *     .. Local Scalars ..
246       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
247       CHARACTER          CHTEMP
248       INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
249      $                   IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
250      $                   LWKMIN, LWKOPT
251       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
252      $                   SMLNUM, TEMP
253       COMPLEX*16         X
254 *     ..
255 *     .. Local Arrays ..
256       LOGICAL            LDUMMA( 1 )
257 *     ..
258 *     .. External Subroutines ..
259       EXTERNAL           DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
260      $                   ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR,
261      $                   ZUNMQR
262 *     ..
263 *     .. External Functions ..
264       LOGICAL            LSAME
265       INTEGER            ILAENV
266       DOUBLE PRECISION   DLAMCH, ZLANGE
267       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
268 *     ..
269 *     .. Intrinsic Functions ..
270       INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
271 *     ..
272 *     .. Statement Functions ..
273       DOUBLE PRECISION   ABS1
274 *     ..
275 *     .. Statement Function definitions ..
276       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
277 *     ..
278 *     .. Executable Statements ..
279 *
280 *     Decode the input arguments
281 *
282       IF( LSAME( JOBVL, 'N' ) ) THEN
283          IJOBVL = 1
284          ILVL = .FALSE.
285       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
286          IJOBVL = 2
287          ILVL = .TRUE.
288       ELSE
289          IJOBVL = -1
290          ILVL = .FALSE.
291       END IF
292 *
293       IF( LSAME( JOBVR, 'N' ) ) THEN
294          IJOBVR = 1
295          ILVR = .FALSE.
296       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
297          IJOBVR = 2
298          ILVR = .TRUE.
299       ELSE
300          IJOBVR = -1
301          ILVR = .FALSE.
302       END IF
303       ILV = ILVL .OR. ILVR
304 *
305 *     Test the input arguments
306 *
307       INFO = 0
308       LQUERY = ( LWORK.EQ.-1 )
309       IF( IJOBVL.LE.0 ) THEN
310          INFO = -1
311       ELSE IF( IJOBVR.LE.0 ) THEN
312          INFO = -2
313       ELSE IF( N.LT.0 ) THEN
314          INFO = -3
315       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
316          INFO = -5
317       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
318          INFO = -7
319       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
320          INFO = -11
321       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
322          INFO = -13
323       END IF
324 *
325 *     Compute workspace
326 *      (Note: Comments in the code beginning "Workspace:" describe the
327 *       minimal amount of workspace needed at that point in the code,
328 *       as well as the preferred amount for good performance.
329 *       NB refers to the optimal block size for the immediately
330 *       following subroutine, as returned by ILAENV. The workspace is
331 *       computed assuming ILO = 1 and IHI = N, the worst case.)
332 *
333       IF( INFO.EQ.0 ) THEN
334          LWKMIN = MAX( 1, 2*N )
335          LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
336          LWKOPT = MAX( LWKOPT, N +
337      $                 N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
338          IF( ILVL ) THEN
339             LWKOPT = MAX( LWKOPT, N +
340      $                    N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) )
341          END IF
342          WORK( 1 ) = LWKOPT
343 *
344          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
345      $      INFO = -15
346       END IF
347 *
348       IF( INFO.NE.0 ) THEN
349          CALL XERBLA( 'ZGGEV ', -INFO )
350          RETURN
351       ELSE IF( LQUERY ) THEN
352          RETURN
353       END IF
354 *
355 *     Quick return if possible
356 *
357       IF( N.EQ.0 )
358      $   RETURN
359 *
360 *     Get machine constants
361 *
362       EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
363       SMLNUM = DLAMCH( 'S' )
364       BIGNUM = ONE / SMLNUM
365       CALL DLABAD( SMLNUM, BIGNUM )
366       SMLNUM = SQRT( SMLNUM ) / EPS
367       BIGNUM = ONE / SMLNUM
368 *
369 *     Scale A if max element outside range [SMLNUM,BIGNUM]
370 *
371       ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
372       ILASCL = .FALSE.
373       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
374          ANRMTO = SMLNUM
375          ILASCL = .TRUE.
376       ELSE IF( ANRM.GT.BIGNUM ) THEN
377          ANRMTO = BIGNUM
378          ILASCL = .TRUE.
379       END IF
380       IF( ILASCL )
381      $   CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
382 *
383 *     Scale B if max element outside range [SMLNUM,BIGNUM]
384 *
385       BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
386       ILBSCL = .FALSE.
387       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
388          BNRMTO = SMLNUM
389          ILBSCL = .TRUE.
390       ELSE IF( BNRM.GT.BIGNUM ) THEN
391          BNRMTO = BIGNUM
392          ILBSCL = .TRUE.
393       END IF
394       IF( ILBSCL )
395      $   CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
396 *
397 *     Permute the matrices A, B to isolate eigenvalues if possible
398 *     (Real Workspace: need 6*N)
399 *
400       ILEFT = 1
401       IRIGHT = N + 1
402       IRWRK = IRIGHT + N
403       CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
404      $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
405 *
406 *     Reduce B to triangular form (QR decomposition of B)
407 *     (Complex Workspace: need N, prefer N*NB)
408 *
409       IROWS = IHI + 1 - ILO
410       IF( ILV ) THEN
411          ICOLS = N + 1 - ILO
412       ELSE
413          ICOLS = IROWS
414       END IF
415       ITAU = 1
416       IWRK = ITAU + IROWS
417       CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
418      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
419 *
420 *     Apply the orthogonal transformation to matrix A
421 *     (Complex Workspace: need N, prefer N*NB)
422 *
423       CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
424      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
425      $             LWORK+1-IWRK, IERR )
426 *
427 *     Initialize VL
428 *     (Complex Workspace: need N, prefer N*NB)
429 *
430       IF( ILVL ) THEN
431          CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
432          IF( IROWS.GT.1 ) THEN
433             CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
434      $                   VL( ILO+1, ILO ), LDVL )
435          END IF
436          CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
437      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
438       END IF
439 *
440 *     Initialize VR
441 *
442       IF( ILVR )
443      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
444 *
445 *     Reduce to generalized Hessenberg form
446 *
447       IF( ILV ) THEN
448 *
449 *        Eigenvectors requested -- work on whole matrix.
450 *
451          CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
452      $                LDVL, VR, LDVR, IERR )
453       ELSE
454          CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
455      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
456       END IF
457 *
458 *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
459 *     Schur form and Schur vectors)
460 *     (Complex Workspace: need N)
461 *     (Real Workspace: need N)
462 *
463       IWRK = ITAU
464       IF( ILV ) THEN
465          CHTEMP = 'S'
466       ELSE
467          CHTEMP = 'E'
468       END IF
469       CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
470      $             ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
471      $             LWORK+1-IWRK, RWORK( IRWRK ), IERR )
472       IF( IERR.NE.0 ) THEN
473          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
474             INFO = IERR
475          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
476             INFO = IERR - N
477          ELSE
478             INFO = N + 1
479          END IF
480          GO TO 70
481       END IF
482 *
483 *     Compute Eigenvectors
484 *     (Real Workspace: need 2*N)
485 *     (Complex Workspace: need 2*N)
486 *
487       IF( ILV ) THEN
488          IF( ILVL ) THEN
489             IF( ILVR ) THEN
490                CHTEMP = 'B'
491             ELSE
492                CHTEMP = 'L'
493             END IF
494          ELSE
495             CHTEMP = 'R'
496          END IF
497 *
498          CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
499      $                VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ),
500      $                IERR )
501          IF( IERR.NE.0 ) THEN
502             INFO = N + 2
503             GO TO 70
504          END IF
505 *
506 *        Undo balancing on VL and VR and normalization
507 *        (Workspace: none needed)
508 *
509          IF( ILVL ) THEN
510             CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
511      $                   RWORK( IRIGHT ), N, VL, LDVL, IERR )
512             DO 30 JC = 1, N
513                TEMP = ZERO
514                DO 10 JR = 1, N
515                   TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
516    10          CONTINUE
517                IF( TEMP.LT.SMLNUM )
518      $            GO TO 30
519                TEMP = ONE / TEMP
520                DO 20 JR = 1, N
521                   VL( JR, JC ) = VL( JR, JC )*TEMP
522    20          CONTINUE
523    30       CONTINUE
524          END IF
525          IF( ILVR ) THEN
526             CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
527      $                   RWORK( IRIGHT ), N, VR, LDVR, IERR )
528             DO 60 JC = 1, N
529                TEMP = ZERO
530                DO 40 JR = 1, N
531                   TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
532    40          CONTINUE
533                IF( TEMP.LT.SMLNUM )
534      $            GO TO 60
535                TEMP = ONE / TEMP
536                DO 50 JR = 1, N
537                   VR( JR, JC ) = VR( JR, JC )*TEMP
538    50          CONTINUE
539    60       CONTINUE
540          END IF
541       END IF
542 *
543 *     Undo scaling if necessary
544 *
545    70 CONTINUE
546 *
547       IF( ILASCL )
548      $   CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
549 *
550       IF( ILBSCL )
551      $   CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
552 *
553       WORK( 1 ) = LWKOPT
554       RETURN
555 *
556 *     End of ZGGEV
557 *
558       END