← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/lapack/single/sgeev.f

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lapack/single/sgeev.f
 
1
*> \brief <b> SGEEV 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 SGEEV + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeev.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeev.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeev.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE SGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
 
22
*                         LDVR, WORK, LWORK, INFO )
 
23
* 
 
24
*       .. Scalar Arguments ..
 
25
*       CHARACTER          JOBVL, JOBVR
 
26
*       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
 
27
*       ..
 
28
*       .. Array Arguments ..
 
29
*       REAL               A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
 
30
*      $                   WI( * ), WORK( * ), WR( * )
 
31
*       ..
 
32
*  
 
33
*
 
34
*> \par Purpose:
 
35
*  =============
 
36
*>
 
37
*> \verbatim
 
38
*>
 
39
*> SGEEV computes for an N-by-N real nonsymmetric matrix A, the
 
40
*> eigenvalues and, optionally, the left and/or right eigenvectors.
 
41
*>
 
42
*> The right eigenvector v(j) of A satisfies
 
43
*>                  A * v(j) = lambda(j) * v(j)
 
44
*> where lambda(j) is its eigenvalue.
 
45
*> The left eigenvector u(j) of A satisfies
 
46
*>               u(j)**H * A = lambda(j) * u(j)**H
 
47
*> where u(j)**H denotes the conjugate-transpose of u(j).
 
48
*>
 
49
*> The computed eigenvectors are normalized to have Euclidean norm
 
50
*> equal to 1 and largest component real.
 
51
*> \endverbatim
 
52
*
 
53
*  Arguments:
 
54
*  ==========
 
55
*
 
56
*> \param[in] JOBVL
 
57
*> \verbatim
 
58
*>          JOBVL is CHARACTER*1
 
59
*>          = 'N': left eigenvectors of A are not computed;
 
60
*>          = 'V': left eigenvectors of A are computed.
 
61
*> \endverbatim
 
62
*>
 
63
*> \param[in] JOBVR
 
64
*> \verbatim
 
65
*>          JOBVR is CHARACTER*1
 
66
*>          = 'N': right eigenvectors of A are not computed;
 
67
*>          = 'V': right eigenvectors of A are computed.
 
68
*> \endverbatim
 
69
*>
 
70
*> \param[in] N
 
71
*> \verbatim
 
72
*>          N is INTEGER
 
73
*>          The order of the matrix A. N >= 0.
 
74
*> \endverbatim
 
75
*>
 
76
*> \param[in,out] A
 
77
*> \verbatim
 
78
*>          A is REAL array, dimension (LDA,N)
 
79
*>          On entry, the N-by-N matrix A.
 
80
*>          On exit, A has been overwritten.
 
81
*> \endverbatim
 
82
*>
 
83
*> \param[in] LDA
 
84
*> \verbatim
 
85
*>          LDA is INTEGER
 
86
*>          The leading dimension of the array A.  LDA >= max(1,N).
 
87
*> \endverbatim
 
88
*>
 
89
*> \param[out] WR
 
90
*> \verbatim
 
91
*>          WR is REAL array, dimension (N)
 
92
*> \endverbatim
 
93
*>
 
94
*> \param[out] WI
 
95
*> \verbatim
 
96
*>          WI is REAL array, dimension (N)
 
97
*>          WR and WI contain the real and imaginary parts,
 
98
*>          respectively, of the computed eigenvalues.  Complex
 
99
*>          conjugate pairs of eigenvalues appear consecutively
 
100
*>          with the eigenvalue having the positive imaginary part
 
101
*>          first.
 
102
*> \endverbatim
 
103
*>
 
104
*> \param[out] VL
 
105
*> \verbatim
 
106
*>          VL is REAL array, dimension (LDVL,N)
 
107
*>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
 
108
*>          after another in the columns of VL, in the same order
 
109
*>          as their eigenvalues.
 
110
*>          If JOBVL = 'N', VL is not referenced.
 
111
*>          If the j-th eigenvalue is real, then u(j) = VL(:,j),
 
112
*>          the j-th column of VL.
 
113
*>          If the j-th and (j+1)-st eigenvalues form a complex
 
114
*>          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
 
115
*>          u(j+1) = VL(:,j) - i*VL(:,j+1).
 
116
*> \endverbatim
 
117
*>
 
118
*> \param[in] LDVL
 
119
*> \verbatim
 
120
*>          LDVL is INTEGER
 
121
*>          The leading dimension of the array VL.  LDVL >= 1; if
 
122
*>          JOBVL = 'V', LDVL >= N.
 
123
*> \endverbatim
 
124
*>
 
125
*> \param[out] VR
 
126
*> \verbatim
 
127
*>          VR is REAL array, dimension (LDVR,N)
 
128
*>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
 
129
*>          after another in the columns of VR, in the same order
 
130
*>          as their eigenvalues.
 
131
*>          If JOBVR = 'N', VR is not referenced.
 
132
*>          If the j-th eigenvalue is real, then v(j) = VR(:,j),
 
133
*>          the j-th column of VR.
 
134
*>          If the j-th and (j+1)-st eigenvalues form a complex
 
135
*>          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
 
136
*>          v(j+1) = VR(:,j) - i*VR(:,j+1).
 
137
*> \endverbatim
 
138
*>
 
139
*> \param[in] LDVR
 
140
*> \verbatim
 
141
*>          LDVR is INTEGER
 
142
*>          The leading dimension of the array VR.  LDVR >= 1; if
 
143
*>          JOBVR = 'V', LDVR >= N.
 
144
*> \endverbatim
 
145
*>
 
146
*> \param[out] WORK
 
147
*> \verbatim
 
148
*>          WORK is REAL array, dimension (MAX(1,LWORK))
 
149
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 
150
*> \endverbatim
 
151
*>
 
152
*> \param[in] LWORK
 
153
*> \verbatim
 
154
*>          LWORK is INTEGER
 
155
*>          The dimension of the array WORK.  LWORK >= max(1,3*N), and
 
156
*>          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
 
157
*>          performance, LWORK must generally be larger.
 
158
*>
 
159
*>          If LWORK = -1, then a workspace query is assumed; the routine
 
160
*>          only calculates the optimal size of the WORK array, returns
 
161
*>          this value as the first entry of the WORK array, and no error
 
162
*>          message related to LWORK is issued by XERBLA.
 
163
*> \endverbatim
 
164
*>
 
165
*> \param[out] INFO
 
166
*> \verbatim
 
167
*>          INFO is INTEGER
 
168
*>          = 0:  successful exit
 
169
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 
170
*>          > 0:  if INFO = i, the QR algorithm failed to compute all the
 
171
*>                eigenvalues, and no eigenvectors have been computed;
 
172
*>                elements i+1:N of WR and WI contain eigenvalues which
 
173
*>                have converged.
 
174
*> \endverbatim
 
175
*
 
176
*  Authors:
 
177
*  ========
 
178
*
 
179
*> \author Univ. of Tennessee 
 
180
*> \author Univ. of California Berkeley 
 
181
*> \author Univ. of Colorado Denver 
 
182
*> \author NAG Ltd. 
 
183
*
 
184
*> \date September 2012
 
185
*
 
186
*> \ingroup realGEeigen
 
187
*
 
188
*  =====================================================================
 
189
      SUBROUTINE SGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
 
190
     $                  LDVR, WORK, LWORK, INFO )
 
191
*
 
192
*  -- LAPACK driver routine (version 3.4.2) --
 
193
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
194
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
195
*     September 2012
 
196
*
 
197
*     .. Scalar Arguments ..
 
198
      CHARACTER          JOBVL, JOBVR
 
199
      INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
 
200
*     ..
 
201
*     .. Array Arguments ..
 
202
      REAL               A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
 
203
     $                   WI( * ), WORK( * ), WR( * )
 
204
*     ..
 
205
*
 
206
*  =====================================================================
 
207
*
 
208
*     .. Parameters ..
 
209
      REAL               ZERO, ONE
 
210
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 
211
*     ..
 
212
*     .. Local Scalars ..
 
213
      LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
 
214
      CHARACTER          SIDE
 
215
      INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
 
216
     $                   MAXWRK, MINWRK, NOUT
 
217
      REAL               ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
 
218
     $                   SN
 
219
*     ..
 
220
*     .. Local Arrays ..
 
221
      LOGICAL            SELECT( 1 )
 
222
      REAL               DUM( 1 )
 
223
*     ..
 
224
*     .. External Subroutines ..
 
225
      EXTERNAL           SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD, SLACPY,
 
226
     $                   SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC,
 
227
     $                   XERBLA
 
228
*     ..
 
229
*     .. External Functions ..
 
230
      LOGICAL            LSAME
 
231
      INTEGER            ILAENV, ISAMAX
 
232
      REAL               SLAMCH, SLANGE, SLAPY2, SNRM2
 
233
      EXTERNAL           LSAME, ILAENV, ISAMAX, SLAMCH, SLANGE, SLAPY2,
 
234
     $                   SNRM2
 
235
*     ..
 
236
*     .. Intrinsic Functions ..
 
237
      INTRINSIC          MAX, SQRT
 
238
*     ..
 
239
*     .. Executable Statements ..
 
240
*
 
241
*     Test the input arguments
 
242
*
 
243
      INFO = 0
 
244
      LQUERY = ( LWORK.EQ.-1 )
 
245
      WANTVL = LSAME( JOBVL, 'V' )
 
246
      WANTVR = LSAME( JOBVR, 'V' )
 
247
      IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
 
248
         INFO = -1
 
249
      ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
 
250
         INFO = -2
 
251
      ELSE IF( N.LT.0 ) THEN
 
252
         INFO = -3
 
253
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 
254
         INFO = -5
 
255
      ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
 
256
         INFO = -9
 
257
      ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
 
258
         INFO = -11
 
259
      END IF
 
260
*
 
261
*     Compute workspace
 
262
*      (Note: Comments in the code beginning "Workspace:" describe the
 
263
*       minimal amount of workspace needed at that point in the code,
 
264
*       as well as the preferred amount for good performance.
 
265
*       NB refers to the optimal block size for the immediately
 
266
*       following subroutine, as returned by ILAENV.
 
267
*       HSWORK refers to the workspace preferred by SHSEQR, as
 
268
*       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
 
269
*       the worst case.)
 
270
*
 
271
      IF( INFO.EQ.0 ) THEN
 
272
         IF( N.EQ.0 ) THEN
 
273
            MINWRK = 1
 
274
            MAXWRK = 1
 
275
         ELSE
 
276
            MAXWRK = 2*N + N*ILAENV( 1, 'SGEHRD', ' ', N, 1, N, 0 )
 
277
            IF( WANTVL ) THEN
 
278
               MINWRK = 4*N
 
279
               MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
 
280
     $                       'SORGHR', ' ', N, 1, N, -1 ) )
 
281
               CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
 
282
     $                WORK, -1, INFO )
 
283
               HSWORK = WORK( 1 )
 
284
               MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
 
285
               MAXWRK = MAX( MAXWRK, 4*N )
 
286
            ELSE IF( WANTVR ) THEN
 
287
               MINWRK = 4*N
 
288
               MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
 
289
     $                       'SORGHR', ' ', N, 1, N, -1 ) )
 
290
               CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
 
291
     $                WORK, -1, INFO )
 
292
               HSWORK = WORK( 1 )
 
293
               MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
 
294
               MAXWRK = MAX( MAXWRK, 4*N )
 
295
            ELSE 
 
296
               MINWRK = 3*N
 
297
               CALL SHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
 
298
     $                WORK, -1, INFO )
 
299
               HSWORK = WORK( 1 )
 
300
               MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
 
301
            END IF
 
302
            MAXWRK = MAX( MAXWRK, MINWRK )
 
303
         END IF
 
304
         WORK( 1 ) = MAXWRK
 
305
*
 
306
         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 
307
            INFO = -13
 
308
         END IF
 
309
      END IF
 
310
*
 
311
      IF( INFO.NE.0 ) THEN
 
312
         CALL XERBLA( 'SGEEV ', -INFO )
 
313
         RETURN
 
314
      ELSE IF( LQUERY ) THEN
 
315
         RETURN
 
316
      END IF
 
317
*
 
318
*     Quick return if possible
 
319
*
 
320
      IF( N.EQ.0 )
 
321
     $   RETURN
 
322
*
 
323
*     Get machine constants
 
324
*
 
325
      EPS = SLAMCH( 'P' )
 
326
      SMLNUM = SLAMCH( 'S' )
 
327
      BIGNUM = ONE / SMLNUM
 
328
      CALL SLABAD( SMLNUM, BIGNUM )
 
329
      SMLNUM = SQRT( SMLNUM ) / EPS
 
330
      BIGNUM = ONE / SMLNUM
 
331
*
 
332
*     Scale A if max element outside range [SMLNUM,BIGNUM]
 
333
*
 
334
      ANRM = SLANGE( 'M', N, N, A, LDA, DUM )
 
335
      SCALEA = .FALSE.
 
336
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 
337
         SCALEA = .TRUE.
 
338
         CSCALE = SMLNUM
 
339
      ELSE IF( ANRM.GT.BIGNUM ) THEN
 
340
         SCALEA = .TRUE.
 
341
         CSCALE = BIGNUM
 
342
      END IF
 
343
      IF( SCALEA )
 
344
     $   CALL SLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
 
345
*
 
346
*     Balance the matrix
 
347
*     (Workspace: need N)
 
348
*
 
349
      IBAL = 1
 
350
      CALL SGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
 
351
*
 
352
*     Reduce to upper Hessenberg form
 
353
*     (Workspace: need 3*N, prefer 2*N+N*NB)
 
354
*
 
355
      ITAU = IBAL + N
 
356
      IWRK = ITAU + N
 
357
      CALL SGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
 
358
     $             LWORK-IWRK+1, IERR )
 
359
*
 
360
      IF( WANTVL ) THEN
 
361
*
 
362
*        Want left eigenvectors
 
363
*        Copy Householder vectors to VL
 
364
*
 
365
         SIDE = 'L'
 
366
         CALL SLACPY( 'L', N, N, A, LDA, VL, LDVL )
 
367
*
 
368
*        Generate orthogonal matrix in VL
 
369
*        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
 
370
*
 
371
         CALL SORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
 
372
     $                LWORK-IWRK+1, IERR )
 
373
*
 
374
*        Perform QR iteration, accumulating Schur vectors in VL
 
375
*        (Workspace: need N+1, prefer N+HSWORK (see comments) )
 
376
*
 
377
         IWRK = ITAU
 
378
         CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
 
379
     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
 
380
*
 
381
         IF( WANTVR ) THEN
 
382
*
 
383
*           Want left and right eigenvectors
 
384
*           Copy Schur vectors to VR
 
385
*
 
386
            SIDE = 'B'
 
387
            CALL SLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
 
388
         END IF
 
389
*
 
390
      ELSE IF( WANTVR ) THEN
 
391
*
 
392
*        Want right eigenvectors
 
393
*        Copy Householder vectors to VR
 
394
*
 
395
         SIDE = 'R'
 
396
         CALL SLACPY( 'L', N, N, A, LDA, VR, LDVR )
 
397
*
 
398
*        Generate orthogonal matrix in VR
 
399
*        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
 
400
*
 
401
         CALL SORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
 
402
     $                LWORK-IWRK+1, IERR )
 
403
*
 
404
*        Perform QR iteration, accumulating Schur vectors in VR
 
405
*        (Workspace: need N+1, prefer N+HSWORK (see comments) )
 
406
*
 
407
         IWRK = ITAU
 
408
         CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
 
409
     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
 
410
*
 
411
      ELSE
 
412
*
 
413
*        Compute eigenvalues only
 
414
*        (Workspace: need N+1, prefer N+HSWORK (see comments) )
 
415
*
 
416
         IWRK = ITAU
 
417
         CALL SHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
 
418
     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
 
419
      END IF
 
420
*
 
421
*     If INFO > 0 from SHSEQR, then quit
 
422
*
 
423
      IF( INFO.GT.0 )
 
424
     $   GO TO 50
 
425
*
 
426
      IF( WANTVL .OR. WANTVR ) THEN
 
427
*
 
428
*        Compute left and/or right eigenvectors
 
429
*        (Workspace: need 4*N)
 
430
*
 
431
         CALL STREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
 
432
     $                N, NOUT, WORK( IWRK ), IERR )
 
433
      END IF
 
434
*
 
435
      IF( WANTVL ) THEN
 
436
*
 
437
*        Undo balancing of left eigenvectors
 
438
*        (Workspace: need N)
 
439
*
 
440
         CALL SGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
 
441
     $                IERR )
 
442
*
 
443
*        Normalize left eigenvectors and make largest component real
 
444
*
 
445
         DO 20 I = 1, N
 
446
            IF( WI( I ).EQ.ZERO ) THEN
 
447
               SCL = ONE / SNRM2( N, VL( 1, I ), 1 )
 
448
               CALL SSCAL( N, SCL, VL( 1, I ), 1 )
 
449
            ELSE IF( WI( I ).GT.ZERO ) THEN
 
450
               SCL = ONE / SLAPY2( SNRM2( N, VL( 1, I ), 1 ),
 
451
     $               SNRM2( N, VL( 1, I+1 ), 1 ) )
 
452
               CALL SSCAL( N, SCL, VL( 1, I ), 1 )
 
453
               CALL SSCAL( N, SCL, VL( 1, I+1 ), 1 )
 
454
               DO 10 K = 1, N
 
455
                  WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
 
456
   10          CONTINUE
 
457
               K = ISAMAX( N, WORK( IWRK ), 1 )
 
458
               CALL SLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
 
459
               CALL SROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
 
460
               VL( K, I+1 ) = ZERO
 
461
            END IF
 
462
   20    CONTINUE
 
463
      END IF
 
464
*
 
465
      IF( WANTVR ) THEN
 
466
*
 
467
*        Undo balancing of right eigenvectors
 
468
*        (Workspace: need N)
 
469
*
 
470
         CALL SGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
 
471
     $                IERR )
 
472
*
 
473
*        Normalize right eigenvectors and make largest component real
 
474
*
 
475
         DO 40 I = 1, N
 
476
            IF( WI( I ).EQ.ZERO ) THEN
 
477
               SCL = ONE / SNRM2( N, VR( 1, I ), 1 )
 
478
               CALL SSCAL( N, SCL, VR( 1, I ), 1 )
 
479
            ELSE IF( WI( I ).GT.ZERO ) THEN
 
480
               SCL = ONE / SLAPY2( SNRM2( N, VR( 1, I ), 1 ),
 
481
     $               SNRM2( N, VR( 1, I+1 ), 1 ) )
 
482
               CALL SSCAL( N, SCL, VR( 1, I ), 1 )
 
483
               CALL SSCAL( N, SCL, VR( 1, I+1 ), 1 )
 
484
               DO 30 K = 1, N
 
485
                  WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
 
486
   30          CONTINUE
 
487
               K = ISAMAX( N, WORK( IWRK ), 1 )
 
488
               CALL SLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
 
489
               CALL SROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
 
490
               VR( K, I+1 ) = ZERO
 
491
            END IF
 
492
   40    CONTINUE
 
493
      END IF
 
494
*
 
495
*     Undo scaling if necessary
 
496
*
 
497
   50 CONTINUE
 
498
      IF( SCALEA ) THEN
 
499
         CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
 
500
     $                MAX( N-INFO, 1 ), IERR )
 
501
         CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
 
502
     $                MAX( N-INFO, 1 ), IERR )
 
503
         IF( INFO.GT.0 ) THEN
 
504
            CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
 
505
     $                   IERR )
 
506
            CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
 
507
     $                   IERR )
 
508
         END IF
 
509
      END IF
 
510
*
 
511
      WORK( 1 ) = MAXWRK
 
512
      RETURN
 
513
*
 
514
*     End of SGEEV
 
515
*
 
516
      END