← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/lapack/zlarzb.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/lapack/zlarzb.f
 
1
      SUBROUTINE ZLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
 
2
     $                   LDV, T, LDT, C, LDC, WORK, LDWORK )
 
3
*
 
4
*  -- LAPACK routine (version 3.0) --
 
5
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 
6
*     Courant Institute, Argonne National Lab, and Rice University
 
7
*     December 1, 1999
 
8
*
 
9
*     .. Scalar Arguments ..
 
10
      CHARACTER          DIRECT, SIDE, STOREV, TRANS
 
11
      INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N
 
12
*     ..
 
13
*     .. Array Arguments ..
 
14
      COMPLEX*16         C( LDC, * ), T( LDT, * ), V( LDV, * ),
 
15
     $                   WORK( LDWORK, * )
 
16
*     ..
 
17
*
 
18
*  Purpose
 
19
*  =======
 
20
*
 
21
*  ZLARZB applies a complex block reflector H or its transpose H**H
 
22
*  to a complex distributed M-by-N  C from the left or the right.
 
23
*
 
24
*  Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
 
25
*
 
26
*  Arguments
 
27
*  =========
 
28
*
 
29
*  SIDE    (input) CHARACTER*1
 
30
*          = 'L': apply H or H' from the Left
 
31
*          = 'R': apply H or H' from the Right
 
32
*
 
33
*  TRANS   (input) CHARACTER*1
 
34
*          = 'N': apply H (No transpose)
 
35
*          = 'C': apply H' (Conjugate transpose)
 
36
*
 
37
*  DIRECT  (input) CHARACTER*1
 
38
*          Indicates how H is formed from a product of elementary
 
39
*          reflectors
 
40
*          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
 
41
*          = 'B': H = H(k) . . . H(2) H(1) (Backward)
 
42
*
 
43
*  STOREV  (input) CHARACTER*1
 
44
*          Indicates how the vectors which define the elementary
 
45
*          reflectors are stored:
 
46
*          = 'C': Columnwise                        (not supported yet)
 
47
*          = 'R': Rowwise
 
48
*
 
49
*  M       (input) INTEGER
 
50
*          The number of rows of the matrix C.
 
51
*
 
52
*  N       (input) INTEGER
 
53
*          The number of columns of the matrix C.
 
54
*
 
55
*  K       (input) INTEGER
 
56
*          The order of the matrix T (= the number of elementary
 
57
*          reflectors whose product defines the block reflector).
 
58
*
 
59
*  L       (input) INTEGER
 
60
*          The number of columns of the matrix V containing the
 
61
*          meaningful part of the Householder reflectors.
 
62
*          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
 
63
*
 
64
*  V       (input) COMPLEX*16 array, dimension (LDV,NV).
 
65
*          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
 
66
*
 
67
*  LDV     (input) INTEGER
 
68
*          The leading dimension of the array V.
 
69
*          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.
 
70
*
 
71
*  T       (input) COMPLEX*16 array, dimension (LDT,K)
 
72
*          The triangular K-by-K matrix T in the representation of the
 
73
*          block reflector.
 
74
*
 
75
*  LDT     (input) INTEGER
 
76
*          The leading dimension of the array T. LDT >= K.
 
77
*
 
78
*  C       (input/output) COMPLEX*16 array, dimension (LDC,N)
 
79
*          On entry, the M-by-N matrix C.
 
80
*          On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
 
81
*
 
82
*  LDC     (input) INTEGER
 
83
*          The leading dimension of the array C. LDC >= max(1,M).
 
84
*
 
85
*  WORK    (workspace) COMPLEX*16 array, dimension (LDWORK,K)
 
86
*
 
87
*  LDWORK  (input) INTEGER
 
88
*          The leading dimension of the array WORK.
 
89
*          If SIDE = 'L', LDWORK >= max(1,N);
 
90
*          if SIDE = 'R', LDWORK >= max(1,M).
 
91
*
 
92
*  Further Details
 
93
*  ===============
 
94
*
 
95
*  Based on contributions by
 
96
*    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
 
97
*
 
98
*  =====================================================================
 
99
*
 
100
*     .. Parameters ..
 
101
      COMPLEX*16         ONE
 
102
      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
 
103
*     ..
 
104
*     .. Local Scalars ..
 
105
      CHARACTER          TRANST
 
106
      INTEGER            I, INFO, J
 
107
*     ..
 
108
*     .. External Functions ..
 
109
      LOGICAL            LSAME
 
110
      EXTERNAL           LSAME
 
111
*     ..
 
112
*     .. External Subroutines ..
 
113
      EXTERNAL           XERBLA, ZCOPY, ZGEMM, ZLACGV, ZTRMM
 
114
*     ..
 
115
*     .. Executable Statements ..
 
116
*
 
117
*     Quick return if possible
 
118
*
 
119
      IF( M.LE.0 .OR. N.LE.0 )
 
120
     $   RETURN
 
121
*
 
122
*     Check for currently supported options
 
123
*
 
124
      INFO = 0
 
125
      IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN
 
126
         INFO = -3
 
127
      ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN
 
128
         INFO = -4
 
129
      END IF
 
130
      IF( INFO.NE.0 ) THEN
 
131
         CALL XERBLA( 'ZLARZB', -INFO )
 
132
         RETURN
 
133
      END IF
 
134
*
 
135
      IF( LSAME( TRANS, 'N' ) ) THEN
 
136
         TRANST = 'C'
 
137
      ELSE
 
138
         TRANST = 'N'
 
139
      END IF
 
140
*
 
141
      IF( LSAME( SIDE, 'L' ) ) THEN
 
142
*
 
143
*        Form  H * C  or  H' * C
 
144
*
 
145
*        W( 1:n, 1:k ) = conjg( C( 1:k, 1:n )' )
 
146
*
 
147
         DO 10 J = 1, K
 
148
            CALL ZCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
 
149
   10    CONTINUE
 
150
*
 
151
*        W( 1:n, 1:k ) = W( 1:n, 1:k ) + ...
 
152
*                        conjg( C( m-l+1:m, 1:n )' ) * V( 1:k, 1:l )'
 
153
*
 
154
         IF( L.GT.0 )
 
155
     $      CALL ZGEMM( 'Transpose', 'Conjugate transpose', N, K, L,
 
156
     $                  ONE, C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK,
 
157
     $                  LDWORK )
 
158
*
 
159
*        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T'  or  W( 1:m, 1:k ) * T
 
160
*
 
161
         CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,
 
162
     $               LDT, WORK, LDWORK )
 
163
*
 
164
*        C( 1:k, 1:n ) = C( 1:k, 1:n ) - conjg( W( 1:n, 1:k )' )
 
165
*
 
166
         DO 30 J = 1, N
 
167
            DO 20 I = 1, K
 
168
               C( I, J ) = C( I, J ) - WORK( J, I )
 
169
   20       CONTINUE
 
170
   30    CONTINUE
 
171
*
 
172
*        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
 
173
*                    conjg( V( 1:k, 1:l )' ) * conjg( W( 1:n, 1:k )' )
 
174
*
 
175
         IF( L.GT.0 )
 
176
     $      CALL ZGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV,
 
177
     $                  WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC )
 
178
*
 
179
      ELSE IF( LSAME( SIDE, 'R' ) ) THEN
 
180
*
 
181
*        Form  C * H  or  C * H'
 
182
*
 
183
*        W( 1:m, 1:k ) = C( 1:m, 1:k )
 
184
*
 
185
         DO 40 J = 1, K
 
186
            CALL ZCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
 
187
   40    CONTINUE
 
188
*
 
189
*        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
 
190
*                        C( 1:m, n-l+1:n ) * conjg( V( 1:k, 1:l )' )
 
191
*
 
192
         IF( L.GT.0 )
 
193
     $      CALL ZGEMM( 'No transpose', 'Transpose', M, K, L, ONE,
 
194
     $                  C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK )
 
195
*
 
196
*        W( 1:m, 1:k ) = W( 1:m, 1:k ) * conjg( T )  or
 
197
*                        W( 1:m, 1:k ) * conjg( T' )
 
198
*
 
199
         DO 50 J = 1, K
 
200
            CALL ZLACGV( K-J+1, T( J, J ), 1 )
 
201
   50    CONTINUE
 
202
         CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T,
 
203
     $               LDT, WORK, LDWORK )
 
204
         DO 60 J = 1, K
 
205
            CALL ZLACGV( K-J+1, T( J, J ), 1 )
 
206
   60    CONTINUE
 
207
*
 
208
*        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
 
209
*
 
210
         DO 80 J = 1, K
 
211
            DO 70 I = 1, M
 
212
               C( I, J ) = C( I, J ) - WORK( I, J )
 
213
   70       CONTINUE
 
214
   80    CONTINUE
 
215
*
 
216
*        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
 
217
*                            W( 1:m, 1:k ) * conjg( V( 1:k, 1:l ) )
 
218
*
 
219
         DO 90 J = 1, L
 
220
            CALL ZLACGV( K, V( 1, J ), 1 )
 
221
   90    CONTINUE
 
222
         IF( L.GT.0 )
 
223
     $      CALL ZGEMM( 'No transpose', 'No transpose', M, L, K, -ONE,
 
224
     $                  WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC )
 
225
         DO 100 J = 1, L
 
226
            CALL ZLACGV( K, V( 1, J ), 1 )
 
227
  100    CONTINUE
 
228
*
 
229
      END IF
 
230
*
 
231
      RETURN
 
232
*
 
233
*     End of ZLARZB
 
234
*
 
235
      END