← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/lapack/ztzrzf.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/ztzrzf.f
 
1
      SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
 
2
*
 
3
*  -- LAPACK routine (version 3.0) --
 
4
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 
5
*     Courant Institute, Argonne National Lab, and Rice University
 
6
*     June 30, 1999
 
7
*
 
8
*     .. Scalar Arguments ..
 
9
      INTEGER            INFO, LDA, LWORK, M, N
 
10
*     ..
 
11
*     .. Array Arguments ..
 
12
      COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
 
13
*     ..
 
14
*
 
15
*  Purpose
 
16
*  =======
 
17
*
 
18
*  ZTZRZF reduces the M-by-N ( M<=N ) complex upper trapezoidal matrix A
 
19
*  to upper triangular form by means of unitary transformations.
 
20
*
 
21
*  The upper trapezoidal matrix A is factored as
 
22
*
 
23
*     A = ( R  0 ) * Z,
 
24
*
 
25
*  where Z is an N-by-N unitary matrix and R is an M-by-M upper
 
26
*  triangular matrix.
 
27
*
 
28
*  Arguments
 
29
*  =========
 
30
*
 
31
*  M       (input) INTEGER
 
32
*          The number of rows of the matrix A.  M >= 0.
 
33
*
 
34
*  N       (input) INTEGER
 
35
*          The number of columns of the matrix A.  N >= 0.
 
36
*
 
37
*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
 
38
*          On entry, the leading M-by-N upper trapezoidal part of the
 
39
*          array A must contain the matrix to be factorized.
 
40
*          On exit, the leading M-by-M upper triangular part of A
 
41
*          contains the upper triangular matrix R, and elements M+1 to
 
42
*          N of the first M rows of A, with the array TAU, represent the
 
43
*          unitary matrix Z as a product of M elementary reflectors.
 
44
*
 
45
*  LDA     (input) INTEGER
 
46
*          The leading dimension of the array A.  LDA >= max(1,M).
 
47
*
 
48
*  TAU     (output) COMPLEX*16 array, dimension (M)
 
49
*          The scalar factors of the elementary reflectors.
 
50
*
 
51
*  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
 
52
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 
53
*
 
54
*  LWORK   (input) INTEGER
 
55
*          The dimension of the array WORK.  LWORK >= max(1,M).
 
56
*          For optimum performance LWORK >= M*NB, where NB is
 
57
*          the optimal blocksize.
 
58
*
 
59
*          If LWORK = -1, then a workspace query is assumed; the routine
 
60
*          only calculates the optimal size of the WORK array, returns
 
61
*          this value as the first entry of the WORK array, and no error
 
62
*          message related to LWORK is issued by XERBLA.
 
63
*
 
64
*  INFO    (output) INTEGER
 
65
*          = 0:  successful exit
 
66
*          < 0:  if INFO = -i, the i-th argument had an illegal value
 
67
*
 
68
*  Further Details
 
69
*  ===============
 
70
*
 
71
*  Based on contributions by
 
72
*    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
 
73
*
 
74
*  The factorization is obtained by Householder's method.  The kth
 
75
*  transformation matrix, Z( k ), which is used to introduce zeros into
 
76
*  the ( m - k + 1 )th row of A, is given in the form
 
77
*
 
78
*     Z( k ) = ( I     0   ),
 
79
*              ( 0  T( k ) )
 
80
*
 
81
*  where
 
82
*
 
83
*     T( k ) = I - tau*u( k )*u( k )',   u( k ) = (   1    ),
 
84
*                                                 (   0    )
 
85
*                                                 ( z( k ) )
 
86
*
 
87
*  tau is a scalar and z( k ) is an ( n - m ) element vector.
 
88
*  tau and z( k ) are chosen to annihilate the elements of the kth row
 
89
*  of X.
 
90
*
 
91
*  The scalar tau is returned in the kth element of TAU and the vector
 
92
*  u( k ) in the kth row of A, such that the elements of z( k ) are
 
93
*  in  a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in
 
94
*  the upper triangular part of A.
 
95
*
 
96
*  Z is given by
 
97
*
 
98
*     Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).
 
99
*
 
100
*  =====================================================================
 
101
*
 
102
*     .. Parameters ..
 
103
      COMPLEX*16         ZERO
 
104
      PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 
105
*     ..
 
106
*     .. Local Scalars ..
 
107
      LOGICAL            LQUERY
 
108
      INTEGER            I, IB, IWS, KI, KK, LDWORK, LWKOPT, M1, MU, NB,
 
109
     $                   NBMIN, NX
 
110
*     ..
 
111
*     .. External Subroutines ..
 
112
      EXTERNAL           XERBLA, ZLARZB, ZLARZT, ZLATRZ
 
113
*     ..
 
114
*     .. Intrinsic Functions ..
 
115
      INTRINSIC          MAX, MIN
 
116
*     ..
 
117
*     .. External Functions ..
 
118
      INTEGER            ILAENV
 
119
      EXTERNAL           ILAENV
 
120
*     ..
 
121
*     .. Executable Statements ..
 
122
*
 
123
*     Test the input arguments
 
124
*
 
125
      INFO = 0
 
126
      LQUERY = ( LWORK.EQ.-1 )
 
127
      IF( M.LT.0 ) THEN
 
128
         INFO = -1
 
129
      ELSE IF( N.LT.M ) THEN
 
130
         INFO = -2
 
131
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 
132
         INFO = -4
 
133
      ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
 
134
         INFO = -7
 
135
      END IF
 
136
*
 
137
      IF( INFO.EQ.0 ) THEN
 
138
*
 
139
*        Determine the block size.
 
140
*
 
141
         NB = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
 
142
         LWKOPT = M*NB
 
143
         WORK( 1 ) = LWKOPT
 
144
      END IF
 
145
*
 
146
      IF( INFO.NE.0 ) THEN
 
147
         CALL XERBLA( 'ZTZRZF', -INFO )
 
148
         RETURN
 
149
      ELSE IF( LQUERY ) THEN
 
150
         RETURN
 
151
      END IF
 
152
*
 
153
*     Quick return if possible
 
154
*
 
155
      IF( M.EQ.0 ) THEN
 
156
         WORK( 1 ) = 1
 
157
         RETURN
 
158
      ELSE IF( M.EQ.N ) THEN
 
159
         DO 10 I = 1, N
 
160
            TAU( I ) = ZERO
 
161
   10    CONTINUE
 
162
         WORK( 1 ) = 1
 
163
         RETURN
 
164
      END IF
 
165
*
 
166
      NBMIN = 2
 
167
      NX = 1
 
168
      IWS = M
 
169
      IF( NB.GT.1 .AND. NB.LT.M ) THEN
 
170
*
 
171
*        Determine when to cross over from blocked to unblocked code.
 
172
*
 
173
         NX = MAX( 0, ILAENV( 3, 'ZGERQF', ' ', M, N, -1, -1 ) )
 
174
         IF( NX.LT.M ) THEN
 
175
*
 
176
*           Determine if workspace is large enough for blocked code.
 
177
*
 
178
            LDWORK = M
 
179
            IWS = LDWORK*NB
 
180
            IF( LWORK.LT.IWS ) THEN
 
181
*
 
182
*              Not enough workspace to use optimal NB:  reduce NB and
 
183
*              determine the minimum value of NB.
 
184
*
 
185
               NB = LWORK / LDWORK
 
186
               NBMIN = MAX( 2, ILAENV( 2, 'ZGERQF', ' ', M, N, -1,
 
187
     $                 -1 ) )
 
188
            END IF
 
189
         END IF
 
190
      END IF
 
191
*
 
192
      IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN
 
193
*
 
194
*        Use blocked code initially.
 
195
*        The last kk rows are handled by the block method.
 
196
*
 
197
         M1 = MIN( M+1, N )
 
198
         KI = ( ( M-NX-1 ) / NB )*NB
 
199
         KK = MIN( M, KI+NB )
 
200
*
 
201
         DO 20 I = M - KK + KI + 1, M - KK + 1, -NB
 
202
            IB = MIN( M-I+1, NB )
 
203
*
 
204
*           Compute the TZ factorization of the current block
 
205
*           A(i:i+ib-1,i:n)
 
206
*
 
207
            CALL ZLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ),
 
208
     $                   WORK )
 
209
            IF( I.GT.1 ) THEN
 
210
*
 
211
*              Form the triangular factor of the block reflector
 
212
*              H = H(i+ib-1) . . . H(i+1) H(i)
 
213
*
 
214
               CALL ZLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ),
 
215
     $                      LDA, TAU( I ), WORK, LDWORK )
 
216
*
 
217
*              Apply H to A(1:i-1,i:n) from the right
 
218
*
 
219
               CALL ZLARZB( 'Right', 'No transpose', 'Backward',
 
220
     $                      'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ),
 
221
     $                      LDA, WORK, LDWORK, A( 1, I ), LDA,
 
222
     $                      WORK( IB+1 ), LDWORK )
 
223
            END IF
 
224
   20    CONTINUE
 
225
         MU = I + NB - 1
 
226
      ELSE
 
227
         MU = M
 
228
      END IF
 
229
*
 
230
*     Use unblocked code to factor the last or only block
 
231
*
 
232
      IF( MU.GT.0 )
 
233
     $   CALL ZLATRZ( MU, N, N-M, A, LDA, TAU, WORK )
 
234
*
 
235
      WORK( 1 ) = LWKOPT
 
236
*
 
237
      RETURN
 
238
*
 
239
*     End of ZTZRZF
 
240
*
 
241
      END