← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/lapack/zgeqp3.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/zgeqp3.f
 
1
      SUBROUTINE ZGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK,
 
2
     $                   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
      INTEGER            INFO, LDA, LWORK, M, N
 
11
*     ..
 
12
*     .. Array Arguments ..
 
13
      INTEGER            JPVT( * )
 
14
      DOUBLE PRECISION   RWORK( * )
 
15
      COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
 
16
*     ..
 
17
*
 
18
*  Purpose
 
19
*  =======
 
20
*
 
21
*  ZGEQP3 computes a QR factorization with column pivoting of a
 
22
*  matrix A:  A*P = Q*R  using Level 3 BLAS.
 
23
*
 
24
*  Arguments
 
25
*  =========
 
26
*
 
27
*  M       (input) INTEGER
 
28
*          The number of rows of the matrix A. M >= 0.
 
29
*
 
30
*  N       (input) INTEGER
 
31
*          The number of columns of the matrix A.  N >= 0.
 
32
*
 
33
*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
 
34
*          On entry, the M-by-N matrix A.
 
35
*          On exit, the upper triangle of the array contains the
 
36
*          min(M,N)-by-N upper trapezoidal matrix R; the elements below
 
37
*          the diagonal, together with the array TAU, represent the
 
38
*          unitary matrix Q as a product of min(M,N) elementary
 
39
*          reflectors.
 
40
*
 
41
*  LDA     (input) INTEGER
 
42
*          The leading dimension of the array A. LDA >= max(1,M).
 
43
*
 
44
*  JPVT    (input/output) INTEGER array, dimension (N)
 
45
*          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
 
46
*          to the front of A*P (a leading column); if JPVT(J)=0,
 
47
*          the J-th column of A is a free column.
 
48
*          On exit, if JPVT(J)=K, then the J-th column of A*P was the
 
49
*          the K-th column of A.
 
50
*
 
51
*  TAU     (output) COMPLEX*16 array, dimension (min(M,N))
 
52
*          The scalar factors of the elementary reflectors.
 
53
*
 
54
*  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
 
55
*          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
 
56
*
 
57
*  LWORK   (input) INTEGER
 
58
*          The dimension of the array WORK. LWORK >= N+1.
 
59
*          For optimal performance LWORK >= ( N+1 )*NB, where NB
 
60
*          is the optimal blocksize.
 
61
*
 
62
*          If LWORK = -1, then a workspace query is assumed; the routine
 
63
*          only calculates the optimal size of the WORK array, returns
 
64
*          this value as the first entry of the WORK array, and no error
 
65
*          message related to LWORK is issued by XERBLA.
 
66
*
 
67
*  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
 
68
*
 
69
*  INFO    (output) INTEGER
 
70
*          = 0: successful exit.
 
71
*          < 0: if INFO = -i, the i-th argument had an illegal value.
 
72
*
 
73
*  Further Details
 
74
*  ===============
 
75
*
 
76
*  The matrix Q is represented as a product of elementary reflectors
 
77
*
 
78
*     Q = H(1) H(2) . . . H(k), where k = min(m,n).
 
79
*
 
80
*  Each H(i) has the form
 
81
*
 
82
*     H(i) = I - tau * v * v'
 
83
*
 
84
*  where tau is a real/complex scalar, and v is a real/complex vector
 
85
*  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
 
86
*  A(i+1:m,i), and tau in TAU(i).
 
87
*
 
88
*  Based on contributions by
 
89
*    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
 
90
*    X. Sun, Computer Science Dept., Duke University, USA
 
91
*
 
92
*  =====================================================================
 
93
*
 
94
*     .. Parameters ..
 
95
      INTEGER            INB, INBMIN, IXOVER
 
96
      PARAMETER          ( INB = 1, INBMIN = 2, IXOVER = 3 )
 
97
*     ..
 
98
*     .. Local Scalars ..
 
99
      LOGICAL            LQUERY
 
100
      INTEGER            FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB,
 
101
     $                   NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN
 
102
*     ..
 
103
*     .. External Subroutines ..
 
104
      EXTERNAL           XERBLA, ZGEQRF, ZLAQP2, ZLAQPS, ZSWAP, ZUNMQR
 
105
*     ..
 
106
*     .. External Functions ..
 
107
      INTEGER            ILAENV
 
108
      DOUBLE PRECISION   DZNRM2
 
109
      EXTERNAL           ILAENV, DZNRM2
 
110
*     ..
 
111
*     .. Intrinsic Functions ..
 
112
      INTRINSIC          INT, MAX, MIN
 
113
*     ..
 
114
*     .. Executable Statements ..
 
115
*
 
116
      IWS = N + 1
 
117
      MINMN = MIN( M, N )
 
118
*
 
119
*     Test input arguments
 
120
*     ====================
 
121
*
 
122
      INFO = 0
 
123
      NB = ILAENV( INB, 'ZGEQRF', ' ', M, N, -1, -1 )
 
124
      LWKOPT = ( N+1 )*NB
 
125
      WORK( 1 ) = LWKOPT
 
126
      LQUERY = ( LWORK.EQ.-1 )
 
127
      IF( M.LT.0 ) THEN
 
128
         INFO = -1
 
129
      ELSE IF( N.LT.0 ) THEN
 
130
         INFO = -2
 
131
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 
132
         INFO = -4
 
133
      ELSE IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
 
134
         INFO = -8
 
135
      END IF
 
136
      IF( INFO.NE.0 ) THEN
 
137
         CALL XERBLA( 'ZGEQP3', -INFO )
 
138
         RETURN
 
139
      ELSE IF( LQUERY ) THEN
 
140
         RETURN
 
141
      END IF
 
142
*
 
143
*     Quick return if possible.
 
144
*
 
145
      IF( MINMN.EQ.0 ) THEN
 
146
         WORK( 1 ) = 1
 
147
         RETURN
 
148
      END IF
 
149
*
 
150
*     Move initial columns up front.
 
151
*
 
152
      NFXD = 1
 
153
      DO 10 J = 1, N
 
154
         IF( JPVT( J ).NE.0 ) THEN
 
155
            IF( J.NE.NFXD ) THEN
 
156
               CALL ZSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
 
157
               JPVT( J ) = JPVT( NFXD )
 
158
               JPVT( NFXD ) = J
 
159
            ELSE
 
160
               JPVT( J ) = J
 
161
            END IF
 
162
            NFXD = NFXD + 1
 
163
         ELSE
 
164
            JPVT( J ) = J
 
165
         END IF
 
166
   10 CONTINUE
 
167
      NFXD = NFXD - 1
 
168
*
 
169
*     Factorize fixed columns
 
170
*     =======================
 
171
*
 
172
*     Compute the QR factorization of fixed columns and update
 
173
*     remaining columns.
 
174
*
 
175
      IF( NFXD.GT.0 ) THEN
 
176
         NA = MIN( M, NFXD )
 
177
*CC      CALL ZGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
 
178
         CALL ZGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
 
179
         IWS = MAX( IWS, INT( WORK( 1 ) ) )
 
180
         IF( NA.LT.N ) THEN
 
181
*CC         CALL ZUNM2R( 'Left', 'Conjugate Transpose', M, N-NA,
 
182
*CC  $                   NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK,
 
183
*CC  $                   INFO )
 
184
            CALL ZUNMQR( 'Left', 'Conjugate Transpose', M, N-NA, NA, A,
 
185
     $                   LDA, TAU, A( 1, NA+1 ), LDA, WORK, LWORK,
 
186
     $                   INFO )
 
187
            IWS = MAX( IWS, INT( WORK( 1 ) ) )
 
188
         END IF
 
189
      END IF
 
190
*
 
191
*     Factorize free columns
 
192
*     ======================
 
193
*
 
194
      IF( NFXD.LT.MINMN ) THEN
 
195
*
 
196
         SM = M - NFXD
 
197
         SN = N - NFXD
 
198
         SMINMN = MINMN - NFXD
 
199
*
 
200
*        Determine the block size.
 
201
*
 
202
         NB = ILAENV( INB, 'ZGEQRF', ' ', SM, SN, -1, -1 )
 
203
         NBMIN = 2
 
204
         NX = 0
 
205
*
 
206
         IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
 
207
*
 
208
*           Determine when to cross over from blocked to unblocked code.
 
209
*
 
210
            NX = MAX( 0, ILAENV( IXOVER, 'ZGEQRF', ' ', SM, SN, -1,
 
211
     $           -1 ) )
 
212
*
 
213
*
 
214
            IF( NX.LT.SMINMN ) THEN
 
215
*
 
216
*              Determine if workspace is large enough for blocked code.
 
217
*
 
218
               MINWS = ( SN+1 )*NB
 
219
               IWS = MAX( IWS, MINWS )
 
220
               IF( LWORK.LT.MINWS ) THEN
 
221
*
 
222
*                 Not enough workspace to use optimal NB: Reduce NB and
 
223
*                 determine the minimum value of NB.
 
224
*
 
225
                  NB = LWORK / ( SN+1 )
 
226
                  NBMIN = MAX( 2, ILAENV( INBMIN, 'ZGEQRF', ' ', SM, SN,
 
227
     $                    -1, -1 ) )
 
228
*
 
229
*
 
230
               END IF
 
231
            END IF
 
232
         END IF
 
233
*
 
234
*        Initialize partial column norms. The first N elements of work
 
235
*        store the exact column norms.
 
236
*
 
237
         DO 20 J = NFXD + 1, N
 
238
            RWORK( J ) = DZNRM2( SM, A( NFXD+1, J ), 1 )
 
239
            RWORK( N+J ) = RWORK( J )
 
240
   20    CONTINUE
 
241
*
 
242
         IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
 
243
     $       ( NX.LT.SMINMN ) ) THEN
 
244
*
 
245
*           Use blocked code initially.
 
246
*
 
247
            J = NFXD + 1
 
248
*
 
249
*           Compute factorization: while loop.
 
250
*
 
251
*
 
252
            TOPBMN = MINMN - NX
 
253
   30       CONTINUE
 
254
            IF( J.LE.TOPBMN ) THEN
 
255
               JB = MIN( NB, TOPBMN-J+1 )
 
256
*
 
257
*              Factorize JB columns among columns J:N.
 
258
*
 
259
               CALL ZLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
 
260
     $                      JPVT( J ), TAU( J ), RWORK( J ),
 
261
     $                      RWORK( N+J ), WORK( 1 ), WORK( JB+1 ),
 
262
     $                      N-J+1 )
 
263
*
 
264
               J = J + FJB
 
265
               GO TO 30
 
266
            END IF
 
267
         ELSE
 
268
            J = NFXD + 1
 
269
         END IF
 
270
*
 
271
*        Use unblocked code to factor the last or only block.
 
272
*
 
273
*
 
274
         IF( J.LE.MINMN )
 
275
     $      CALL ZLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
 
276
     $                   TAU( J ), RWORK( J ), RWORK( N+J ), WORK( 1 ) )
 
277
*
 
278
      END IF
 
279
*
 
280
      WORK( 1 ) = IWS
 
281
      RETURN
 
282
*
 
283
*     End of ZGEQP3
 
284
*
 
285
      END