← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/linalg/dlaed0.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/dlaed0.f
 
1
*> \brief \b DLAED0 used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method.
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*> \htmlonly
 
9
*> Download DLAED0 + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed0.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed0.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed0.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
 
22
*                          WORK, IWORK, INFO )
 
23
* 
 
24
*       .. Scalar Arguments ..
 
25
*       INTEGER            ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
 
26
*       ..
 
27
*       .. Array Arguments ..
 
28
*       INTEGER            IWORK( * )
 
29
*       DOUBLE PRECISION   D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
 
30
*      $                   WORK( * )
 
31
*       ..
 
32
*  
 
33
*
 
34
*> \par Purpose:
 
35
*  =============
 
36
*>
 
37
*> \verbatim
 
38
*>
 
39
*> DLAED0 computes all eigenvalues and corresponding eigenvectors of a
 
40
*> symmetric tridiagonal matrix using the divide and conquer method.
 
41
*> \endverbatim
 
42
*
 
43
*  Arguments:
 
44
*  ==========
 
45
*
 
46
*> \param[in] ICOMPQ
 
47
*> \verbatim
 
48
*>          ICOMPQ is INTEGER
 
49
*>          = 0:  Compute eigenvalues only.
 
50
*>          = 1:  Compute eigenvectors of original dense symmetric matrix
 
51
*>                also.  On entry, Q contains the orthogonal matrix used
 
52
*>                to reduce the original matrix to tridiagonal form.
 
53
*>          = 2:  Compute eigenvalues and eigenvectors of tridiagonal
 
54
*>                matrix.
 
55
*> \endverbatim
 
56
*>
 
57
*> \param[in] QSIZ
 
58
*> \verbatim
 
59
*>          QSIZ is INTEGER
 
60
*>         The dimension of the orthogonal matrix used to reduce
 
61
*>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
 
62
*> \endverbatim
 
63
*>
 
64
*> \param[in] N
 
65
*> \verbatim
 
66
*>          N is INTEGER
 
67
*>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
 
68
*> \endverbatim
 
69
*>
 
70
*> \param[in,out] D
 
71
*> \verbatim
 
72
*>          D is DOUBLE PRECISION array, dimension (N)
 
73
*>         On entry, the main diagonal of the tridiagonal matrix.
 
74
*>         On exit, its eigenvalues.
 
75
*> \endverbatim
 
76
*>
 
77
*> \param[in] E
 
78
*> \verbatim
 
79
*>          E is DOUBLE PRECISION array, dimension (N-1)
 
80
*>         The off-diagonal elements of the tridiagonal matrix.
 
81
*>         On exit, E has been destroyed.
 
82
*> \endverbatim
 
83
*>
 
84
*> \param[in,out] Q
 
85
*> \verbatim
 
86
*>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
 
87
*>         On entry, Q must contain an N-by-N orthogonal matrix.
 
88
*>         If ICOMPQ = 0    Q is not referenced.
 
89
*>         If ICOMPQ = 1    On entry, Q is a subset of the columns of the
 
90
*>                          orthogonal matrix used to reduce the full
 
91
*>                          matrix to tridiagonal form corresponding to
 
92
*>                          the subset of the full matrix which is being
 
93
*>                          decomposed at this time.
 
94
*>         If ICOMPQ = 2    On entry, Q will be the identity matrix.
 
95
*>                          On exit, Q contains the eigenvectors of the
 
96
*>                          tridiagonal matrix.
 
97
*> \endverbatim
 
98
*>
 
99
*> \param[in] LDQ
 
100
*> \verbatim
 
101
*>          LDQ is INTEGER
 
102
*>         The leading dimension of the array Q.  If eigenvectors are
 
103
*>         desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.
 
104
*> \endverbatim
 
105
*>
 
106
*> \param[out] QSTORE
 
107
*> \verbatim
 
108
*>          QSTORE is DOUBLE PRECISION array, dimension (LDQS, N)
 
109
*>         Referenced only when ICOMPQ = 1.  Used to store parts of
 
110
*>         the eigenvector matrix when the updating matrix multiplies
 
111
*>         take place.
 
112
*> \endverbatim
 
113
*>
 
114
*> \param[in] LDQS
 
115
*> \verbatim
 
116
*>          LDQS is INTEGER
 
117
*>         The leading dimension of the array QSTORE.  If ICOMPQ = 1,
 
118
*>         then  LDQS >= max(1,N).  In any case,  LDQS >= 1.
 
119
*> \endverbatim
 
120
*>
 
121
*> \param[out] WORK
 
122
*> \verbatim
 
123
*>          WORK is DOUBLE PRECISION array,
 
124
*>         If ICOMPQ = 0 or 1, the dimension of WORK must be at least
 
125
*>                     1 + 3*N + 2*N*lg N + 3*N**2
 
126
*>                     ( lg( N ) = smallest integer k
 
127
*>                                 such that 2^k >= N )
 
128
*>         If ICOMPQ = 2, the dimension of WORK must be at least
 
129
*>                     4*N + N**2.
 
130
*> \endverbatim
 
131
*>
 
132
*> \param[out] IWORK
 
133
*> \verbatim
 
134
*>          IWORK is INTEGER array,
 
135
*>         If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
 
136
*>                        6 + 6*N + 5*N*lg N.
 
137
*>                        ( lg( N ) = smallest integer k
 
138
*>                                    such that 2^k >= N )
 
139
*>         If ICOMPQ = 2, the dimension of IWORK must be at least
 
140
*>                        3 + 5*N.
 
141
*> \endverbatim
 
142
*>
 
143
*> \param[out] INFO
 
144
*> \verbatim
 
145
*>          INFO is INTEGER
 
146
*>          = 0:  successful exit.
 
147
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 
148
*>          > 0:  The algorithm failed to compute an eigenvalue while
 
149
*>                working on the submatrix lying in rows and columns
 
150
*>                INFO/(N+1) through mod(INFO,N+1).
 
151
*> \endverbatim
 
152
*
 
153
*  Authors:
 
154
*  ========
 
155
*
 
156
*> \author Univ. of Tennessee 
 
157
*> \author Univ. of California Berkeley 
 
158
*> \author Univ. of Colorado Denver 
 
159
*> \author NAG Ltd. 
 
160
*
 
161
*> \date September 2012
 
162
*
 
163
*> \ingroup auxOTHERcomputational
 
164
*
 
165
*> \par Contributors:
 
166
*  ==================
 
167
*>
 
168
*> Jeff Rutter, Computer Science Division, University of California
 
169
*> at Berkeley, USA
 
170
*
 
171
*  =====================================================================
 
172
      SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
 
173
     $                   WORK, IWORK, INFO )
 
174
*
 
175
*  -- LAPACK computational routine (version 3.4.2) --
 
176
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
177
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
178
*     September 2012
 
179
*
 
180
*     .. Scalar Arguments ..
 
181
      INTEGER            ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
 
182
*     ..
 
183
*     .. Array Arguments ..
 
184
      INTEGER            IWORK( * )
 
185
      DOUBLE PRECISION   D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
 
186
     $                   WORK( * )
 
187
*     ..
 
188
*
 
189
*  =====================================================================
 
190
*
 
191
*     .. Parameters ..
 
192
      DOUBLE PRECISION   ZERO, ONE, TWO
 
193
      PARAMETER          ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 )
 
194
*     ..
 
195
*     .. Local Scalars ..
 
196
      INTEGER            CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
 
197
     $                   IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
 
198
     $                   J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
 
199
     $                   SPM2, SUBMAT, SUBPBS, TLVLS
 
200
      DOUBLE PRECISION   TEMP
 
201
*     ..
 
202
*     .. External Subroutines ..
 
203
      EXTERNAL           DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR,
 
204
     $                   XERBLA
 
205
*     ..
 
206
*     .. External Functions ..
 
207
      INTEGER            ILAENV
 
208
      EXTERNAL           ILAENV
 
209
*     ..
 
210
*     .. Intrinsic Functions ..
 
211
      INTRINSIC          ABS, DBLE, INT, LOG, MAX
 
212
*     ..
 
213
*     .. Executable Statements ..
 
214
*
 
215
*     Test the input parameters.
 
216
*
 
217
      INFO = 0
 
218
*
 
219
      IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
 
220
         INFO = -1
 
221
      ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
 
222
         INFO = -2
 
223
      ELSE IF( N.LT.0 ) THEN
 
224
         INFO = -3
 
225
      ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
 
226
         INFO = -7
 
227
      ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
 
228
         INFO = -9
 
229
      END IF
 
230
      IF( INFO.NE.0 ) THEN
 
231
         CALL XERBLA( 'DLAED0', -INFO )
 
232
         RETURN
 
233
      END IF
 
234
*
 
235
*     Quick return if possible
 
236
*
 
237
      IF( N.EQ.0 )
 
238
     $   RETURN
 
239
*
 
240
      SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
 
241
*
 
242
*     Determine the size and placement of the submatrices, and save in
 
243
*     the leading elements of IWORK.
 
244
*
 
245
      IWORK( 1 ) = N
 
246
      SUBPBS = 1
 
247
      TLVLS = 0
 
248
   10 CONTINUE
 
249
      IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
 
250
         DO 20 J = SUBPBS, 1, -1
 
251
            IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
 
252
            IWORK( 2*J-1 ) = IWORK( J ) / 2
 
253
   20    CONTINUE
 
254
         TLVLS = TLVLS + 1
 
255
         SUBPBS = 2*SUBPBS
 
256
         GO TO 10
 
257
      END IF
 
258
      DO 30 J = 2, SUBPBS
 
259
         IWORK( J ) = IWORK( J ) + IWORK( J-1 )
 
260
   30 CONTINUE
 
261
*
 
262
*     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
 
263
*     using rank-1 modifications (cuts).
 
264
*
 
265
      SPM1 = SUBPBS - 1
 
266
      DO 40 I = 1, SPM1
 
267
         SUBMAT = IWORK( I ) + 1
 
268
         SMM1 = SUBMAT - 1
 
269
         D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
 
270
         D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
 
271
   40 CONTINUE
 
272
*
 
273
      INDXQ = 4*N + 3
 
274
      IF( ICOMPQ.NE.2 ) THEN
 
275
*
 
276
*        Set up workspaces for eigenvalues only/accumulate new vectors
 
277
*        routine
 
278
*
 
279
         TEMP = LOG( DBLE( N ) ) / LOG( TWO )
 
280
         LGN = INT( TEMP )
 
281
         IF( 2**LGN.LT.N )
 
282
     $      LGN = LGN + 1
 
283
         IF( 2**LGN.LT.N )
 
284
     $      LGN = LGN + 1
 
285
         IPRMPT = INDXQ + N + 1
 
286
         IPERM = IPRMPT + N*LGN
 
287
         IQPTR = IPERM + N*LGN
 
288
         IGIVPT = IQPTR + N + 2
 
289
         IGIVCL = IGIVPT + N*LGN
 
290
*
 
291
         IGIVNM = 1
 
292
         IQ = IGIVNM + 2*N*LGN
 
293
         IWREM = IQ + N**2 + 1
 
294
*
 
295
*        Initialize pointers
 
296
*
 
297
         DO 50 I = 0, SUBPBS
 
298
            IWORK( IPRMPT+I ) = 1
 
299
            IWORK( IGIVPT+I ) = 1
 
300
   50    CONTINUE
 
301
         IWORK( IQPTR ) = 1
 
302
      END IF
 
303
*
 
304
*     Solve each submatrix eigenproblem at the bottom of the divide and
 
305
*     conquer tree.
 
306
*
 
307
      CURR = 0
 
308
      DO 70 I = 0, SPM1
 
309
         IF( I.EQ.0 ) THEN
 
310
            SUBMAT = 1
 
311
            MATSIZ = IWORK( 1 )
 
312
         ELSE
 
313
            SUBMAT = IWORK( I ) + 1
 
314
            MATSIZ = IWORK( I+1 ) - IWORK( I )
 
315
         END IF
 
316
         IF( ICOMPQ.EQ.2 ) THEN
 
317
            CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
 
318
     $                   Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
 
319
            IF( INFO.NE.0 )
 
320
     $         GO TO 130
 
321
         ELSE
 
322
            CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
 
323
     $                   WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
 
324
     $                   INFO )
 
325
            IF( INFO.NE.0 )
 
326
     $         GO TO 130
 
327
            IF( ICOMPQ.EQ.1 ) THEN
 
328
               CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
 
329
     $                     Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
 
330
     $                     CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
 
331
     $                     LDQS )
 
332
            END IF
 
333
            IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
 
334
            CURR = CURR + 1
 
335
         END IF
 
336
         K = 1
 
337
         DO 60 J = SUBMAT, IWORK( I+1 )
 
338
            IWORK( INDXQ+J ) = K
 
339
            K = K + 1
 
340
   60    CONTINUE
 
341
   70 CONTINUE
 
342
*
 
343
*     Successively merge eigensystems of adjacent submatrices
 
344
*     into eigensystem for the corresponding larger matrix.
 
345
*
 
346
*     while ( SUBPBS > 1 )
 
347
*
 
348
      CURLVL = 1
 
349
   80 CONTINUE
 
350
      IF( SUBPBS.GT.1 ) THEN
 
351
         SPM2 = SUBPBS - 2
 
352
         DO 90 I = 0, SPM2, 2
 
353
            IF( I.EQ.0 ) THEN
 
354
               SUBMAT = 1
 
355
               MATSIZ = IWORK( 2 )
 
356
               MSD2 = IWORK( 1 )
 
357
               CURPRB = 0
 
358
            ELSE
 
359
               SUBMAT = IWORK( I ) + 1
 
360
               MATSIZ = IWORK( I+2 ) - IWORK( I )
 
361
               MSD2 = MATSIZ / 2
 
362
               CURPRB = CURPRB + 1
 
363
            END IF
 
364
*
 
365
*     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
 
366
*     into an eigensystem of size MATSIZ.
 
367
*     DLAED1 is used only for the full eigensystem of a tridiagonal
 
368
*     matrix.
 
369
*     DLAED7 handles the cases in which eigenvalues only or eigenvalues
 
370
*     and eigenvectors of a full symmetric matrix (which was reduced to
 
371
*     tridiagonal form) are desired.
 
372
*
 
373
            IF( ICOMPQ.EQ.2 ) THEN
 
374
               CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
 
375
     $                      LDQ, IWORK( INDXQ+SUBMAT ),
 
376
     $                      E( SUBMAT+MSD2-1 ), MSD2, WORK,
 
377
     $                      IWORK( SUBPBS+1 ), INFO )
 
378
            ELSE
 
379
               CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
 
380
     $                      D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
 
381
     $                      IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
 
382
     $                      MSD2, WORK( IQ ), IWORK( IQPTR ),
 
383
     $                      IWORK( IPRMPT ), IWORK( IPERM ),
 
384
     $                      IWORK( IGIVPT ), IWORK( IGIVCL ),
 
385
     $                      WORK( IGIVNM ), WORK( IWREM ),
 
386
     $                      IWORK( SUBPBS+1 ), INFO )
 
387
            END IF
 
388
            IF( INFO.NE.0 )
 
389
     $         GO TO 130
 
390
            IWORK( I / 2+1 ) = IWORK( I+2 )
 
391
   90    CONTINUE
 
392
         SUBPBS = SUBPBS / 2
 
393
         CURLVL = CURLVL + 1
 
394
         GO TO 80
 
395
      END IF
 
396
*
 
397
*     end while
 
398
*
 
399
*     Re-merge the eigenvalues/vectors which were deflated at the final
 
400
*     merge step.
 
401
*
 
402
      IF( ICOMPQ.EQ.1 ) THEN
 
403
         DO 100 I = 1, N
 
404
            J = IWORK( INDXQ+I )
 
405
            WORK( I ) = D( J )
 
406
            CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
 
407
  100    CONTINUE
 
408
         CALL DCOPY( N, WORK, 1, D, 1 )
 
409
      ELSE IF( ICOMPQ.EQ.2 ) THEN
 
410
         DO 110 I = 1, N
 
411
            J = IWORK( INDXQ+I )
 
412
            WORK( I ) = D( J )
 
413
            CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
 
414
  110    CONTINUE
 
415
         CALL DCOPY( N, WORK, 1, D, 1 )
 
416
         CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
 
417
      ELSE
 
418
         DO 120 I = 1, N
 
419
            J = IWORK( INDXQ+I )
 
420
            WORK( I ) = D( J )
 
421
  120    CONTINUE
 
422
         CALL DCOPY( N, WORK, 1, D, 1 )
 
423
      END IF
 
424
      GO TO 140
 
425
*
 
426
  130 CONTINUE
 
427
      INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
 
428
*
 
429
  140 CONTINUE
 
430
      RETURN
 
431
*
 
432
*     End of DLAED0
 
433
*
 
434
      END