← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/linalg/dstedc.f

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/linalg/dstedc.f
 
1
*> \brief \b DSTEBZ
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*> \htmlonly
 
9
*> Download DSTEDC + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstedc.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstedc.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstedc.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
 
22
*                          LIWORK, INFO )
 
23
* 
 
24
*       .. Scalar Arguments ..
 
25
*       CHARACTER          COMPZ
 
26
*       INTEGER            INFO, LDZ, LIWORK, LWORK, N
 
27
*       ..
 
28
*       .. Array Arguments ..
 
29
*       INTEGER            IWORK( * )
 
30
*       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
 
31
*       ..
 
32
*  
 
33
*
 
34
*> \par Purpose:
 
35
*  =============
 
36
*>
 
37
*> \verbatim
 
38
*>
 
39
*> DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
 
40
*> symmetric tridiagonal matrix using the divide and conquer method.
 
41
*> The eigenvectors of a full or band real symmetric matrix can also be
 
42
*> found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
 
43
*> matrix to tridiagonal form.
 
44
*>
 
45
*> This code makes very mild assumptions about floating point
 
46
*> arithmetic. It will work on machines with a guard digit in
 
47
*> add/subtract, or on those binary machines without guard digits
 
48
*> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
 
49
*> It could conceivably fail on hexadecimal or decimal machines
 
50
*> without guard digits, but we know of none.  See DLAED3 for details.
 
51
*> \endverbatim
 
52
*
 
53
*  Arguments:
 
54
*  ==========
 
55
*
 
56
*> \param[in] COMPZ
 
57
*> \verbatim
 
58
*>          COMPZ is CHARACTER*1
 
59
*>          = 'N':  Compute eigenvalues only.
 
60
*>          = 'I':  Compute eigenvectors of tridiagonal matrix also.
 
61
*>          = 'V':  Compute eigenvectors of original dense symmetric
 
62
*>                  matrix also.  On entry, Z contains the orthogonal
 
63
*>                  matrix used to reduce the original matrix to
 
64
*>                  tridiagonal form.
 
65
*> \endverbatim
 
66
*>
 
67
*> \param[in] N
 
68
*> \verbatim
 
69
*>          N is INTEGER
 
70
*>          The dimension of the symmetric tridiagonal matrix.  N >= 0.
 
71
*> \endverbatim
 
72
*>
 
73
*> \param[in,out] D
 
74
*> \verbatim
 
75
*>          D is DOUBLE PRECISION array, dimension (N)
 
76
*>          On entry, the diagonal elements of the tridiagonal matrix.
 
77
*>          On exit, if INFO = 0, the eigenvalues in ascending order.
 
78
*> \endverbatim
 
79
*>
 
80
*> \param[in,out] E
 
81
*> \verbatim
 
82
*>          E is DOUBLE PRECISION array, dimension (N-1)
 
83
*>          On entry, the subdiagonal elements of the tridiagonal matrix.
 
84
*>          On exit, E has been destroyed.
 
85
*> \endverbatim
 
86
*>
 
87
*> \param[in,out] Z
 
88
*> \verbatim
 
89
*>          Z is DOUBLE PRECISION array, dimension (LDZ,N)
 
90
*>          On entry, if COMPZ = 'V', then Z contains the orthogonal
 
91
*>          matrix used in the reduction to tridiagonal form.
 
92
*>          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
 
93
*>          orthonormal eigenvectors of the original symmetric matrix,
 
94
*>          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
 
95
*>          of the symmetric tridiagonal matrix.
 
96
*>          If  COMPZ = 'N', then Z is not referenced.
 
97
*> \endverbatim
 
98
*>
 
99
*> \param[in] LDZ
 
100
*> \verbatim
 
101
*>          LDZ is INTEGER
 
102
*>          The leading dimension of the array Z.  LDZ >= 1.
 
103
*>          If eigenvectors are desired, then LDZ >= max(1,N).
 
104
*> \endverbatim
 
105
*>
 
106
*> \param[out] WORK
 
107
*> \verbatim
 
108
*>          WORK is DOUBLE PRECISION array,
 
109
*>                                         dimension (LWORK)
 
110
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 
111
*> \endverbatim
 
112
*>
 
113
*> \param[in] LWORK
 
114
*> \verbatim
 
115
*>          LWORK is INTEGER
 
116
*>          The dimension of the array WORK.
 
117
*>          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
 
118
*>          If COMPZ = 'V' and N > 1 then LWORK must be at least
 
119
*>                         ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
 
120
*>                         where lg( N ) = smallest integer k such
 
121
*>                         that 2**k >= N.
 
122
*>          If COMPZ = 'I' and N > 1 then LWORK must be at least
 
123
*>                         ( 1 + 4*N + N**2 ).
 
124
*>          Note that for COMPZ = 'I' or 'V', then if N is less than or
 
125
*>          equal to the minimum divide size, usually 25, then LWORK need
 
126
*>          only be max(1,2*(N-1)).
 
127
*>
 
128
*>          If LWORK = -1, then a workspace query is assumed; the routine
 
129
*>          only calculates the optimal size of the WORK array, returns
 
130
*>          this value as the first entry of the WORK array, and no error
 
131
*>          message related to LWORK is issued by XERBLA.
 
132
*> \endverbatim
 
133
*>
 
134
*> \param[out] IWORK
 
135
*> \verbatim
 
136
*>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
 
137
*>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 
138
*> \endverbatim
 
139
*>
 
140
*> \param[in] LIWORK
 
141
*> \verbatim
 
142
*>          LIWORK is INTEGER
 
143
*>          The dimension of the array IWORK.
 
144
*>          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
 
145
*>          If COMPZ = 'V' and N > 1 then LIWORK must be at least
 
146
*>                         ( 6 + 6*N + 5*N*lg N ).
 
147
*>          If COMPZ = 'I' and N > 1 then LIWORK must be at least
 
148
*>                         ( 3 + 5*N ).
 
149
*>          Note that for COMPZ = 'I' or 'V', then if N is less than or
 
150
*>          equal to the minimum divide size, usually 25, then LIWORK
 
151
*>          need only be 1.
 
152
*>
 
153
*>          If LIWORK = -1, then a workspace query is assumed; the
 
154
*>          routine only calculates the optimal size of the IWORK array,
 
155
*>          returns this value as the first entry of the IWORK array, and
 
156
*>          no error message related to LIWORK is issued by XERBLA.
 
157
*> \endverbatim
 
158
*>
 
159
*> \param[out] INFO
 
160
*> \verbatim
 
161
*>          INFO is INTEGER
 
162
*>          = 0:  successful exit.
 
163
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 
164
*>          > 0:  The algorithm failed to compute an eigenvalue while
 
165
*>                working on the submatrix lying in rows and columns
 
166
*>                INFO/(N+1) through mod(INFO,N+1).
 
167
*> \endverbatim
 
168
*
 
169
*  Authors:
 
170
*  ========
 
171
*
 
172
*> \author Univ. of Tennessee 
 
173
*> \author Univ. of California Berkeley 
 
174
*> \author Univ. of Colorado Denver 
 
175
*> \author NAG Ltd. 
 
176
*
 
177
*> \date November 2011
 
178
*
 
179
*> \ingroup auxOTHERcomputational
 
180
*
 
181
*> \par Contributors:
 
182
*  ==================
 
183
*>
 
184
*> Jeff Rutter, Computer Science Division, University of California
 
185
*> at Berkeley, USA \n
 
186
*>  Modified by Francoise Tisseur, University of Tennessee
 
187
*>
 
188
*  =====================================================================
 
189
      SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
 
190
     $                   LIWORK, INFO )
 
191
*
 
192
*  -- LAPACK computational routine (version 3.4.0) --
 
193
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
194
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
195
*     November 2011
 
196
*
 
197
*     .. Scalar Arguments ..
 
198
      CHARACTER          COMPZ
 
199
      INTEGER            INFO, LDZ, LIWORK, LWORK, N
 
200
*     ..
 
201
*     .. Array Arguments ..
 
202
      INTEGER            IWORK( * )
 
203
      DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
 
204
*     ..
 
205
*
 
206
*  =====================================================================
 
207
*
 
208
*     .. Parameters ..
 
209
      DOUBLE PRECISION   ZERO, ONE, TWO
 
210
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
 
211
*     ..
 
212
*     .. Local Scalars ..
 
213
      LOGICAL            LQUERY
 
214
      INTEGER            FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
 
215
     $                   LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
 
216
      DOUBLE PRECISION   EPS, ORGNRM, P, TINY
 
217
*     ..
 
218
*     .. External Functions ..
 
219
      LOGICAL            LSAME
 
220
      INTEGER            ILAENV
 
221
      DOUBLE PRECISION   DLAMCH, DLANST
 
222
      EXTERNAL           LSAME, ILAENV, DLAMCH, DLANST
 
223
*     ..
 
224
*     .. External Subroutines ..
 
225
      EXTERNAL           DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
 
226
     $                   DSTEQR, DSTERF, DSWAP, XERBLA
 
227
*     ..
 
228
*     .. Intrinsic Functions ..
 
229
      INTRINSIC          ABS, DBLE, INT, LOG, MAX, MOD, SQRT
 
230
*     ..
 
231
*     .. Executable Statements ..
 
232
*
 
233
*     Test the input parameters.
 
234
*
 
235
      INFO = 0
 
236
      LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
 
237
*
 
238
      IF( LSAME( COMPZ, 'N' ) ) THEN
 
239
         ICOMPZ = 0
 
240
      ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
 
241
         ICOMPZ = 1
 
242
      ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
 
243
         ICOMPZ = 2
 
244
      ELSE
 
245
         ICOMPZ = -1
 
246
      END IF
 
247
      IF( ICOMPZ.LT.0 ) THEN
 
248
         INFO = -1
 
249
      ELSE IF( N.LT.0 ) THEN
 
250
         INFO = -2
 
251
      ELSE IF( ( LDZ.LT.1 ) .OR.
 
252
     $         ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
 
253
         INFO = -6
 
254
      END IF
 
255
*
 
256
      IF( INFO.EQ.0 ) THEN
 
257
*
 
258
*        Compute the workspace requirements
 
259
*
 
260
         SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
 
261
         IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
 
262
            LIWMIN = 1
 
263
            LWMIN = 1
 
264
         ELSE IF( N.LE.SMLSIZ ) THEN
 
265
            LIWMIN = 1
 
266
            LWMIN = 2*( N - 1 )
 
267
         ELSE
 
268
            LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
 
269
            IF( 2**LGN.LT.N )
 
270
     $         LGN = LGN + 1
 
271
            IF( 2**LGN.LT.N )
 
272
     $         LGN = LGN + 1
 
273
            IF( ICOMPZ.EQ.1 ) THEN
 
274
               LWMIN = 1 + 3*N + 2*N*LGN + 4*N**2
 
275
               LIWMIN = 6 + 6*N + 5*N*LGN
 
276
            ELSE IF( ICOMPZ.EQ.2 ) THEN
 
277
               LWMIN = 1 + 4*N + N**2
 
278
               LIWMIN = 3 + 5*N
 
279
            END IF
 
280
         END IF
 
281
         WORK( 1 ) = LWMIN
 
282
         IWORK( 1 ) = LIWMIN
 
283
*
 
284
         IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
 
285
            INFO = -8
 
286
         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
 
287
            INFO = -10
 
288
         END IF
 
289
      END IF
 
290
*
 
291
      IF( INFO.NE.0 ) THEN
 
292
         CALL XERBLA( 'DSTEDC', -INFO )
 
293
         RETURN
 
294
      ELSE IF (LQUERY) THEN
 
295
         RETURN
 
296
      END IF
 
297
*
 
298
*     Quick return if possible
 
299
*
 
300
      IF( N.EQ.0 )
 
301
     $   RETURN
 
302
      IF( N.EQ.1 ) THEN
 
303
         IF( ICOMPZ.NE.0 )
 
304
     $      Z( 1, 1 ) = ONE
 
305
         RETURN
 
306
      END IF
 
307
*
 
308
*     If the following conditional clause is removed, then the routine
 
309
*     will use the Divide and Conquer routine to compute only the
 
310
*     eigenvalues, which requires (3N + 3N**2) real workspace and
 
311
*     (2 + 5N + 2N lg(N)) integer workspace.
 
312
*     Since on many architectures DSTERF is much faster than any other
 
313
*     algorithm for finding eigenvalues only, it is used here
 
314
*     as the default. If the conditional clause is removed, then
 
315
*     information on the size of workspace needs to be changed.
 
316
*
 
317
*     If COMPZ = 'N', use DSTERF to compute the eigenvalues.
 
318
*
 
319
      IF( ICOMPZ.EQ.0 ) THEN
 
320
         CALL DSTERF( N, D, E, INFO )
 
321
         GO TO 50
 
322
      END IF
 
323
*
 
324
*     If N is smaller than the minimum divide size (SMLSIZ+1), then
 
325
*     solve the problem with another solver.
 
326
*
 
327
      IF( N.LE.SMLSIZ ) THEN
 
328
*
 
329
         CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
 
330
*
 
331
      ELSE
 
332
*
 
333
*        If COMPZ = 'V', the Z matrix must be stored elsewhere for later
 
334
*        use.
 
335
*
 
336
         IF( ICOMPZ.EQ.1 ) THEN
 
337
            STOREZ = 1 + N*N
 
338
         ELSE
 
339
            STOREZ = 1
 
340
         END IF
 
341
*
 
342
         IF( ICOMPZ.EQ.2 ) THEN
 
343
            CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
 
344
         END IF
 
345
*
 
346
*        Scale.
 
347
*
 
348
         ORGNRM = DLANST( 'M', N, D, E )
 
349
         IF( ORGNRM.EQ.ZERO )
 
350
     $      GO TO 50
 
351
*
 
352
         EPS = DLAMCH( 'Epsilon' )
 
353
*
 
354
         START = 1
 
355
*
 
356
*        while ( START <= N )
 
357
*
 
358
   10    CONTINUE
 
359
         IF( START.LE.N ) THEN
 
360
*
 
361
*           Let FINISH be the position of the next subdiagonal entry
 
362
*           such that E( FINISH ) <= TINY or FINISH = N if no such
 
363
*           subdiagonal exists.  The matrix identified by the elements
 
364
*           between START and FINISH constitutes an independent
 
365
*           sub-problem.
 
366
*
 
367
            FINISH = START
 
368
   20       CONTINUE
 
369
            IF( FINISH.LT.N ) THEN
 
370
               TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
 
371
     $                    SQRT( ABS( D( FINISH+1 ) ) )
 
372
               IF( ABS( E( FINISH ) ).GT.TINY ) THEN
 
373
                  FINISH = FINISH + 1
 
374
                  GO TO 20
 
375
               END IF
 
376
            END IF
 
377
*
 
378
*           (Sub) Problem determined.  Compute its size and solve it.
 
379
*
 
380
            M = FINISH - START + 1
 
381
            IF( M.EQ.1 ) THEN
 
382
               START = FINISH + 1
 
383
               GO TO 10
 
384
            END IF
 
385
            IF( M.GT.SMLSIZ ) THEN
 
386
*
 
387
*              Scale.
 
388
*
 
389
               ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
 
390
               CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
 
391
     $                      INFO )
 
392
               CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
 
393
     $                      M-1, INFO )
 
394
*
 
395
               IF( ICOMPZ.EQ.1 ) THEN
 
396
                  STRTRW = 1
 
397
               ELSE
 
398
                  STRTRW = START
 
399
               END IF
 
400
               CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
 
401
     $                      Z( STRTRW, START ), LDZ, WORK( 1 ), N,
 
402
     $                      WORK( STOREZ ), IWORK, INFO )
 
403
               IF( INFO.NE.0 ) THEN
 
404
                  INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
 
405
     $                   MOD( INFO, ( M+1 ) ) + START - 1
 
406
                  GO TO 50
 
407
               END IF
 
408
*
 
409
*              Scale back.
 
410
*
 
411
               CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
 
412
     $                      INFO )
 
413
*
 
414
            ELSE
 
415
               IF( ICOMPZ.EQ.1 ) THEN
 
416
*
 
417
*                 Since QR won't update a Z matrix which is larger than
 
418
*                 the length of D, we must solve the sub-problem in a
 
419
*                 workspace and then multiply back into Z.
 
420
*
 
421
                  CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
 
422
     $                         WORK( M*M+1 ), INFO )
 
423
                  CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
 
424
     $                         WORK( STOREZ ), N )
 
425
                  CALL DGEMM( 'N', 'N', N, M, M, ONE,
 
426
     $                        WORK( STOREZ ), N, WORK, M, ZERO,
 
427
     $                        Z( 1, START ), LDZ )
 
428
               ELSE IF( ICOMPZ.EQ.2 ) THEN
 
429
                  CALL DSTEQR( 'I', M, D( START ), E( START ),
 
430
     $                         Z( START, START ), LDZ, WORK, INFO )
 
431
               ELSE
 
432
                  CALL DSTERF( M, D( START ), E( START ), INFO )
 
433
               END IF
 
434
               IF( INFO.NE.0 ) THEN
 
435
                  INFO = START*( N+1 ) + FINISH
 
436
                  GO TO 50
 
437
               END IF
 
438
            END IF
 
439
*
 
440
            START = FINISH + 1
 
441
            GO TO 10
 
442
         END IF
 
443
*
 
444
*        endwhile
 
445
*
 
446
*        If the problem split any number of times, then the eigenvalues
 
447
*        will not be properly ordered.  Here we permute the eigenvalues
 
448
*        (and the associated eigenvectors) into ascending order.
 
449
*
 
450
         IF( M.NE.N ) THEN
 
451
            IF( ICOMPZ.EQ.0 ) THEN
 
452
*
 
453
*              Use Quick Sort
 
454
*
 
455
               CALL DLASRT( 'I', N, D, INFO )
 
456
*
 
457
            ELSE
 
458
*
 
459
*              Use Selection Sort to minimize swaps of eigenvectors
 
460
*
 
461
               DO 40 II = 2, N
 
462
                  I = II - 1
 
463
                  K = I
 
464
                  P = D( I )
 
465
                  DO 30 J = II, N
 
466
                     IF( D( J ).LT.P ) THEN
 
467
                        K = J
 
468
                        P = D( J )
 
469
                     END IF
 
470
   30             CONTINUE
 
471
                  IF( K.NE.I ) THEN
 
472
                     D( K ) = D( I )
 
473
                     D( I ) = P
 
474
                     CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
 
475
                  END IF
 
476
   40          CONTINUE
 
477
            END IF
 
478
         END IF
 
479
      END IF
 
480
*
 
481
   50 CONTINUE
 
482
      WORK( 1 ) = LWMIN
 
483
      IWORK( 1 ) = LIWMIN
 
484
*
 
485
      RETURN
 
486
*
 
487
*     End of DSTEDC
 
488
*
 
489
      END