← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/linalg/dlaeda.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/dlaeda.f
 
1
*> \brief \b DLAEDA used by sstedc. Computes the Z vector determining the rank-one modification of the diagonal matrix. Used when the original matrix is dense.
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*> \htmlonly
 
9
*> Download DLAEDA + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaeda.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaeda.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaeda.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
 
22
*                          GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
 
23
* 
 
24
*       .. Scalar Arguments ..
 
25
*       INTEGER            CURLVL, CURPBM, INFO, N, TLVLS
 
26
*       ..
 
27
*       .. Array Arguments ..
 
28
*       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
 
29
*      $                   PRMPTR( * ), QPTR( * )
 
30
*       DOUBLE PRECISION   GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
 
31
*       ..
 
32
*  
 
33
*
 
34
*> \par Purpose:
 
35
*  =============
 
36
*>
 
37
*> \verbatim
 
38
*>
 
39
*> DLAEDA computes the Z vector corresponding to the merge step in the
 
40
*> CURLVLth step of the merge process with TLVLS steps for the CURPBMth
 
41
*> problem.
 
42
*> \endverbatim
 
43
*
 
44
*  Arguments:
 
45
*  ==========
 
46
*
 
47
*> \param[in] N
 
48
*> \verbatim
 
49
*>          N is INTEGER
 
50
*>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
 
51
*> \endverbatim
 
52
*>
 
53
*> \param[in] TLVLS
 
54
*> \verbatim
 
55
*>          TLVLS is INTEGER
 
56
*>         The total number of merging levels in the overall divide and
 
57
*>         conquer tree.
 
58
*> \endverbatim
 
59
*>
 
60
*> \param[in] CURLVL
 
61
*> \verbatim
 
62
*>          CURLVL is INTEGER
 
63
*>         The current level in the overall merge routine,
 
64
*>         0 <= curlvl <= tlvls.
 
65
*> \endverbatim
 
66
*>
 
67
*> \param[in] CURPBM
 
68
*> \verbatim
 
69
*>          CURPBM is INTEGER
 
70
*>         The current problem in the current level in the overall
 
71
*>         merge routine (counting from upper left to lower right).
 
72
*> \endverbatim
 
73
*>
 
74
*> \param[in] PRMPTR
 
75
*> \verbatim
 
76
*>          PRMPTR is INTEGER array, dimension (N lg N)
 
77
*>         Contains a list of pointers which indicate where in PERM a
 
78
*>         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
 
79
*>         indicates the size of the permutation and incidentally the
 
80
*>         size of the full, non-deflated problem.
 
81
*> \endverbatim
 
82
*>
 
83
*> \param[in] PERM
 
84
*> \verbatim
 
85
*>          PERM is INTEGER array, dimension (N lg N)
 
86
*>         Contains the permutations (from deflation and sorting) to be
 
87
*>         applied to each eigenblock.
 
88
*> \endverbatim
 
89
*>
 
90
*> \param[in] GIVPTR
 
91
*> \verbatim
 
92
*>          GIVPTR is INTEGER array, dimension (N lg N)
 
93
*>         Contains a list of pointers which indicate where in GIVCOL a
 
94
*>         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
 
95
*>         indicates the number of Givens rotations.
 
96
*> \endverbatim
 
97
*>
 
98
*> \param[in] GIVCOL
 
99
*> \verbatim
 
100
*>          GIVCOL is INTEGER array, dimension (2, N lg N)
 
101
*>         Each pair of numbers indicates a pair of columns to take place
 
102
*>         in a Givens rotation.
 
103
*> \endverbatim
 
104
*>
 
105
*> \param[in] GIVNUM
 
106
*> \verbatim
 
107
*>          GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)
 
108
*>         Each number indicates the S value to be used in the
 
109
*>         corresponding Givens rotation.
 
110
*> \endverbatim
 
111
*>
 
112
*> \param[in] Q
 
113
*> \verbatim
 
114
*>          Q is DOUBLE PRECISION array, dimension (N**2)
 
115
*>         Contains the square eigenblocks from previous levels, the
 
116
*>         starting positions for blocks are given by QPTR.
 
117
*> \endverbatim
 
118
*>
 
119
*> \param[in] QPTR
 
120
*> \verbatim
 
121
*>          QPTR is INTEGER array, dimension (N+2)
 
122
*>         Contains a list of pointers which indicate where in Q an
 
123
*>         eigenblock is stored.  SQRT( QPTR(i+1) - QPTR(i) ) indicates
 
124
*>         the size of the block.
 
125
*> \endverbatim
 
126
*>
 
127
*> \param[out] Z
 
128
*> \verbatim
 
129
*>          Z is DOUBLE PRECISION array, dimension (N)
 
130
*>         On output this vector contains the updating vector (the last
 
131
*>         row of the first sub-eigenvector matrix and the first row of
 
132
*>         the second sub-eigenvector matrix).
 
133
*> \endverbatim
 
134
*>
 
135
*> \param[out] ZTEMP
 
136
*> \verbatim
 
137
*>          ZTEMP is DOUBLE PRECISION array, dimension (N)
 
138
*> \endverbatim
 
139
*>
 
140
*> \param[out] INFO
 
141
*> \verbatim
 
142
*>          INFO is INTEGER
 
143
*>          = 0:  successful exit.
 
144
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 
145
*> \endverbatim
 
146
*
 
147
*  Authors:
 
148
*  ========
 
149
*
 
150
*> \author Univ. of Tennessee 
 
151
*> \author Univ. of California Berkeley 
 
152
*> \author Univ. of Colorado Denver 
 
153
*> \author NAG Ltd. 
 
154
*
 
155
*> \date September 2012
 
156
*
 
157
*> \ingroup auxOTHERcomputational
 
158
*
 
159
*> \par Contributors:
 
160
*  ==================
 
161
*>
 
162
*> Jeff Rutter, Computer Science Division, University of California
 
163
*> at Berkeley, USA
 
164
*
 
165
*  =====================================================================
 
166
      SUBROUTINE DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
 
167
     $                   GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO )
 
168
*
 
169
*  -- LAPACK computational routine (version 3.4.2) --
 
170
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
171
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
172
*     September 2012
 
173
*
 
174
*     .. Scalar Arguments ..
 
175
      INTEGER            CURLVL, CURPBM, INFO, N, TLVLS
 
176
*     ..
 
177
*     .. Array Arguments ..
 
178
      INTEGER            GIVCOL( 2, * ), GIVPTR( * ), PERM( * ),
 
179
     $                   PRMPTR( * ), QPTR( * )
 
180
      DOUBLE PRECISION   GIVNUM( 2, * ), Q( * ), Z( * ), ZTEMP( * )
 
181
*     ..
 
182
*
 
183
*  =====================================================================
 
184
*
 
185
*     .. Parameters ..
 
186
      DOUBLE PRECISION   ZERO, HALF, ONE
 
187
      PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 )
 
188
*     ..
 
189
*     .. Local Scalars ..
 
190
      INTEGER            BSIZ1, BSIZ2, CURR, I, K, MID, PSIZ1, PSIZ2,
 
191
     $                   PTR, ZPTR1
 
192
*     ..
 
193
*     .. External Subroutines ..
 
194
      EXTERNAL           DCOPY, DGEMV, DROT, XERBLA
 
195
*     ..
 
196
*     .. Intrinsic Functions ..
 
197
      INTRINSIC          DBLE, INT, SQRT
 
198
*     ..
 
199
*     .. Executable Statements ..
 
200
*
 
201
*     Test the input parameters.
 
202
*
 
203
      INFO = 0
 
204
*
 
205
      IF( N.LT.0 ) THEN
 
206
         INFO = -1
 
207
      END IF
 
208
      IF( INFO.NE.0 ) THEN
 
209
         CALL XERBLA( 'DLAEDA', -INFO )
 
210
         RETURN
 
211
      END IF
 
212
*
 
213
*     Quick return if possible
 
214
*
 
215
      IF( N.EQ.0 )
 
216
     $   RETURN
 
217
*
 
218
*     Determine location of first number in second half.
 
219
*
 
220
      MID = N / 2 + 1
 
221
*
 
222
*     Gather last/first rows of appropriate eigenblocks into center of Z
 
223
*
 
224
      PTR = 1
 
225
*
 
226
*     Determine location of lowest level subproblem in the full storage
 
227
*     scheme
 
228
*
 
229
      CURR = PTR + CURPBM*2**CURLVL + 2**( CURLVL-1 ) - 1
 
230
*
 
231
*     Determine size of these matrices.  We add HALF to the value of
 
232
*     the SQRT in case the machine underestimates one of these square
 
233
*     roots.
 
234
*
 
235
      BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
 
236
      BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+1 ) ) ) )
 
237
      DO 10 K = 1, MID - BSIZ1 - 1
 
238
         Z( K ) = ZERO
 
239
   10 CONTINUE
 
240
      CALL DCOPY( BSIZ1, Q( QPTR( CURR )+BSIZ1-1 ), BSIZ1,
 
241
     $            Z( MID-BSIZ1 ), 1 )
 
242
      CALL DCOPY( BSIZ2, Q( QPTR( CURR+1 ) ), BSIZ2, Z( MID ), 1 )
 
243
      DO 20 K = MID + BSIZ2, N
 
244
         Z( K ) = ZERO
 
245
   20 CONTINUE
 
246
*
 
247
*     Loop through remaining levels 1 -> CURLVL applying the Givens
 
248
*     rotations and permutation and then multiplying the center matrices
 
249
*     against the current Z.
 
250
*
 
251
      PTR = 2**TLVLS + 1
 
252
      DO 70 K = 1, CURLVL - 1
 
253
         CURR = PTR + CURPBM*2**( CURLVL-K ) + 2**( CURLVL-K-1 ) - 1
 
254
         PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
 
255
         PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
 
256
         ZPTR1 = MID - PSIZ1
 
257
*
 
258
*       Apply Givens at CURR and CURR+1
 
259
*
 
260
         DO 30 I = GIVPTR( CURR ), GIVPTR( CURR+1 ) - 1
 
261
            CALL DROT( 1, Z( ZPTR1+GIVCOL( 1, I )-1 ), 1,
 
262
     $                 Z( ZPTR1+GIVCOL( 2, I )-1 ), 1, GIVNUM( 1, I ),
 
263
     $                 GIVNUM( 2, I ) )
 
264
   30    CONTINUE
 
265
         DO 40 I = GIVPTR( CURR+1 ), GIVPTR( CURR+2 ) - 1
 
266
            CALL DROT( 1, Z( MID-1+GIVCOL( 1, I ) ), 1,
 
267
     $                 Z( MID-1+GIVCOL( 2, I ) ), 1, GIVNUM( 1, I ),
 
268
     $                 GIVNUM( 2, I ) )
 
269
   40    CONTINUE
 
270
         PSIZ1 = PRMPTR( CURR+1 ) - PRMPTR( CURR )
 
271
         PSIZ2 = PRMPTR( CURR+2 ) - PRMPTR( CURR+1 )
 
272
         DO 50 I = 0, PSIZ1 - 1
 
273
            ZTEMP( I+1 ) = Z( ZPTR1+PERM( PRMPTR( CURR )+I )-1 )
 
274
   50    CONTINUE
 
275
         DO 60 I = 0, PSIZ2 - 1
 
276
            ZTEMP( PSIZ1+I+1 ) = Z( MID+PERM( PRMPTR( CURR+1 )+I )-1 )
 
277
   60    CONTINUE
 
278
*
 
279
*        Multiply Blocks at CURR and CURR+1
 
280
*
 
281
*        Determine size of these matrices.  We add HALF to the value of
 
282
*        the SQRT in case the machine underestimates one of these
 
283
*        square roots.
 
284
*
 
285
         BSIZ1 = INT( HALF+SQRT( DBLE( QPTR( CURR+1 )-QPTR( CURR ) ) ) )
 
286
         BSIZ2 = INT( HALF+SQRT( DBLE( QPTR( CURR+2 )-QPTR( CURR+
 
287
     $           1 ) ) ) )
 
288
         IF( BSIZ1.GT.0 ) THEN
 
289
            CALL DGEMV( 'T', BSIZ1, BSIZ1, ONE, Q( QPTR( CURR ) ),
 
290
     $                  BSIZ1, ZTEMP( 1 ), 1, ZERO, Z( ZPTR1 ), 1 )
 
291
         END IF
 
292
         CALL DCOPY( PSIZ1-BSIZ1, ZTEMP( BSIZ1+1 ), 1, Z( ZPTR1+BSIZ1 ),
 
293
     $               1 )
 
294
         IF( BSIZ2.GT.0 ) THEN
 
295
            CALL DGEMV( 'T', BSIZ2, BSIZ2, ONE, Q( QPTR( CURR+1 ) ),
 
296
     $                  BSIZ2, ZTEMP( PSIZ1+1 ), 1, ZERO, Z( MID ), 1 )
 
297
         END IF
 
298
         CALL DCOPY( PSIZ2-BSIZ2, ZTEMP( PSIZ1+BSIZ2+1 ), 1,
 
299
     $               Z( MID+BSIZ2 ), 1 )
 
300
*
 
301
         PTR = PTR + 2**( TLVLS-K )
 
302
   70 CONTINUE
 
303
*
 
304
      RETURN
 
305
*
 
306
*     End of DLAEDA
 
307
*
 
308
      END