← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/lapack/double/dlarfb.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/double/dlarfb.f
 
1
*> \brief \b DLARFB applies a block reflector or its transpose to a general rectangular matrix.
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*> \htmlonly
 
9
*> Download DLARFB + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfb.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfb.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfb.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
 
22
*                          T, LDT, C, LDC, WORK, LDWORK )
 
23
* 
 
24
*       .. Scalar Arguments ..
 
25
*       CHARACTER          DIRECT, SIDE, STOREV, TRANS
 
26
*       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
 
27
*       ..
 
28
*       .. Array Arguments ..
 
29
*       DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),
 
30
*      $                   WORK( LDWORK, * )
 
31
*       ..
 
32
*  
 
33
*
 
34
*> \par Purpose:
 
35
*  =============
 
36
*>
 
37
*> \verbatim
 
38
*>
 
39
*> DLARFB applies a real block reflector H or its transpose H**T to a
 
40
*> real m by n matrix C, from either the left or the right.
 
41
*> \endverbatim
 
42
*
 
43
*  Arguments:
 
44
*  ==========
 
45
*
 
46
*> \param[in] SIDE
 
47
*> \verbatim
 
48
*>          SIDE is CHARACTER*1
 
49
*>          = 'L': apply H or H**T from the Left
 
50
*>          = 'R': apply H or H**T from the Right
 
51
*> \endverbatim
 
52
*>
 
53
*> \param[in] TRANS
 
54
*> \verbatim
 
55
*>          TRANS is CHARACTER*1
 
56
*>          = 'N': apply H (No transpose)
 
57
*>          = 'T': apply H**T (Transpose)
 
58
*> \endverbatim
 
59
*>
 
60
*> \param[in] DIRECT
 
61
*> \verbatim
 
62
*>          DIRECT is CHARACTER*1
 
63
*>          Indicates how H is formed from a product of elementary
 
64
*>          reflectors
 
65
*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
 
66
*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
 
67
*> \endverbatim
 
68
*>
 
69
*> \param[in] STOREV
 
70
*> \verbatim
 
71
*>          STOREV is CHARACTER*1
 
72
*>          Indicates how the vectors which define the elementary
 
73
*>          reflectors are stored:
 
74
*>          = 'C': Columnwise
 
75
*>          = 'R': Rowwise
 
76
*> \endverbatim
 
77
*>
 
78
*> \param[in] M
 
79
*> \verbatim
 
80
*>          M is INTEGER
 
81
*>          The number of rows of the matrix C.
 
82
*> \endverbatim
 
83
*>
 
84
*> \param[in] N
 
85
*> \verbatim
 
86
*>          N is INTEGER
 
87
*>          The number of columns of the matrix C.
 
88
*> \endverbatim
 
89
*>
 
90
*> \param[in] K
 
91
*> \verbatim
 
92
*>          K is INTEGER
 
93
*>          The order of the matrix T (= the number of elementary
 
94
*>          reflectors whose product defines the block reflector).
 
95
*> \endverbatim
 
96
*>
 
97
*> \param[in] V
 
98
*> \verbatim
 
99
*>          V is DOUBLE PRECISION array, dimension
 
100
*>                                (LDV,K) if STOREV = 'C'
 
101
*>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
 
102
*>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
 
103
*>          The matrix V. See Further Details.
 
104
*> \endverbatim
 
105
*>
 
106
*> \param[in] LDV
 
107
*> \verbatim
 
108
*>          LDV is INTEGER
 
109
*>          The leading dimension of the array V.
 
110
*>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
 
111
*>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
 
112
*>          if STOREV = 'R', LDV >= K.
 
113
*> \endverbatim
 
114
*>
 
115
*> \param[in] T
 
116
*> \verbatim
 
117
*>          T is DOUBLE PRECISION array, dimension (LDT,K)
 
118
*>          The triangular k by k matrix T in the representation of the
 
119
*>          block reflector.
 
120
*> \endverbatim
 
121
*>
 
122
*> \param[in] LDT
 
123
*> \verbatim
 
124
*>          LDT is INTEGER
 
125
*>          The leading dimension of the array T. LDT >= K.
 
126
*> \endverbatim
 
127
*>
 
128
*> \param[in,out] C
 
129
*> \verbatim
 
130
*>          C is DOUBLE PRECISION array, dimension (LDC,N)
 
131
*>          On entry, the m by n matrix C.
 
132
*>          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.
 
133
*> \endverbatim
 
134
*>
 
135
*> \param[in] LDC
 
136
*> \verbatim
 
137
*>          LDC is INTEGER
 
138
*>          The leading dimension of the array C. LDC >= max(1,M).
 
139
*> \endverbatim
 
140
*>
 
141
*> \param[out] WORK
 
142
*> \verbatim
 
143
*>          WORK is DOUBLE PRECISION array, dimension (LDWORK,K)
 
144
*> \endverbatim
 
145
*>
 
146
*> \param[in] LDWORK
 
147
*> \verbatim
 
148
*>          LDWORK is INTEGER
 
149
*>          The leading dimension of the array WORK.
 
150
*>          If SIDE = 'L', LDWORK >= max(1,N);
 
151
*>          if SIDE = 'R', LDWORK >= max(1,M).
 
152
*> \endverbatim
 
153
*
 
154
*  Authors:
 
155
*  ========
 
156
*
 
157
*> \author Univ. of Tennessee 
 
158
*> \author Univ. of California Berkeley 
 
159
*> \author Univ. of Colorado Denver 
 
160
*> \author NAG Ltd. 
 
161
*
 
162
*> \date September 2012
 
163
*
 
164
*> \ingroup doubleOTHERauxiliary
 
165
*
 
166
*> \par Further Details:
 
167
*  =====================
 
168
*>
 
169
*> \verbatim
 
170
*>
 
171
*>  The shape of the matrix V and the storage of the vectors which define
 
172
*>  the H(i) is best illustrated by the following example with n = 5 and
 
173
*>  k = 3. The elements equal to 1 are not stored; the corresponding
 
174
*>  array elements are modified but restored on exit. The rest of the
 
175
*>  array is not used.
 
176
*>
 
177
*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
 
178
*>
 
179
*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
 
180
*>                   ( v1  1    )                     (     1 v2 v2 v2 )
 
181
*>                   ( v1 v2  1 )                     (        1 v3 v3 )
 
182
*>                   ( v1 v2 v3 )
 
183
*>                   ( v1 v2 v3 )
 
184
*>
 
185
*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
 
186
*>
 
187
*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
 
188
*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
 
189
*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
 
190
*>                   (     1 v3 )
 
191
*>                   (        1 )
 
192
*> \endverbatim
 
193
*>
 
194
*  =====================================================================
1
195
      SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
2
196
     $                   T, LDT, C, LDC, WORK, LDWORK )
3
197
*
4
 
*  -- LAPACK auxiliary routine (version 2.0) --
5
 
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6
 
*     Courant Institute, Argonne National Lab, and Rice University
7
 
*     February 29, 1992
 
198
*  -- LAPACK auxiliary routine (version 3.4.2) --
 
199
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
200
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
201
*     September 2012
8
202
*
9
203
*     .. Scalar Arguments ..
10
204
      CHARACTER          DIRECT, SIDE, STOREV, TRANS
15
209
     $                   WORK( LDWORK, * )
16
210
*     ..
17
211
*
18
 
c
19
 
* $Id: dlarfb.f 19697 2010-10-29 16:57:34Z d3y133 $
20
 
c
21
 
*  Purpose
22
 
*  =======
23
 
*
24
 
*  DLARFB applies a real block reflector H or its transpose H' to a
25
 
*  real m by n matrix C, from either the left or the right.
26
 
*
27
 
*  Arguments
28
 
*  =========
29
 
*
30
 
*  SIDE    (input) CHARACTER*1
31
 
*          = 'L': apply H or H' from the Left
32
 
*          = 'R': apply H or H' from the Right
33
 
*
34
 
*  TRANS   (input) CHARACTER*1
35
 
*          = 'N': apply H (No transpose)
36
 
*          = 'T': apply H' (Transpose)
37
 
*
38
 
*  DIRECT  (input) CHARACTER*1
39
 
*          Indicates how H is formed from a product of elementary
40
 
*          reflectors
41
 
*          = 'F': H = H(1) H(2) . . . H(k) (Forward)
42
 
*          = 'B': H = H(k) . . . H(2) H(1) (Backward)
43
 
*
44
 
*  STOREV  (input) CHARACTER*1
45
 
*          Indicates how the vectors which define the elementary
46
 
*          reflectors are stored:
47
 
*          = 'C': Columnwise
48
 
*          = 'R': Rowwise
49
 
*
50
 
*  M       (input) INTEGER
51
 
*          The number of rows of the matrix C.
52
 
*
53
 
*  N       (input) INTEGER
54
 
*          The number of columns of the matrix C.
55
 
*
56
 
*  K       (input) INTEGER
57
 
*          The order of the matrix T (= the number of elementary
58
 
*          reflectors whose product defines the block reflector).
59
 
*
60
 
*  V       (input) DOUBLE PRECISION array, dimension
61
 
*                                (LDV,K) if STOREV = 'C'
62
 
*                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
63
 
*                                (LDV,N) if STOREV = 'R' and SIDE = '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' and SIDE = 'L', LDV >= max(1,M);
69
 
*          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
70
 
*          if STOREV = 'R', LDV >= K.
71
 
*
72
 
*  T       (input) DOUBLE PRECISION array, dimension (LDT,K)
73
 
*          The triangular k by k matrix T in the representation of the
74
 
*          block reflector.
75
 
*
76
 
*  LDT     (input) INTEGER
77
 
*          The leading dimension of the array T. LDT >= K.
78
 
*
79
 
*  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
80
 
*          On entry, the m by n matrix C.
81
 
*          On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
82
 
*
83
 
*  LDC     (input) INTEGER
84
 
*          The leading dimension of the array C. LDA >= max(1,M).
85
 
*
86
 
*  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
87
 
*
88
 
*  LDWORK  (input) INTEGER
89
 
*          The leading dimension of the array WORK.
90
 
*          If SIDE = 'L', LDWORK >= max(1,N);
91
 
*          if SIDE = 'R', LDWORK >= max(1,M).
92
 
*
93
212
*  =====================================================================
94
213
*
95
214
*     .. Parameters ..
98
217
*     ..
99
218
*     .. Local Scalars ..
100
219
      CHARACTER          TRANST
101
 
      INTEGER            I, J
 
220
      INTEGER            I, J, LASTV, LASTC, lastv2
102
221
*     ..
103
222
*     .. External Functions ..
104
223
      LOGICAL            LSAME
105
 
      EXTERNAL           LSAME
 
224
      INTEGER            ILADLR, ILADLC
 
225
      EXTERNAL           LSAME, ILADLR, ILADLC
106
226
*     ..
107
227
*     .. External Subroutines ..
108
228
      EXTERNAL           DCOPY, DGEMM, DTRMM
130
250
*
131
251
            IF( LSAME( SIDE, 'L' ) ) THEN
132
252
*
133
 
*              Form  H * C  or  H' * C  where  C = ( C1 )
134
 
*                                                  ( C2 )
135
 
*
136
 
*              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)
137
 
*
138
 
*              W := C1'
 
253
*              Form  H * C  or  H**T * C  where  C = ( C1 )
 
254
*                                                    ( C2 )
 
255
*
 
256
               LASTV = MAX( K, ILADLR( M, K, V, LDV ) )
 
257
               LASTC = ILADLC( LASTV, N, C, LDC )
 
258
*
 
259
*              W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
 
260
*
 
261
*              W := C1**T
139
262
*
140
263
               DO 10 J = 1, K
141
 
                  CALL DCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
 
264
                  CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
142
265
   10          CONTINUE
143
266
*
144
267
*              W := W * V1
145
268
*
146
 
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
147
 
     $                     K, ONE, V, LDV, WORK, LDWORK )
148
 
               IF( M.GT.K ) THEN
149
 
*
150
 
*                 W := W + C2'*V2
151
 
*
152
 
                  CALL DGEMM( 'Transpose', 'No transpose', N, K, M-K,
153
 
     $                        ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV,
154
 
     $                        ONE, WORK, LDWORK )
155
 
               END IF
156
 
*
157
 
*              W := W * T'  or  W * T
158
 
*
159
 
               CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
160
 
     $                     ONE, T, LDT, WORK, LDWORK )
161
 
*
162
 
*              C := C - V * W'
163
 
*
164
 
               IF( M.GT.K ) THEN
165
 
*
166
 
*                 C2 := C2 - V2 * W'
167
 
*
168
 
                  CALL DGEMM( 'No transpose', 'Transpose', M-K, N, K,
169
 
     $                        -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE,
170
 
     $                        C( K+1, 1 ), LDC )
171
 
               END IF
172
 
*
173
 
*              W := W * V1'
174
 
*
175
 
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', N, K,
176
 
     $                     ONE, V, LDV, WORK, LDWORK )
177
 
*
178
 
*              C1 := C1 - W'
 
269
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
 
270
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
 
271
               IF( LASTV.GT.K ) THEN
 
272
*
 
273
*                 W := W + C2**T *V2
 
274
*
 
275
                  CALL DGEMM( 'Transpose', 'No transpose',
 
276
     $                 LASTC, K, LASTV-K,
 
277
     $                 ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV,
 
278
     $                 ONE, WORK, LDWORK )
 
279
               END IF
 
280
*
 
281
*              W := W * T**T  or  W * T
 
282
*
 
283
               CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
 
284
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 
285
*
 
286
*              C := C - V * W**T
 
287
*
 
288
               IF( LASTV.GT.K ) THEN
 
289
*
 
290
*                 C2 := C2 - V2 * W**T
 
291
*
 
292
                  CALL DGEMM( 'No transpose', 'Transpose',
 
293
     $                 LASTV-K, LASTC, K,
 
294
     $                 -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE,
 
295
     $                 C( K+1, 1 ), LDC )
 
296
               END IF
 
297
*
 
298
*              W := W * V1**T
 
299
*
 
300
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
 
301
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
 
302
*
 
303
*              C1 := C1 - W**T
179
304
*
180
305
               DO 30 J = 1, K
181
 
                  DO 20 I = 1, N
 
306
                  DO 20 I = 1, LASTC
182
307
                     C( J, I ) = C( J, I ) - WORK( I, J )
183
308
   20             CONTINUE
184
309
   30          CONTINUE
185
310
*
186
311
            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
187
312
*
188
 
*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 
313
*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
 
314
*
 
315
               LASTV = MAX( K, ILADLR( N, K, V, LDV ) )
 
316
               LASTC = ILADLR( M, LASTV, C, LDC )
189
317
*
190
318
*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
191
319
*
192
320
*              W := C1
193
321
*
194
322
               DO 40 J = 1, K
195
 
                  CALL DCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
 
323
                  CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
196
324
   40          CONTINUE
197
325
*
198
326
*              W := W * V1
199
327
*
200
 
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
201
 
     $                     K, ONE, V, LDV, WORK, LDWORK )
202
 
               IF( N.GT.K ) THEN
 
328
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
 
329
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
 
330
               IF( LASTV.GT.K ) THEN
203
331
*
204
332
*                 W := W + C2 * V2
205
333
*
206
 
                  CALL DGEMM( 'No transpose', 'No transpose', M, K, N-K,
207
 
     $                        ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
208
 
     $                        ONE, WORK, LDWORK )
209
 
               END IF
210
 
*
211
 
*              W := W * T  or  W * T'
212
 
*
213
 
               CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
214
 
     $                     ONE, T, LDT, WORK, LDWORK )
215
 
*
216
 
*              C := C - W * V'
217
 
*
218
 
               IF( N.GT.K ) THEN
219
 
*
220
 
*                 C2 := C2 - W * V2'
221
 
*
222
 
                  CALL DGEMM( 'No transpose', 'Transpose', M, N-K, K,
223
 
     $                        -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE,
224
 
     $                        C( 1, K+1 ), LDC )
225
 
               END IF
226
 
*
227
 
*              W := W * V1'
228
 
*
229
 
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', M, K,
230
 
     $                     ONE, V, LDV, WORK, LDWORK )
 
334
                  CALL DGEMM( 'No transpose', 'No transpose',
 
335
     $                 LASTC, K, LASTV-K,
 
336
     $                 ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
 
337
     $                 ONE, WORK, LDWORK )
 
338
               END IF
 
339
*
 
340
*              W := W * T  or  W * T**T
 
341
*
 
342
               CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
 
343
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 
344
*
 
345
*              C := C - W * V**T
 
346
*
 
347
               IF( LASTV.GT.K ) THEN
 
348
*
 
349
*                 C2 := C2 - W * V2**T
 
350
*
 
351
                  CALL DGEMM( 'No transpose', 'Transpose',
 
352
     $                 LASTC, LASTV-K, K,
 
353
     $                 -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE,
 
354
     $                 C( 1, K+1 ), LDC )
 
355
               END IF
 
356
*
 
357
*              W := W * V1**T
 
358
*
 
359
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
 
360
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
231
361
*
232
362
*              C1 := C1 - W
233
363
*
234
364
               DO 60 J = 1, K
235
 
                  DO 50 I = 1, M
 
365
                  DO 50 I = 1, LASTC
236
366
                     C( I, J ) = C( I, J ) - WORK( I, J )
237
367
   50             CONTINUE
238
368
   60          CONTINUE
246
376
*
247
377
            IF( LSAME( SIDE, 'L' ) ) THEN
248
378
*
249
 
*              Form  H * C  or  H' * C  where  C = ( C1 )
250
 
*                                                  ( C2 )
251
 
*
252
 
*              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)
253
 
*
254
 
*              W := C2'
 
379
*              Form  H * C  or  H**T * C  where  C = ( C1 )
 
380
*                                                    ( C2 )
 
381
*
 
382
               LASTC = ILADLC( M, N, C, LDC )
 
383
*
 
384
*              W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
 
385
*
 
386
*              W := C2**T
255
387
*
256
388
               DO 70 J = 1, K
257
 
                  CALL DCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
 
389
                  CALL DCOPY( LASTC, C( M-K+J, 1 ), LDC,
 
390
     $                 WORK( 1, J ), 1 )
258
391
   70          CONTINUE
259
392
*
260
393
*              W := W * V2
261
394
*
262
 
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
263
 
     $                     K, ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
264
 
               IF( M.GT.K ) THEN
265
 
*
266
 
*                 W := W + C1'*V1
267
 
*
268
 
                  CALL DGEMM( 'Transpose', 'No transpose', N, K, M-K,
269
 
     $                        ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
270
 
               END IF
271
 
*
272
 
*              W := W * T'  or  W * T
273
 
*
274
 
               CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
275
 
     $                     ONE, T, LDT, WORK, LDWORK )
276
 
*
277
 
*              C := C - V * W'
278
 
*
279
 
               IF( M.GT.K ) THEN
280
 
*
281
 
*                 C1 := C1 - V1 * W'
282
 
*
283
 
                  CALL DGEMM( 'No transpose', 'Transpose', M-K, N, K,
284
 
     $                        -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC )
285
 
               END IF
286
 
*
287
 
*              W := W * V2'
288
 
*
289
 
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', N, K,
290
 
     $                     ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
291
 
*
292
 
*              C2 := C2 - W'
 
395
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
 
396
     $              LASTC, K, ONE, V( M-K+1, 1 ), LDV,
 
397
     $              WORK, LDWORK )
 
398
               IF( M.GT.K ) THEN
 
399
*
 
400
*                 W := W + C1**T*V1
 
401
*
 
402
                  CALL DGEMM( 'Transpose', 'No transpose',
 
403
     $                 LASTC, K, M-K, ONE, C, LDC, V, LDV,
 
404
     $                 ONE, WORK, LDWORK )
 
405
               END IF
 
406
*
 
407
*              W := W * T**T  or  W * T
 
408
*
 
409
               CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
 
410
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 
411
*
 
412
*              C := C - V * W**T
 
413
*
 
414
               IF( M.GT.K ) THEN
 
415
*
 
416
*                 C1 := C1 - V1 * W**T
 
417
*
 
418
                  CALL DGEMM( 'No transpose', 'Transpose',
 
419
     $                 M-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
 
420
     $                 ONE, C, LDC )
 
421
               END IF
 
422
*
 
423
*              W := W * V2**T
 
424
*
 
425
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
 
426
     $              LASTC, K, ONE, V( M-K+1, 1 ), LDV,
 
427
     $              WORK, LDWORK )
 
428
*
 
429
*              C2 := C2 - W**T
293
430
*
294
431
               DO 90 J = 1, K
295
 
                  DO 80 I = 1, N
296
 
                     C( M-K+J, I ) = C( M-K+J, I ) - WORK( I, J )
 
432
                  DO 80 I = 1, LASTC
 
433
                     C( M-K+J, I ) = C( M-K+J, I ) - WORK(I, J)
297
434
   80             CONTINUE
298
435
   90          CONTINUE
299
436
*
300
437
            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
301
438
*
302
 
*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 
439
*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
 
440
*
 
441
               LASTC = ILADLR( M, N, C, LDC )
303
442
*
304
443
*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
305
444
*
306
445
*              W := C2
307
446
*
308
447
               DO 100 J = 1, K
309
 
                  CALL DCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
 
448
                  CALL DCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
310
449
  100          CONTINUE
311
450
*
312
451
*              W := W * V2
313
452
*
314
 
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
315
 
     $                     K, ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
 
453
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
 
454
     $              LASTC, K, ONE, V( N-K+1, 1 ), LDV,
 
455
     $              WORK, LDWORK )
316
456
               IF( N.GT.K ) THEN
317
457
*
318
458
*                 W := W + C1 * V1
319
459
*
320
 
                  CALL DGEMM( 'No transpose', 'No transpose', M, K, N-K,
321
 
     $                        ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
 
460
                  CALL DGEMM( 'No transpose', 'No transpose',
 
461
     $                 LASTC, K, N-K, ONE, C, LDC, V, LDV,
 
462
     $                 ONE, WORK, LDWORK )
322
463
               END IF
323
464
*
324
 
*              W := W * T  or  W * T'
325
 
*
326
 
               CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
327
 
     $                     ONE, T, LDT, WORK, LDWORK )
328
 
*
329
 
*              C := C - W * V'
 
465
*              W := W * T  or  W * T**T
 
466
*
 
467
               CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
 
468
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 
469
*
 
470
*              C := C - W * V**T
330
471
*
331
472
               IF( N.GT.K ) THEN
332
473
*
333
 
*                 C1 := C1 - W * V1'
 
474
*                 C1 := C1 - W * V1**T
334
475
*
335
 
                  CALL DGEMM( 'No transpose', 'Transpose', M, N-K, K,
336
 
     $                        -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
 
476
                  CALL DGEMM( 'No transpose', 'Transpose',
 
477
     $                 LASTC, N-K, K, -ONE, WORK, LDWORK, V, LDV,
 
478
     $                 ONE, C, LDC )
337
479
               END IF
338
480
*
339
 
*              W := W * V2'
 
481
*              W := W * V2**T
340
482
*
341
 
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', M, K,
342
 
     $                     ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
 
483
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
 
484
     $              LASTC, K, ONE, V( N-K+1, 1 ), LDV,
 
485
     $              WORK, LDWORK )
343
486
*
344
487
*              C2 := C2 - W
345
488
*
346
489
               DO 120 J = 1, K
347
 
                  DO 110 I = 1, M
348
 
                     C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
 
490
                  DO 110 I = 1, LASTC
 
491
                     C( I, N-K+J ) = C( I, N-K+J ) - WORK(I, J)
349
492
  110             CONTINUE
350
493
  120          CONTINUE
351
494
            END IF
360
503
*
361
504
            IF( LSAME( SIDE, 'L' ) ) THEN
362
505
*
363
 
*              Form  H * C  or  H' * C  where  C = ( C1 )
364
 
*                                                  ( C2 )
365
 
*
366
 
*              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)
367
 
*
368
 
*              W := C1'
 
506
*              Form  H * C  or  H**T * C  where  C = ( C1 )
 
507
*                                                    ( C2 )
 
508
*
 
509
               LASTV = MAX( K, ILADLC( K, M, V, LDV ) )
 
510
               LASTC = ILADLC( LASTV, N, C, LDC )
 
511
*
 
512
*              W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
 
513
*
 
514
*              W := C1**T
369
515
*
370
516
               DO 130 J = 1, K
371
 
                  CALL DCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
 
517
                  CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
372
518
  130          CONTINUE
373
519
*
374
 
*              W := W * V1'
375
 
*
376
 
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', N, K,
377
 
     $                     ONE, V, LDV, WORK, LDWORK )
378
 
               IF( M.GT.K ) THEN
379
 
*
380
 
*                 W := W + C2'*V2'
381
 
*
382
 
                  CALL DGEMM( 'Transpose', 'Transpose', N, K, M-K, ONE,
383
 
     $                        C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, ONE,
384
 
     $                        WORK, LDWORK )
 
520
*              W := W * V1**T
 
521
*
 
522
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
 
523
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
 
524
               IF( LASTV.GT.K ) THEN
 
525
*
 
526
*                 W := W + C2**T*V2**T
 
527
*
 
528
                  CALL DGEMM( 'Transpose', 'Transpose',
 
529
     $                 LASTC, K, LASTV-K,
 
530
     $                 ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV,
 
531
     $                 ONE, WORK, LDWORK )
385
532
               END IF
386
533
*
387
 
*              W := W * T'  or  W * T
388
 
*
389
 
               CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
390
 
     $                     ONE, T, LDT, WORK, LDWORK )
391
 
*
392
 
*              C := C - V' * W'
393
 
*
394
 
               IF( M.GT.K ) THEN
395
 
*
396
 
*                 C2 := C2 - V2' * W'
397
 
*
398
 
                  CALL DGEMM( 'Transpose', 'Transpose', M-K, N, K, -ONE,
399
 
     $                        V( 1, K+1 ), LDV, WORK, LDWORK, ONE,
400
 
     $                        C( K+1, 1 ), LDC )
 
534
*              W := W * T**T  or  W * T
 
535
*
 
536
               CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
 
537
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 
538
*
 
539
*              C := C - V**T * W**T
 
540
*
 
541
               IF( LASTV.GT.K ) THEN
 
542
*
 
543
*                 C2 := C2 - V2**T * W**T
 
544
*
 
545
                  CALL DGEMM( 'Transpose', 'Transpose',
 
546
     $                 LASTV-K, LASTC, K,
 
547
     $                 -ONE, V( 1, K+1 ), LDV, WORK, LDWORK,
 
548
     $                 ONE, C( K+1, 1 ), LDC )
401
549
               END IF
402
550
*
403
551
*              W := W * V1
404
552
*
405
 
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
406
 
     $                     K, ONE, V, LDV, WORK, LDWORK )
 
553
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
 
554
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
407
555
*
408
 
*              C1 := C1 - W'
 
556
*              C1 := C1 - W**T
409
557
*
410
558
               DO 150 J = 1, K
411
 
                  DO 140 I = 1, N
 
559
                  DO 140 I = 1, LASTC
412
560
                     C( J, I ) = C( J, I ) - WORK( I, J )
413
561
  140             CONTINUE
414
562
  150          CONTINUE
415
563
*
416
564
            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
417
565
*
418
 
*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
419
 
*
420
 
*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)
 
566
*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
 
567
*
 
568
               LASTV = MAX( K, ILADLC( K, N, V, LDV ) )
 
569
               LASTC = ILADLR( M, LASTV, C, LDC )
 
570
*
 
571
*              W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
421
572
*
422
573
*              W := C1
423
574
*
424
575
               DO 160 J = 1, K
425
 
                  CALL DCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
 
576
                  CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
426
577
  160          CONTINUE
427
578
*
428
 
*              W := W * V1'
429
 
*
430
 
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', M, K,
431
 
     $                     ONE, V, LDV, WORK, LDWORK )
432
 
               IF( N.GT.K ) THEN
433
 
*
434
 
*                 W := W + C2 * V2'
435
 
*
436
 
                  CALL DGEMM( 'No transpose', 'Transpose', M, K, N-K,
437
 
     $                        ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV,
438
 
     $                        ONE, WORK, LDWORK )
 
579
*              W := W * V1**T
 
580
*
 
581
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
 
582
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
 
583
               IF( LASTV.GT.K ) THEN
 
584
*
 
585
*                 W := W + C2 * V2**T
 
586
*
 
587
                  CALL DGEMM( 'No transpose', 'Transpose',
 
588
     $                 LASTC, K, LASTV-K,
 
589
     $                 ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV,
 
590
     $                 ONE, WORK, LDWORK )
439
591
               END IF
440
592
*
441
 
*              W := W * T  or  W * T'
 
593
*              W := W * T  or  W * T**T
442
594
*
443
 
               CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
444
 
     $                     ONE, T, LDT, WORK, LDWORK )
 
595
               CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
 
596
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
445
597
*
446
598
*              C := C - W * V
447
599
*
448
 
               IF( N.GT.K ) THEN
 
600
               IF( LASTV.GT.K ) THEN
449
601
*
450
602
*                 C2 := C2 - W * V2
451
603
*
452
 
                  CALL DGEMM( 'No transpose', 'No transpose', M, N-K, K,
453
 
     $                        -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, ONE,
454
 
     $                        C( 1, K+1 ), LDC )
 
604
                  CALL DGEMM( 'No transpose', 'No transpose',
 
605
     $                 LASTC, LASTV-K, K,
 
606
     $                 -ONE, WORK, LDWORK, V( 1, K+1 ), LDV,
 
607
     $                 ONE, C( 1, K+1 ), LDC )
455
608
               END IF
456
609
*
457
610
*              W := W * V1
458
611
*
459
 
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
460
 
     $                     K, ONE, V, LDV, WORK, LDWORK )
 
612
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
 
613
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
461
614
*
462
615
*              C1 := C1 - W
463
616
*
464
617
               DO 180 J = 1, K
465
 
                  DO 170 I = 1, M
 
618
                  DO 170 I = 1, LASTC
466
619
                     C( I, J ) = C( I, J ) - WORK( I, J )
467
620
  170             CONTINUE
468
621
  180          CONTINUE
476
629
*
477
630
            IF( LSAME( SIDE, 'L' ) ) THEN
478
631
*
479
 
*              Form  H * C  or  H' * C  where  C = ( C1 )
480
 
*                                                  ( C2 )
481
 
*
482
 
*              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)
483
 
*
484
 
*              W := C2'
 
632
*              Form  H * C  or  H**T * C  where  C = ( C1 )
 
633
*                                                    ( C2 )
 
634
*
 
635
               LASTC = ILADLC( M, N, C, LDC )
 
636
*
 
637
*              W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
 
638
*
 
639
*              W := C2**T
485
640
*
486
641
               DO 190 J = 1, K
487
 
                  CALL DCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
 
642
                  CALL DCOPY( LASTC, C( M-K+J, 1 ), LDC,
 
643
     $                 WORK( 1, J ), 1 )
488
644
  190          CONTINUE
489
645
*
490
 
*              W := W * V2'
 
646
*              W := W * V2**T
491
647
*
492
 
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', N, K,
493
 
     $                     ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
 
648
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
 
649
     $              LASTC, K, ONE, V( 1, M-K+1 ), LDV,
 
650
     $              WORK, LDWORK )
494
651
               IF( M.GT.K ) THEN
495
652
*
496
 
*                 W := W + C1'*V1'
 
653
*                 W := W + C1**T * V1**T
497
654
*
498
 
                  CALL DGEMM( 'Transpose', 'Transpose', N, K, M-K, ONE,
499
 
     $                        C, LDC, V, LDV, ONE, WORK, LDWORK )
 
655
                  CALL DGEMM( 'Transpose', 'Transpose',
 
656
     $                 LASTC, K, M-K, ONE, C, LDC, V, LDV,
 
657
     $                 ONE, WORK, LDWORK )
500
658
               END IF
501
659
*
502
 
*              W := W * T'  or  W * T
503
 
*
504
 
               CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
505
 
     $                     ONE, T, LDT, WORK, LDWORK )
506
 
*
507
 
*              C := C - V' * W'
 
660
*              W := W * T**T  or  W * T
 
661
*
 
662
               CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
 
663
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
 
664
*
 
665
*              C := C - V**T * W**T
508
666
*
509
667
               IF( M.GT.K ) THEN
510
668
*
511
 
*                 C1 := C1 - V1' * W'
 
669
*                 C1 := C1 - V1**T * W**T
512
670
*
513
 
                  CALL DGEMM( 'Transpose', 'Transpose', M-K, N, K, -ONE,
514
 
     $                        V, LDV, WORK, LDWORK, ONE, C, LDC )
 
671
                  CALL DGEMM( 'Transpose', 'Transpose',
 
672
     $                 M-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
 
673
     $                 ONE, C, LDC )
515
674
               END IF
516
675
*
517
676
*              W := W * V2
518
677
*
519
 
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
520
 
     $                     K, ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
 
678
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
 
679
     $              LASTC, K, ONE, V( 1, M-K+1 ), LDV,
 
680
     $              WORK, LDWORK )
521
681
*
522
 
*              C2 := C2 - W'
 
682
*              C2 := C2 - W**T
523
683
*
524
684
               DO 210 J = 1, K
525
 
                  DO 200 I = 1, N
526
 
                     C( M-K+J, I ) = C( M-K+J, I ) - WORK( I, J )
 
685
                  DO 200 I = 1, LASTC
 
686
                     C( M-K+J, I ) = C( M-K+J, I ) - WORK(I, J)
527
687
  200             CONTINUE
528
688
  210          CONTINUE
529
689
*
530
690
            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
531
691
*
532
 
*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
533
 
*
534
 
*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)
 
692
*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
 
693
*
 
694
               LASTC = ILADLR( M, N, C, LDC )
 
695
*
 
696
*              W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
535
697
*
536
698
*              W := C2
537
699
*
538
700
               DO 220 J = 1, K
539
 
                  CALL DCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
 
701
                  CALL DCOPY( LASTC, C( 1, N-K+J ), 1,
 
702
     $                 WORK( 1, J ), 1 )
540
703
  220          CONTINUE
541
704
*
542
 
*              W := W * V2'
 
705
*              W := W * V2**T
543
706
*
544
 
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', M, K,
545
 
     $                     ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
 
707
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
 
708
     $              LASTC, K, ONE, V( 1, N-K+1 ), LDV,
 
709
     $              WORK, LDWORK )
546
710
               IF( N.GT.K ) THEN
547
711
*
548
 
*                 W := W + C1 * V1'
 
712
*                 W := W + C1 * V1**T
549
713
*
550
 
                  CALL DGEMM( 'No transpose', 'Transpose', M, K, N-K,
551
 
     $                        ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
 
714
                  CALL DGEMM( 'No transpose', 'Transpose',
 
715
     $                 LASTC, K, N-K, ONE, C, LDC, V, LDV,
 
716
     $                 ONE, WORK, LDWORK )
552
717
               END IF
553
718
*
554
 
*              W := W * T  or  W * T'
 
719
*              W := W * T  or  W * T**T
555
720
*
556
 
               CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
557
 
     $                     ONE, T, LDT, WORK, LDWORK )
 
721
               CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
 
722
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
558
723
*
559
724
*              C := C - W * V
560
725
*
562
727
*
563
728
*                 C1 := C1 - W * V1
564
729
*
565
 
                  CALL DGEMM( 'No transpose', 'No transpose', M, N-K, K,
566
 
     $                        -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
 
730
                  CALL DGEMM( 'No transpose', 'No transpose',
 
731
     $                 LASTC, N-K, K, -ONE, WORK, LDWORK, V, LDV,
 
732
     $                 ONE, C, LDC )
567
733
               END IF
568
734
*
569
735
*              W := W * V2
570
736
*
571
 
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
572
 
     $                     K, ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
 
737
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
 
738
     $              LASTC, K, ONE, V( 1, N-K+1 ), LDV,
 
739
     $              WORK, LDWORK )
573
740
*
574
741
*              C1 := C1 - W
575
742
*
576
743
               DO 240 J = 1, K
577
 
                  DO 230 I = 1, M
578
 
                     C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
 
744
                  DO 230 I = 1, LASTC
 
745
                     C( I, N-K+J ) = C( I, N-K+J ) - WORK(I, J)
579
746
  230             CONTINUE
580
747
  240          CONTINUE
581
748
*