← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/lapack/single/clarft.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/clarft.f
 
1
*> \brief \b CLARFT forms the triangular factor T of a block reflector H = I - vtvH
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*> \htmlonly
 
9
*> Download CLARFT + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarft.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarft.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarft.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
 
22
* 
 
23
*       .. Scalar Arguments ..
 
24
*       CHARACTER          DIRECT, STOREV
 
25
*       INTEGER            K, LDT, LDV, N
 
26
*       ..
 
27
*       .. Array Arguments ..
 
28
*       COMPLEX            T( LDT, * ), TAU( * ), V( LDV, * )
 
29
*       ..
 
30
*  
 
31
*
 
32
*> \par Purpose:
 
33
*  =============
 
34
*>
 
35
*> \verbatim
 
36
*>
 
37
*> CLARFT forms the triangular factor T of a complex block reflector H
 
38
*> of order n, which is defined as a product of k elementary reflectors.
 
39
*>
 
40
*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
 
41
*>
 
42
*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
 
43
*>
 
44
*> If STOREV = 'C', the vector which defines the elementary reflector
 
45
*> H(i) is stored in the i-th column of the array V, and
 
46
*>
 
47
*>    H  =  I - V * T * V**H
 
48
*>
 
49
*> If STOREV = 'R', the vector which defines the elementary reflector
 
50
*> H(i) is stored in the i-th row of the array V, and
 
51
*>
 
52
*>    H  =  I - V**H * T * V
 
53
*> \endverbatim
 
54
*
 
55
*  Arguments:
 
56
*  ==========
 
57
*
 
58
*> \param[in] DIRECT
 
59
*> \verbatim
 
60
*>          DIRECT is CHARACTER*1
 
61
*>          Specifies the order in which the elementary reflectors are
 
62
*>          multiplied to form the block reflector:
 
63
*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
 
64
*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
 
65
*> \endverbatim
 
66
*>
 
67
*> \param[in] STOREV
 
68
*> \verbatim
 
69
*>          STOREV is CHARACTER*1
 
70
*>          Specifies how the vectors which define the elementary
 
71
*>          reflectors are stored (see also Further Details):
 
72
*>          = 'C': columnwise
 
73
*>          = 'R': rowwise
 
74
*> \endverbatim
 
75
*>
 
76
*> \param[in] N
 
77
*> \verbatim
 
78
*>          N is INTEGER
 
79
*>          The order of the block reflector H. N >= 0.
 
80
*> \endverbatim
 
81
*>
 
82
*> \param[in] K
 
83
*> \verbatim
 
84
*>          K is INTEGER
 
85
*>          The order of the triangular factor T (= the number of
 
86
*>          elementary reflectors). K >= 1.
 
87
*> \endverbatim
 
88
*>
 
89
*> \param[in] V
 
90
*> \verbatim
 
91
*>          V is COMPLEX array, dimension
 
92
*>                               (LDV,K) if STOREV = 'C'
 
93
*>                               (LDV,N) if STOREV = 'R'
 
94
*>          The matrix V. See further details.
 
95
*> \endverbatim
 
96
*>
 
97
*> \param[in] LDV
 
98
*> \verbatim
 
99
*>          LDV is INTEGER
 
100
*>          The leading dimension of the array V.
 
101
*>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
 
102
*> \endverbatim
 
103
*>
 
104
*> \param[in] TAU
 
105
*> \verbatim
 
106
*>          TAU is COMPLEX array, dimension (K)
 
107
*>          TAU(i) must contain the scalar factor of the elementary
 
108
*>          reflector H(i).
 
109
*> \endverbatim
 
110
*>
 
111
*> \param[out] T
 
112
*> \verbatim
 
113
*>          T is COMPLEX array, dimension (LDT,K)
 
114
*>          The k by k triangular factor T of the block reflector.
 
115
*>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
 
116
*>          lower triangular. The rest of the array is not used.
 
117
*> \endverbatim
 
118
*>
 
119
*> \param[in] LDT
 
120
*> \verbatim
 
121
*>          LDT is INTEGER
 
122
*>          The leading dimension of the array T. LDT >= K.
 
123
*> \endverbatim
 
124
*
 
125
*  Authors:
 
126
*  ========
 
127
*
 
128
*> \author Univ. of Tennessee 
 
129
*> \author Univ. of California Berkeley 
 
130
*> \author Univ. of Colorado Denver 
 
131
*> \author NAG Ltd. 
 
132
*
 
133
*> \date September 2012
 
134
*
 
135
*> \ingroup complexOTHERauxiliary
 
136
*
 
137
*> \par Further Details:
 
138
*  =====================
 
139
*>
 
140
*> \verbatim
 
141
*>
 
142
*>  The shape of the matrix V and the storage of the vectors which define
 
143
*>  the H(i) is best illustrated by the following example with n = 5 and
 
144
*>  k = 3. The elements equal to 1 are not stored.
 
145
*>
 
146
*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
 
147
*>
 
148
*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
 
149
*>                   ( v1  1    )                     (     1 v2 v2 v2 )
 
150
*>                   ( v1 v2  1 )                     (        1 v3 v3 )
 
151
*>                   ( v1 v2 v3 )
 
152
*>                   ( v1 v2 v3 )
 
153
*>
 
154
*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
 
155
*>
 
156
*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
 
157
*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
 
158
*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
 
159
*>                   (     1 v3 )
 
160
*>                   (        1 )
 
161
*> \endverbatim
 
162
*>
 
163
*  =====================================================================
1
164
      SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
2
165
*
3
 
*  -- LAPACK auxiliary routine (version 2.0) --
4
 
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5
 
*     Courant Institute, Argonne National Lab, and Rice University
6
 
*     September 30, 1994
 
166
*  -- LAPACK auxiliary routine (version 3.4.2) --
 
167
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
168
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
169
*     September 2012
7
170
*
8
171
*     .. Scalar Arguments ..
9
172
      CHARACTER          DIRECT, STOREV
13
176
      COMPLEX            T( LDT, * ), TAU( * ), V( LDV, * )
14
177
*     ..
15
178
*
16
 
c
17
 
* $Id: clarft.f 19697 2010-10-29 16:57:34Z d3y133 $
18
 
c
19
 
*  Purpose
20
 
*  =======
21
 
*
22
 
*  CLARFT forms the triangular factor T of a complex block reflector H
23
 
*  of order n, which is defined as a product of k elementary reflectors.
24
 
*
25
 
*  If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
26
 
*
27
 
*  If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
28
 
*
29
 
*  If STOREV = 'C', the vector which defines the elementary reflector
30
 
*  H(i) is stored in the i-th column of the array V, and
31
 
*
32
 
*     H  =  I - V * T * V'
33
 
*
34
 
*  If STOREV = 'R', the vector which defines the elementary reflector
35
 
*  H(i) is stored in the i-th row of the array V, and
36
 
*
37
 
*     H  =  I - V' * T * V
38
 
*
39
 
*  Arguments
40
 
*  =========
41
 
*
42
 
*  DIRECT  (input) CHARACTER*1
43
 
*          Specifies the order in which the elementary reflectors are
44
 
*          multiplied to form the block reflector:
45
 
*          = 'F': H = H(1) H(2) . . . H(k) (Forward)
46
 
*          = 'B': H = H(k) . . . H(2) H(1) (Backward)
47
 
*
48
 
*  STOREV  (input) CHARACTER*1
49
 
*          Specifies how the vectors which define the elementary
50
 
*          reflectors are stored (see also Further Details):
51
 
*          = 'C': columnwise
52
 
*          = 'R': rowwise
53
 
*
54
 
*  N       (input) INTEGER
55
 
*          The order of the block reflector H. N >= 0.
56
 
*
57
 
*  K       (input) INTEGER
58
 
*          The order of the triangular factor T (= the number of
59
 
*          elementary reflectors). K >= 1.
60
 
*
61
 
*  V       (input/output) COMPLEX array, dimension
62
 
*                               (LDV,K) if STOREV = 'C'
63
 
*                               (LDV,N) if STOREV = 'R'
64
 
*          The matrix V. See further details.
65
 
*
66
 
*  LDV     (input) INTEGER
67
 
*          The leading dimension of the array V.
68
 
*          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
69
 
*
70
 
*  TAU     (input) COMPLEX array, dimension (K)
71
 
*          TAU(i) must contain the scalar factor of the elementary
72
 
*          reflector H(i).
73
 
*
74
 
*  T       (output) COMPLEX array, dimension (LDT,K)
75
 
*          The k by k triangular factor T of the block reflector.
76
 
*          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
77
 
*          lower triangular. The rest of the array is not used.
78
 
*
79
 
*  LDT     (input) INTEGER
80
 
*          The leading dimension of the array T. LDT >= K.
81
 
*
82
 
*  Further Details
83
 
*  ===============
84
 
*
85
 
*  The shape of the matrix V and the storage of the vectors which define
86
 
*  the H(i) is best illustrated by the following example with n = 5 and
87
 
*  k = 3. The elements equal to 1 are not stored; the corresponding
88
 
*  array elements are modified but restored on exit. The rest of the
89
 
*  array is not used.
90
 
*
91
 
*  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
92
 
*
93
 
*               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
94
 
*                   ( v1  1    )                     (     1 v2 v2 v2 )
95
 
*                   ( v1 v2  1 )                     (        1 v3 v3 )
96
 
*                   ( v1 v2 v3 )
97
 
*                   ( v1 v2 v3 )
98
 
*
99
 
*  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
100
 
*
101
 
*               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
102
 
*                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
103
 
*                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
104
 
*                   (     1 v3 )
105
 
*                   (        1 )
106
 
*
107
179
*  =====================================================================
108
180
*
109
181
*     .. Parameters ..
112
184
     $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
113
185
*     ..
114
186
*     .. Local Scalars ..
115
 
      INTEGER            I, J
116
 
      COMPLEX            VII
 
187
      INTEGER            I, J, PREVLASTV, LASTV
117
188
*     ..
118
189
*     .. External Subroutines ..
119
190
      EXTERNAL           CGEMV, CLACGV, CTRMV
130
201
     $   RETURN
131
202
*
132
203
      IF( LSAME( DIRECT, 'F' ) ) THEN
133
 
         DO 20 I = 1, K
 
204
         PREVLASTV = N
 
205
         DO I = 1, K
 
206
            PREVLASTV = MAX( PREVLASTV, I )
134
207
            IF( TAU( I ).EQ.ZERO ) THEN
135
208
*
136
209
*              H(i)  =  I
137
210
*
138
 
               DO 10 J = 1, I
 
211
               DO J = 1, I
139
212
                  T( J, I ) = ZERO
140
 
   10          CONTINUE
 
213
               END DO
141
214
            ELSE
142
215
*
143
216
*              general case
144
217
*
145
 
               VII = V( I, I )
146
 
               V( I, I ) = ONE
147
218
               IF( LSAME( STOREV, 'C' ) ) THEN
148
 
*
149
 
*                 T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i)
150
 
*
151
 
                  CALL CGEMV( 'Conjugate transpose', N-I+1, I-1,
152
 
     $                        -TAU( I ), V( I, 1 ), LDV, V( I, I ), 1,
153
 
     $                        ZERO, T( 1, I ), 1 )
 
219
*                 Skip any trailing zeros.
 
220
                  DO LASTV = N, I+1, -1
 
221
                     IF( V( LASTV, I ).NE.ZERO ) EXIT
 
222
                  END DO
 
223
                  DO J = 1, I-1
 
224
                     T( J, I ) = -TAU( I ) * CONJG( V( I , J ) )
 
225
                  END DO                     
 
226
                  J = MIN( LASTV, PREVLASTV )
 
227
*
 
228
*                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i)
 
229
*
 
230
                  CALL CGEMV( 'Conjugate transpose', J-I, I-1,
 
231
     $                        -TAU( I ), V( I+1, 1 ), LDV, 
 
232
     $                        V( I+1, I ), 1,
 
233
     $                        ONE, T( 1, I ), 1 )
154
234
               ELSE
155
 
*
156
 
*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)'
157
 
*
158
 
                  IF( I.LT.N )
159
 
     $               CALL CLACGV( N-I, V( I, I+1 ), LDV )
160
 
                  CALL CGEMV( 'No transpose', I-1, N-I+1, -TAU( I ),
161
 
     $                        V( 1, I ), LDV, V( I, I ), LDV, ZERO,
162
 
     $                        T( 1, I ), 1 )
163
 
                  IF( I.LT.N )
164
 
     $               CALL CLACGV( N-I, V( I, I+1 ), LDV )
 
235
*                 Skip any trailing zeros.
 
236
                  DO LASTV = N, I+1, -1
 
237
                     IF( V( I, LASTV ).NE.ZERO ) EXIT
 
238
                  END DO
 
239
                  DO J = 1, I-1
 
240
                     T( J, I ) = -TAU( I ) * V( J , I )
 
241
                  END DO                     
 
242
                  J = MIN( LASTV, PREVLASTV )
 
243
*
 
244
*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H
 
245
*
 
246
                  CALL CGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ),
 
247
     $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV,
 
248
     $                        ONE, T( 1, I ), LDT )                  
165
249
               END IF
166
 
               V( I, I ) = VII
167
250
*
168
251
*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
169
252
*
170
253
               CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
171
254
     $                     LDT, T( 1, I ), 1 )
172
255
               T( I, I ) = TAU( I )
 
256
               IF( I.GT.1 ) THEN
 
257
                  PREVLASTV = MAX( PREVLASTV, LASTV )
 
258
               ELSE
 
259
                  PREVLASTV = LASTV
 
260
               END IF
173
261
            END IF
174
 
   20    CONTINUE
 
262
         END DO
175
263
      ELSE
176
 
         DO 40 I = K, 1, -1
 
264
         PREVLASTV = 1
 
265
         DO I = K, 1, -1
177
266
            IF( TAU( I ).EQ.ZERO ) THEN
178
267
*
179
268
*              H(i)  =  I
180
269
*
181
 
               DO 30 J = I, K
 
270
               DO J = I, K
182
271
                  T( J, I ) = ZERO
183
 
   30          CONTINUE
 
272
               END DO
184
273
            ELSE
185
274
*
186
275
*              general case
187
276
*
188
277
               IF( I.LT.K ) THEN
189
278
                  IF( LSAME( STOREV, 'C' ) ) THEN
190
 
                     VII = V( N-K+I, I )
191
 
                     V( N-K+I, I ) = ONE
192
 
*
193
 
*                    T(i+1:k,i) :=
194
 
*                            - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i)
195
 
*
196
 
                     CALL CGEMV( 'Conjugate transpose', N-K+I, K-I,
197
 
     $                           -TAU( I ), V( 1, I+1 ), LDV, V( 1, I ),
198
 
     $                           1, ZERO, T( I+1, I ), 1 )
199
 
                     V( N-K+I, I ) = VII
 
279
*                    Skip any leading zeros.
 
280
                     DO LASTV = 1, I-1
 
281
                        IF( V( LASTV, I ).NE.ZERO ) EXIT
 
282
                     END DO
 
283
                     DO J = I+1, K
 
284
                        T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) )
 
285
                     END DO                        
 
286
                     J = MAX( LASTV, PREVLASTV )
 
287
*
 
288
*                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i)
 
289
*
 
290
                     CALL CGEMV( 'Conjugate transpose', N-K+I-J, K-I,
 
291
     $                           -TAU( I ), V( J, I+1 ), LDV, V( J, I ),
 
292
     $                           1, ONE, T( I+1, I ), 1 )
200
293
                  ELSE
201
 
                     VII = V( I, N-K+I )
202
 
                     V( I, N-K+I ) = ONE
203
 
*
204
 
*                    T(i+1:k,i) :=
205
 
*                            - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)'
206
 
*
207
 
                     CALL CLACGV( N-K+I-1, V( I, 1 ), LDV )
208
 
                     CALL CGEMV( 'No transpose', K-I, N-K+I, -TAU( I ),
209
 
     $                           V( I+1, 1 ), LDV, V( I, 1 ), LDV, ZERO,
210
 
     $                           T( I+1, I ), 1 )
211
 
                     CALL CLACGV( N-K+I-1, V( I, 1 ), LDV )
212
 
                     V( I, N-K+I ) = VII
 
294
*                    Skip any leading zeros.
 
295
                     DO LASTV = 1, I-1
 
296
                        IF( V( I, LASTV ).NE.ZERO ) EXIT
 
297
                     END DO
 
298
                     DO J = I+1, K
 
299
                        T( J, I ) = -TAU( I ) * V( J, N-K+I )
 
300
                     END DO                      
 
301
                     J = MAX( LASTV, PREVLASTV )
 
302
*
 
303
*                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H
 
304
*
 
305
                     CALL CGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ),
 
306
     $                           V( I+1, J ), LDV, V( I, J ), LDV,
 
307
     $                           ONE, T( I+1, I ), LDT )                     
213
308
                  END IF
214
309
*
215
310
*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
216
311
*
217
312
                  CALL CTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
218
313
     $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
 
314
                  IF( I.GT.1 ) THEN
 
315
                     PREVLASTV = MIN( PREVLASTV, LASTV )
 
316
                  ELSE
 
317
                     PREVLASTV = LASTV
 
318
                  END IF
219
319
               END IF
220
320
               T( I, I ) = TAU( I )
221
321
            END IF
222
 
   40    CONTINUE
 
322
         END DO
223
323
      END IF
224
324
      RETURN
225
325
*