← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/lapack/dormqr.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

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