← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/lapack/zunmlq.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/zunmlq.f
 
1
      SUBROUTINE ZUNMLQ( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
 
2
     $                   WORK, LWORK, INFO )
 
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
*     June 30, 1999
 
8
*
 
9
*     .. Scalar Arguments ..
 
10
      CHARACTER          SIDE, TRANS
 
11
      INTEGER            INFO, K, LDA, LDC, LWORK, M, N
 
12
*     ..
 
13
*     .. Array Arguments ..
 
14
      COMPLEX*16         A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
 
15
*     ..
 
16
*
 
17
*  Purpose
 
18
*  =======
 
19
*
 
20
*  ZUNMLQ overwrites the general complex M-by-N matrix C with
 
21
*
 
22
*                  SIDE = 'L'     SIDE = 'R'
 
23
*  TRANS = 'N':      Q * C          C * Q
 
24
*  TRANS = 'C':      Q**H * C       C * Q**H
 
25
*
 
26
*  where Q is a complex unitary matrix defined as the product of k
 
27
*  elementary reflectors
 
28
*
 
29
*        Q = H(k)' . . . H(2)' H(1)'
 
30
*
 
31
*  as returned by ZGELQF. Q is of order M if SIDE = 'L' and of order N
 
32
*  if SIDE = 'R'.
 
33
*
 
34
*  Arguments
 
35
*  =========
 
36
*
 
37
*  SIDE    (input) CHARACTER*1
 
38
*          = 'L': apply Q or Q**H from the Left;
 
39
*          = 'R': apply Q or Q**H from the Right.
 
40
*
 
41
*  TRANS   (input) CHARACTER*1
 
42
*          = 'N':  No transpose, apply Q;
 
43
*          = 'C':  Conjugate transpose, apply Q**H.
 
44
*
 
45
*  M       (input) INTEGER
 
46
*          The number of rows of the matrix C. M >= 0.
 
47
*
 
48
*  N       (input) INTEGER
 
49
*          The number of columns of the matrix C. N >= 0.
 
50
*
 
51
*  K       (input) INTEGER
 
52
*          The number of elementary reflectors whose product defines
 
53
*          the matrix Q.
 
54
*          If SIDE = 'L', M >= K >= 0;
 
55
*          if SIDE = 'R', N >= K >= 0.
 
56
*
 
57
*  A       (input) COMPLEX*16 array, dimension
 
58
*                               (LDA,M) if SIDE = 'L',
 
59
*                               (LDA,N) if SIDE = 'R'
 
60
*          The i-th row must contain the vector which defines the
 
61
*          elementary reflector H(i), for i = 1,2,...,k, as returned by
 
62
*          ZGELQF in the first k rows of its array argument A.
 
63
*          A is modified by the routine but restored on exit.
 
64
*
 
65
*  LDA     (input) INTEGER
 
66
*          The leading dimension of the array A. LDA >= max(1,K).
 
67
*
 
68
*  TAU     (input) COMPLEX*16 array, dimension (K)
 
69
*          TAU(i) must contain the scalar factor of the elementary
 
70
*          reflector H(i), as returned by ZGELQF.
 
71
*
 
72
*  C       (input/output) COMPLEX*16 array, dimension (LDC,N)
 
73
*          On entry, the M-by-N matrix C.
 
74
*          On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
 
75
*
 
76
*  LDC     (input) INTEGER
 
77
*          The leading dimension of the array C. LDC >= max(1,M).
 
78
*
 
79
*  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
 
80
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 
81
*
 
82
*  LWORK   (input) INTEGER
 
83
*          The dimension of the array WORK.
 
84
*          If SIDE = 'L', LWORK >= max(1,N);
 
85
*          if SIDE = 'R', LWORK >= max(1,M).
 
86
*          For optimum performance LWORK >= N*NB if SIDE 'L', and
 
87
*          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
 
88
*          blocksize.
 
89
*
 
90
*          If LWORK = -1, then a workspace query is assumed; the routine
 
91
*          only calculates the optimal size of the WORK array, returns
 
92
*          this value as the first entry of the WORK array, and no error
 
93
*          message related to LWORK is issued by XERBLA.
 
94
*
 
95
*  INFO    (output) INTEGER
 
96
*          = 0:  successful exit
 
97
*          < 0:  if INFO = -i, the i-th argument had an illegal value
 
98
*
 
99
*  =====================================================================
 
100
*
 
101
*     .. Parameters ..
 
102
      INTEGER            NBMAX, LDT
 
103
      PARAMETER          ( NBMAX = 64, LDT = NBMAX+1 )
 
104
*     ..
 
105
*     .. Local Scalars ..
 
106
      LOGICAL            LEFT, LQUERY, NOTRAN
 
107
      CHARACTER          TRANST
 
108
      INTEGER            I, I1, I2, I3, IB, IC, IINFO, IWS, JC, LDWORK,
 
109
     $                   LWKOPT, MI, NB, NBMIN, NI, NQ, NW
 
110
*     ..
 
111
*     .. Local Arrays ..
 
112
      COMPLEX*16         T( LDT, NBMAX )
 
113
*     ..
 
114
*     .. External Functions ..
 
115
      LOGICAL            LSAME
 
116
      INTEGER            ILAENV
 
117
      EXTERNAL           LSAME, ILAENV
 
118
*     ..
 
119
*     .. External Subroutines ..
 
120
      EXTERNAL           XERBLA, ZLARFB, ZLARFT, ZUNML2
 
121
*     ..
 
122
*     .. Intrinsic Functions ..
 
123
      INTRINSIC          MAX, MIN
 
124
*     ..
 
125
*     .. Executable Statements ..
 
126
*
 
127
*     Test the input arguments
 
128
*
 
129
      INFO = 0
 
130
      LEFT = LSAME( SIDE, 'L' )
 
131
      NOTRAN = LSAME( TRANS, 'N' )
 
132
      LQUERY = ( LWORK.EQ.-1 )
 
133
*
 
134
*     NQ is the order of Q and NW is the minimum dimension of WORK
 
135
*
 
136
      IF( LEFT ) THEN
 
137
         NQ = M
 
138
         NW = N
 
139
      ELSE
 
140
         NQ = N
 
141
         NW = M
 
142
      END IF
 
143
      IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
 
144
         INFO = -1
 
145
      ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
 
146
         INFO = -2
 
147
      ELSE IF( M.LT.0 ) THEN
 
148
         INFO = -3
 
149
      ELSE IF( N.LT.0 ) THEN
 
150
         INFO = -4
 
151
      ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
 
152
         INFO = -5
 
153
      ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
 
154
         INFO = -7
 
155
      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
 
156
         INFO = -10
 
157
      ELSE IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
 
158
         INFO = -12
 
159
      END IF
 
160
*
 
161
      IF( INFO.EQ.0 ) THEN
 
162
*
 
163
*        Determine the block size.  NB may be at most NBMAX, where NBMAX
 
164
*        is used to define the local array T.
 
165
*
 
166
         NB = MIN( NBMAX, ILAENV( 1, 'ZUNMLQ', SIDE // TRANS, M, N, K,
 
167
     $        -1 ) )
 
168
         LWKOPT = MAX( 1, NW )*NB
 
169
         WORK( 1 ) = LWKOPT
 
170
      END IF
 
171
*
 
172
      IF( INFO.NE.0 ) THEN
 
173
         CALL XERBLA( 'ZUNMLQ', -INFO )
 
174
         RETURN
 
175
      ELSE IF( LQUERY ) THEN
 
176
         RETURN
 
177
      END IF
 
178
*
 
179
*     Quick return if possible
 
180
*
 
181
      IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) THEN
 
182
         WORK( 1 ) = 1
 
183
         RETURN
 
184
      END IF
 
185
*
 
186
      NBMIN = 2
 
187
      LDWORK = NW
 
188
      IF( NB.GT.1 .AND. NB.LT.K ) THEN
 
189
         IWS = NW*NB
 
190
         IF( LWORK.LT.IWS ) THEN
 
191
            NB = LWORK / LDWORK
 
192
            NBMIN = MAX( 2, ILAENV( 2, 'ZUNMLQ', SIDE // TRANS, M, N, K,
 
193
     $              -1 ) )
 
194
         END IF
 
195
      ELSE
 
196
         IWS = NW
 
197
      END IF
 
198
*
 
199
      IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN
 
200
*
 
201
*        Use unblocked code
 
202
*
 
203
         CALL ZUNML2( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK,
 
204
     $                IINFO )
 
205
      ELSE
 
206
*
 
207
*        Use blocked code
 
208
*
 
209
         IF( ( LEFT .AND. NOTRAN ) .OR.
 
210
     $       ( .NOT.LEFT .AND. .NOT.NOTRAN ) ) THEN
 
211
            I1 = 1
 
212
            I2 = K
 
213
            I3 = NB
 
214
         ELSE
 
215
            I1 = ( ( K-1 ) / NB )*NB + 1
 
216
            I2 = 1
 
217
            I3 = -NB
 
218
         END IF
 
219
*
 
220
         IF( LEFT ) THEN
 
221
            NI = N
 
222
            JC = 1
 
223
         ELSE
 
224
            MI = M
 
225
            IC = 1
 
226
         END IF
 
227
*
 
228
         IF( NOTRAN ) THEN
 
229
            TRANST = 'C'
 
230
         ELSE
 
231
            TRANST = 'N'
 
232
         END IF
 
233
*
 
234
         DO 10 I = I1, I2, I3
 
235
            IB = MIN( NB, K-I+1 )
 
236
*
 
237
*           Form the triangular factor of the block reflector
 
238
*           H = H(i) H(i+1) . . . H(i+ib-1)
 
239
*
 
240
            CALL ZLARFT( 'Forward', 'Rowwise', NQ-I+1, IB, A( I, I ),
 
241
     $                   LDA, TAU( I ), T, LDT )
 
242
            IF( LEFT ) THEN
 
243
*
 
244
*              H or H' is applied to C(i:m,1:n)
 
245
*
 
246
               MI = M - I + 1
 
247
               IC = I
 
248
            ELSE
 
249
*
 
250
*              H or H' is applied to C(1:m,i:n)
 
251
*
 
252
               NI = N - I + 1
 
253
               JC = I
 
254
            END IF
 
255
*
 
256
*           Apply H or H'
 
257
*
 
258
            CALL ZLARFB( SIDE, TRANST, 'Forward', 'Rowwise', MI, NI, IB,
 
259
     $                   A( I, I ), LDA, T, LDT, C( IC, JC ), LDC, WORK,
 
260
     $                   LDWORK )
 
261
   10    CONTINUE
 
262
      END IF
 
263
      WORK( 1 ) = LWKOPT
 
264
      RETURN
 
265
*
 
266
*     End of ZUNMLQ
 
267
*
 
268
      END