LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Collaboration diagram for complex16:

Functions

subroutine zgeqpf (M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO)
 ZGEQPF More...
 
subroutine zgebak (JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
 ZGEBAK More...
 
subroutine zgebal (JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
 ZGEBAL More...
 
subroutine zgebd2 (M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO)
 ZGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm. More...
 
subroutine zgebrd (M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
 ZGEBRD More...
 
subroutine zgecon (NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
 ZGECON More...
 
subroutine zgeequ (M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
 ZGEEQU More...
 
subroutine zgeequb (M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
 ZGEEQUB More...
 
subroutine zgehd2 (N, ILO, IHI, A, LDA, TAU, WORK, INFO)
 ZGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm. More...
 
subroutine zgehrd (N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
 ZGEHRD More...
 
subroutine zgelq2 (M, N, A, LDA, TAU, WORK, INFO)
 ZGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine zgelqf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 ZGELQF More...
 
subroutine zgemqrt (SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
 ZGEMQRT More...
 
subroutine zgeql2 (M, N, A, LDA, TAU, WORK, INFO)
 ZGEQL2 computes the QL factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine zgeqlf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 ZGEQLF More...
 
subroutine zgeqp3 (M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
 ZGEQP3 More...
 
subroutine zgeqr2 (M, N, A, LDA, TAU, WORK, INFO)
 ZGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine zgeqr2p (M, N, A, LDA, TAU, WORK, INFO)
 ZGEQR2P computes the QR factorization of a general rectangular matrix with non-negative diagonal elements using an unblocked algorithm. More...
 
subroutine zgeqrf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 ZGEQRF More...
 
subroutine zgeqrfp (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 ZGEQRFP More...
 
subroutine zgeqrt (M, N, NB, A, LDA, T, LDT, WORK, INFO)
 ZGEQRT More...
 
subroutine zgeqrt2 (M, N, A, LDA, T, LDT, INFO)
 ZGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY representation of Q. More...
 
recursive subroutine zgeqrt3 (M, N, A, LDA, T, LDT, INFO)
 ZGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact WY representation of Q. More...
 
subroutine zgerfs (TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
 ZGERFS More...
 
subroutine zgerfsx (TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO)
 ZGERFSX More...
 
subroutine zgerq2 (M, N, A, LDA, TAU, WORK, INFO)
 ZGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine zgerqf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 ZGERQF More...
 
subroutine zgetf2 (M, N, A, LDA, IPIV, INFO)
 ZGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row interchanges (unblocked algorithm). More...
 
subroutine zgetrf (M, N, A, LDA, IPIV, INFO)
 ZGETRF More...
 
recursive subroutine zgetrf2 (M, N, A, LDA, IPIV, INFO)
 ZGETRF2 More...
 
subroutine zgetri (N, A, LDA, IPIV, WORK, LWORK, INFO)
 ZGETRI More...
 
subroutine zgetrs (TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 ZGETRS More...
 
subroutine zhgeqz (JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
 ZHGEQZ More...
 
subroutine zla_geamv (TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
 ZLA_GEAMV computes a matrix-vector product using a general matrix to calculate error bounds. More...
 
double precision function zla_gercond_c (TRANS, N, A, LDA, AF, LDAF, IPIV, C, CAPPLY, INFO, WORK, RWORK)
 ZLA_GERCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general matrices. More...
 
double precision function zla_gercond_x (TRANS, N, A, LDA, AF, LDAF, IPIV, X, INFO, WORK, RWORK)
 ZLA_GERCOND_X computes the infinity norm condition number of op(A)*diag(x) for general matrices. More...
 
subroutine zla_gerfsx_extended (PREC_TYPE, TRANS_TYPE, N, NRHS, A, LDA, AF, LDAF, IPIV, COLEQU, C, B, LDB, Y, LDY, BERR_OUT, N_NORMS, ERRS_N, ERRS_C, RES, AYB, DY, Y_TAIL, RCOND, ITHRESH, RTHRESH, DZ_UB, IGNORE_CWISE, INFO)
 ZLA_GERFSX_EXTENDED More...
 
double precision function zla_gerpvgrw (N, NCOLS, A, LDA, AF, LDAF)
 ZLA_GERPVGRW multiplies a square real matrix by a complex matrix. More...
 
subroutine ztgevc (SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
 ZTGEVC More...
 
subroutine ztgexc (WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
 ZTGEXC More...
 

Detailed Description

This is the group of complex16 computational functions for GE matrices

Function Documentation

subroutine zgebak ( character  JOB,
character  SIDE,
integer  N,
integer  ILO,
integer  IHI,
double precision, dimension( * )  SCALE,
integer  M,
complex*16, dimension( ldv, * )  V,
integer  LDV,
integer  INFO 
)

ZGEBAK

Download ZGEBAK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEBAK forms the right or left eigenvectors of a complex general
 matrix by backward transformation on the computed eigenvectors of the
 balanced matrix output by ZGEBAL.
Parameters
[in]JOB
          JOB is CHARACTER*1
          Specifies the type of backward transformation required:
          = 'N', do nothing, return immediately;
          = 'P', do backward transformation for permutation only;
          = 'S', do backward transformation for scaling only;
          = 'B', do backward transformations for both permutation and
                 scaling.
          JOB must be the same as the argument JOB supplied to ZGEBAL.
[in]SIDE
          SIDE is CHARACTER*1
          = 'R':  V contains right eigenvectors;
          = 'L':  V contains left eigenvectors.
[in]N
          N is INTEGER
          The number of rows of the matrix V.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER
          The integers ILO and IHI determined by ZGEBAL.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
[in]SCALE
          SCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutation and scaling factors, as returned
          by ZGEBAL.
[in]M
          M is INTEGER
          The number of columns of the matrix V.  M >= 0.
[in,out]V
          V is COMPLEX*16 array, dimension (LDV,M)
          On entry, the matrix of right or left eigenvectors to be
          transformed, as returned by ZHSEIN or ZTREVC.
          On exit, V is overwritten by the transformed eigenvectors.
[in]LDV
          LDV is INTEGER
          The leading dimension of the array V. LDV >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 133 of file zgebak.f.

133 *
134 * -- LAPACK computational routine (version 3.4.0) --
135 * -- LAPACK is a software package provided by Univ. of Tennessee, --
136 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
137 * November 2011
138 *
139 * .. Scalar Arguments ..
140  CHARACTER job, side
141  INTEGER ihi, ilo, info, ldv, m, n
142 * ..
143 * .. Array Arguments ..
144  DOUBLE PRECISION scale( * )
145  COMPLEX*16 v( ldv, * )
146 * ..
147 *
148 * =====================================================================
149 *
150 * .. Parameters ..
151  DOUBLE PRECISION one
152  parameter( one = 1.0d+0 )
153 * ..
154 * .. Local Scalars ..
155  LOGICAL leftv, rightv
156  INTEGER i, ii, k
157  DOUBLE PRECISION s
158 * ..
159 * .. External Functions ..
160  LOGICAL lsame
161  EXTERNAL lsame
162 * ..
163 * .. External Subroutines ..
164  EXTERNAL xerbla, zdscal, zswap
165 * ..
166 * .. Intrinsic Functions ..
167  INTRINSIC max, min
168 * ..
169 * .. Executable Statements ..
170 *
171 * Decode and Test the input parameters
172 *
173  rightv = lsame( side, 'R' )
174  leftv = lsame( side, 'L' )
175 *
176  info = 0
177  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
178  $ .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
179  info = -1
180  ELSE IF( .NOT.rightv .AND. .NOT.leftv ) THEN
181  info = -2
182  ELSE IF( n.LT.0 ) THEN
183  info = -3
184  ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
185  info = -4
186  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
187  info = -5
188  ELSE IF( m.LT.0 ) THEN
189  info = -7
190  ELSE IF( ldv.LT.max( 1, n ) ) THEN
191  info = -9
192  END IF
193  IF( info.NE.0 ) THEN
194  CALL xerbla( 'ZGEBAK', -info )
195  RETURN
196  END IF
197 *
198 * Quick return if possible
199 *
200  IF( n.EQ.0 )
201  $ RETURN
202  IF( m.EQ.0 )
203  $ RETURN
204  IF( lsame( job, 'N' ) )
205  $ RETURN
206 *
207  IF( ilo.EQ.ihi )
208  $ GO TO 30
209 *
210 * Backward balance
211 *
212  IF( lsame( job, 'S' ) .OR. lsame( job, 'B' ) ) THEN
213 *
214  IF( rightv ) THEN
215  DO 10 i = ilo, ihi
216  s = scale( i )
217  CALL zdscal( m, s, v( i, 1 ), ldv )
218  10 CONTINUE
219  END IF
220 *
221  IF( leftv ) THEN
222  DO 20 i = ilo, ihi
223  s = one / scale( i )
224  CALL zdscal( m, s, v( i, 1 ), ldv )
225  20 CONTINUE
226  END IF
227 *
228  END IF
229 *
230 * Backward permutation
231 *
232 * For I = ILO-1 step -1 until 1,
233 * IHI+1 step 1 until N do --
234 *
235  30 CONTINUE
236  IF( lsame( job, 'P' ) .OR. lsame( job, 'B' ) ) THEN
237  IF( rightv ) THEN
238  DO 40 ii = 1, n
239  i = ii
240  IF( i.GE.ilo .AND. i.LE.ihi )
241  $ GO TO 40
242  IF( i.LT.ilo )
243  $ i = ilo - ii
244  k = scale( i )
245  IF( k.EQ.i )
246  $ GO TO 40
247  CALL zswap( m, v( i, 1 ), ldv, v( k, 1 ), ldv )
248  40 CONTINUE
249  END IF
250 *
251  IF( leftv ) THEN
252  DO 50 ii = 1, n
253  i = ii
254  IF( i.GE.ilo .AND. i.LE.ihi )
255  $ GO TO 50
256  IF( i.LT.ilo )
257  $ i = ilo - ii
258  k = scale( i )
259  IF( k.EQ.i )
260  $ GO TO 50
261  CALL zswap( m, v( i, 1 ), ldv, v( k, 1 ), ldv )
262  50 CONTINUE
263  END IF
264  END IF
265 *
266  RETURN
267 *
268 * End of ZGEBAK
269 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgebal ( character  JOB,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer  ILO,
integer  IHI,
double precision, dimension( * )  SCALE,
integer  INFO 
)

ZGEBAL

Download ZGEBAL + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEBAL balances a general complex matrix A.  This involves, first,
 permuting A by a similarity transformation to isolate eigenvalues
 in the first 1 to ILO-1 and last IHI+1 to N elements on the
 diagonal; and second, applying a diagonal similarity transformation
 to rows and columns ILO to IHI to make the rows and columns as
 close in norm as possible.  Both steps are optional.

 Balancing may reduce the 1-norm of the matrix, and improve the
 accuracy of the computed eigenvalues and/or eigenvectors.
Parameters
[in]JOB
          JOB is CHARACTER*1
          Specifies the operations to be performed on A:
          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
                  for i = 1,...,N;
          = 'P':  permute only;
          = 'S':  scale only;
          = 'B':  both permute and scale.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the input matrix A.
          On exit,  A is overwritten by the balanced matrix.
          If JOB = 'N', A is not referenced.
          See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]ILO
 
[out]IHI
          ILO and IHI are set to INTEGER such that on exit
          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
[out]SCALE
          SCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied to
          A.  If P(j) is the index of the row and column interchanged
          with row and column j and D(j) is the scaling factor
          applied to row and column j, then
          SCALE(j) = P(j)    for j = 1,...,ILO-1
                   = D(j)    for j = ILO,...,IHI
                   = P(j)    for j = IHI+1,...,N.
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Further Details:
  The permutations consist of row and column interchanges which put
  the matrix in the form

             ( T1   X   Y  )
     P A P = (  0   B   Z  )
             (  0   0   T2 )

  where T1 and T2 are upper triangular matrices whose eigenvalues lie
  along the diagonal.  The column indices ILO and IHI mark the starting
  and ending columns of the submatrix B. Balancing consists of applying
  a diagonal similarity transformation inv(D) * B * D to make the
  1-norms of each row of B and its corresponding column nearly equal.
  The output matrix is

     ( T1     X*D          Y    )
     (  0  inv(D)*B*D  inv(D)*Z ).
     (  0      0           T2   )

  Information about the permutations P and the diagonal matrix D is
  returned in the vector SCALE.

  This subroutine is based on the EISPACK routine CBAL.

  Modified by Tzu-Yi Chen, Computer Science Division, University of
    California at Berkeley, USA

Definition at line 162 of file zgebal.f.

162 *
163 * -- LAPACK computational routine (version 3.6.0) --
164 * -- LAPACK is a software package provided by Univ. of Tennessee, --
165 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166 * November 2015
167 *
168 * .. Scalar Arguments ..
169  CHARACTER job
170  INTEGER ihi, ilo, info, lda, n
171 * ..
172 * .. Array Arguments ..
173  DOUBLE PRECISION scale( * )
174  COMPLEX*16 a( lda, * )
175 * ..
176 *
177 * =====================================================================
178 *
179 * .. Parameters ..
180  DOUBLE PRECISION zero, one
181  parameter( zero = 0.0d+0, one = 1.0d+0 )
182  DOUBLE PRECISION sclfac
183  parameter( sclfac = 2.0d+0 )
184  DOUBLE PRECISION factor
185  parameter( factor = 0.95d+0 )
186 * ..
187 * .. Local Scalars ..
188  LOGICAL noconv
189  INTEGER i, ica, iexc, ira, j, k, l, m
190  DOUBLE PRECISION c, ca, f, g, r, ra, s, sfmax1, sfmax2, sfmin1,
191  $ sfmin2
192  COMPLEX*16 cdum
193 * ..
194 * .. External Functions ..
195  LOGICAL disnan, lsame
196  INTEGER izamax
197  DOUBLE PRECISION dlamch, dznrm2
198  EXTERNAL disnan, lsame, izamax, dlamch, dznrm2
199 * ..
200 * .. External Subroutines ..
201  EXTERNAL xerbla, zdscal, zswap
202 * ..
203 * .. Intrinsic Functions ..
204  INTRINSIC abs, dble, dimag, max, min
205 *
206 * Test the input parameters
207 *
208  info = 0
209  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
210  $ .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
211  info = -1
212  ELSE IF( n.LT.0 ) THEN
213  info = -2
214  ELSE IF( lda.LT.max( 1, n ) ) THEN
215  info = -4
216  END IF
217  IF( info.NE.0 ) THEN
218  CALL xerbla( 'ZGEBAL', -info )
219  RETURN
220  END IF
221 *
222  k = 1
223  l = n
224 *
225  IF( n.EQ.0 )
226  $ GO TO 210
227 *
228  IF( lsame( job, 'N' ) ) THEN
229  DO 10 i = 1, n
230  scale( i ) = one
231  10 CONTINUE
232  GO TO 210
233  END IF
234 *
235  IF( lsame( job, 'S' ) )
236  $ GO TO 120
237 *
238 * Permutation to isolate eigenvalues if possible
239 *
240  GO TO 50
241 *
242 * Row and column exchange.
243 *
244  20 CONTINUE
245  scale( m ) = j
246  IF( j.EQ.m )
247  $ GO TO 30
248 *
249  CALL zswap( l, a( 1, j ), 1, a( 1, m ), 1 )
250  CALL zswap( n-k+1, a( j, k ), lda, a( m, k ), lda )
251 *
252  30 CONTINUE
253  GO TO ( 40, 80 )iexc
254 *
255 * Search for rows isolating an eigenvalue and push them down.
256 *
257  40 CONTINUE
258  IF( l.EQ.1 )
259  $ GO TO 210
260  l = l - 1
261 *
262  50 CONTINUE
263  DO 70 j = l, 1, -1
264 *
265  DO 60 i = 1, l
266  IF( i.EQ.j )
267  $ GO TO 60
268  IF( dble( a( j, i ) ).NE.zero .OR. dimag( a( j, i ) ).NE.
269  $ zero )GO TO 70
270  60 CONTINUE
271 *
272  m = l
273  iexc = 1
274  GO TO 20
275  70 CONTINUE
276 *
277  GO TO 90
278 *
279 * Search for columns isolating an eigenvalue and push them left.
280 *
281  80 CONTINUE
282  k = k + 1
283 *
284  90 CONTINUE
285  DO 110 j = k, l
286 *
287  DO 100 i = k, l
288  IF( i.EQ.j )
289  $ GO TO 100
290  IF( dble( a( i, j ) ).NE.zero .OR. dimag( a( i, j ) ).NE.
291  $ zero )GO TO 110
292  100 CONTINUE
293 *
294  m = k
295  iexc = 2
296  GO TO 20
297  110 CONTINUE
298 *
299  120 CONTINUE
300  DO 130 i = k, l
301  scale( i ) = one
302  130 CONTINUE
303 *
304  IF( lsame( job, 'P' ) )
305  $ GO TO 210
306 *
307 * Balance the submatrix in rows K to L.
308 *
309 * Iterative loop for norm reduction
310 *
311  sfmin1 = dlamch( 'S' ) / dlamch( 'P' )
312  sfmax1 = one / sfmin1
313  sfmin2 = sfmin1*sclfac
314  sfmax2 = one / sfmin2
315  140 CONTINUE
316  noconv = .false.
317 *
318  DO 200 i = k, l
319 *
320  c = dznrm2( l-k+1, a( k, i ), 1 )
321  r = dznrm2( l-k+1, a( i, k ), lda )
322  ica = izamax( l, a( 1, i ), 1 )
323  ca = abs( a( ica, i ) )
324  ira = izamax( n-k+1, a( i, k ), lda )
325  ra = abs( a( i, ira+k-1 ) )
326 *
327 * Guard against zero C or R due to underflow.
328 *
329  IF( c.EQ.zero .OR. r.EQ.zero )
330  $ GO TO 200
331  g = r / sclfac
332  f = one
333  s = c + r
334  160 CONTINUE
335  IF( c.GE.g .OR. max( f, c, ca ).GE.sfmax2 .OR.
336  $ min( r, g, ra ).LE.sfmin2 )GO TO 170
337  IF( disnan( c+f+ca+r+g+ra ) ) THEN
338 *
339 * Exit if NaN to avoid infinite loop
340 *
341  info = -3
342  CALL xerbla( 'ZGEBAL', -info )
343  RETURN
344  END IF
345  f = f*sclfac
346  c = c*sclfac
347  ca = ca*sclfac
348  r = r / sclfac
349  g = g / sclfac
350  ra = ra / sclfac
351  GO TO 160
352 *
353  170 CONTINUE
354  g = c / sclfac
355  180 CONTINUE
356  IF( g.LT.r .OR. max( r, ra ).GE.sfmax2 .OR.
357  $ min( f, c, g, ca ).LE.sfmin2 )GO TO 190
358  f = f / sclfac
359  c = c / sclfac
360  g = g / sclfac
361  ca = ca / sclfac
362  r = r*sclfac
363  ra = ra*sclfac
364  GO TO 180
365 *
366 * Now balance.
367 *
368  190 CONTINUE
369  IF( ( c+r ).GE.factor*s )
370  $ GO TO 200
371  IF( f.LT.one .AND. scale( i ).LT.one ) THEN
372  IF( f*scale( i ).LE.sfmin1 )
373  $ GO TO 200
374  END IF
375  IF( f.GT.one .AND. scale( i ).GT.one ) THEN
376  IF( scale( i ).GE.sfmax1 / f )
377  $ GO TO 200
378  END IF
379  g = one / f
380  scale( i ) = scale( i )*f
381  noconv = .true.
382 *
383  CALL zdscal( n-k+1, g, a( i, k ), lda )
384  CALL zdscal( l, f, a( 1, i ), 1 )
385 *
386  200 CONTINUE
387 *
388  IF( noconv )
389  $ GO TO 140
390 *
391  210 CONTINUE
392  ilo = k
393  ihi = l
394 *
395  RETURN
396 *
397 * End of ZGEBAL
398 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
double precision function dznrm2(N, X, INCX)
DZNRM2
Definition: dznrm2.f:56
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgebd2 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
complex*16, dimension( * )  TAUQ,
complex*16, dimension( * )  TAUP,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm.

Download ZGEBD2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEBD2 reduces a complex general m by n matrix A to upper or lower
 real bidiagonal form B by a unitary transformation: Q**H * A * P = B.

 If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
Parameters
[in]M
          M is INTEGER
          The number of rows in the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns in the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the m by n general matrix to be reduced.
          On exit,
          if m >= n, the diagonal and the first superdiagonal are
            overwritten with the upper bidiagonal matrix B; the
            elements below the diagonal, with the array TAUQ, represent
            the unitary matrix Q as a product of elementary
            reflectors, and the elements above the first superdiagonal,
            with the array TAUP, represent the unitary matrix P as
            a product of elementary reflectors;
          if m < n, the diagonal and the first subdiagonal are
            overwritten with the lower bidiagonal matrix B; the
            elements below the first subdiagonal, with the array TAUQ,
            represent the unitary matrix Q as a product of
            elementary reflectors, and the elements above the diagonal,
            with the array TAUP, represent the unitary matrix P as
            a product of elementary reflectors.
          See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]D
          D is DOUBLE PRECISION array, dimension (min(M,N))
          The diagonal elements of the bidiagonal matrix B:
          D(i) = A(i,i).
[out]E
          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
          The off-diagonal elements of the bidiagonal matrix B:
          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
[out]TAUQ
          TAUQ is COMPLEX*16 array dimension (min(M,N))
          The scalar factors of the elementary reflectors which
          represent the unitary matrix Q. See Further Details.
[out]TAUP
          TAUP is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors which
          represent the unitary matrix P. See Further Details.
[out]WORK
          WORK is COMPLEX*16 array, dimension (max(M,N))
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  The matrices Q and P are represented as products of elementary
  reflectors:

  If m >= n,

     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)

  Each H(i) and G(i) has the form:

     H(i) = I - tauq * v * v**H  and G(i) = I - taup * u * u**H

  where tauq and taup are complex scalars, and v and u are complex
  vectors; v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in
  A(i+1:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in
  A(i,i+2:n); tauq is stored in TAUQ(i) and taup in TAUP(i).

  If m < n,

     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)

  Each H(i) and G(i) has the form:

     H(i) = I - tauq * v * v**H  and G(i) = I - taup * u * u**H

  where tauq and taup are complex scalars, v and u are complex vectors;
  v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
  u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
  tauq is stored in TAUQ(i) and taup in TAUP(i).

  The contents of A on exit are illustrated by the following examples:

  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):

    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
    (  v1  v2  v3  v4  v5 )

  where d and e denote diagonal and off-diagonal elements of B, vi
  denotes an element of the vector defining H(i), and ui an element of
  the vector defining G(i).

Definition at line 191 of file zgebd2.f.

191 *
192 * -- LAPACK computational routine (version 3.4.2) --
193 * -- LAPACK is a software package provided by Univ. of Tennessee, --
194 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195 * September 2012
196 *
197 * .. Scalar Arguments ..
198  INTEGER info, lda, m, n
199 * ..
200 * .. Array Arguments ..
201  DOUBLE PRECISION d( * ), e( * )
202  COMPLEX*16 a( lda, * ), taup( * ), tauq( * ), work( * )
203 * ..
204 *
205 * =====================================================================
206 *
207 * .. Parameters ..
208  COMPLEX*16 zero, one
209  parameter( zero = ( 0.0d+0, 0.0d+0 ),
210  $ one = ( 1.0d+0, 0.0d+0 ) )
211 * ..
212 * .. Local Scalars ..
213  INTEGER i
214  COMPLEX*16 alpha
215 * ..
216 * .. External Subroutines ..
217  EXTERNAL xerbla, zlacgv, zlarf, zlarfg
218 * ..
219 * .. Intrinsic Functions ..
220  INTRINSIC dconjg, max, min
221 * ..
222 * .. Executable Statements ..
223 *
224 * Test the input parameters
225 *
226  info = 0
227  IF( m.LT.0 ) THEN
228  info = -1
229  ELSE IF( n.LT.0 ) THEN
230  info = -2
231  ELSE IF( lda.LT.max( 1, m ) ) THEN
232  info = -4
233  END IF
234  IF( info.LT.0 ) THEN
235  CALL xerbla( 'ZGEBD2', -info )
236  RETURN
237  END IF
238 *
239  IF( m.GE.n ) THEN
240 *
241 * Reduce to upper bidiagonal form
242 *
243  DO 10 i = 1, n
244 *
245 * Generate elementary reflector H(i) to annihilate A(i+1:m,i)
246 *
247  alpha = a( i, i )
248  CALL zlarfg( m-i+1, alpha, a( min( i+1, m ), i ), 1,
249  $ tauq( i ) )
250  d( i ) = alpha
251  a( i, i ) = one
252 *
253 * Apply H(i)**H to A(i:m,i+1:n) from the left
254 *
255  IF( i.LT.n )
256  $ CALL zlarf( 'Left', m-i+1, n-i, a( i, i ), 1,
257  $ dconjg( tauq( i ) ), a( i, i+1 ), lda, work )
258  a( i, i ) = d( i )
259 *
260  IF( i.LT.n ) THEN
261 *
262 * Generate elementary reflector G(i) to annihilate
263 * A(i,i+2:n)
264 *
265  CALL zlacgv( n-i, a( i, i+1 ), lda )
266  alpha = a( i, i+1 )
267  CALL zlarfg( n-i, alpha, a( i, min( i+2, n ) ), lda,
268  $ taup( i ) )
269  e( i ) = alpha
270  a( i, i+1 ) = one
271 *
272 * Apply G(i) to A(i+1:m,i+1:n) from the right
273 *
274  CALL zlarf( 'Right', m-i, n-i, a( i, i+1 ), lda,
275  $ taup( i ), a( i+1, i+1 ), lda, work )
276  CALL zlacgv( n-i, a( i, i+1 ), lda )
277  a( i, i+1 ) = e( i )
278  ELSE
279  taup( i ) = zero
280  END IF
281  10 CONTINUE
282  ELSE
283 *
284 * Reduce to lower bidiagonal form
285 *
286  DO 20 i = 1, m
287 *
288 * Generate elementary reflector G(i) to annihilate A(i,i+1:n)
289 *
290  CALL zlacgv( n-i+1, a( i, i ), lda )
291  alpha = a( i, i )
292  CALL zlarfg( n-i+1, alpha, a( i, min( i+1, n ) ), lda,
293  $ taup( i ) )
294  d( i ) = alpha
295  a( i, i ) = one
296 *
297 * Apply G(i) to A(i+1:m,i:n) from the right
298 *
299  IF( i.LT.m )
300  $ CALL zlarf( 'Right', m-i, n-i+1, a( i, i ), lda,
301  $ taup( i ), a( i+1, i ), lda, work )
302  CALL zlacgv( n-i+1, a( i, i ), lda )
303  a( i, i ) = d( i )
304 *
305  IF( i.LT.m ) THEN
306 *
307 * Generate elementary reflector H(i) to annihilate
308 * A(i+2:m,i)
309 *
310  alpha = a( i+1, i )
311  CALL zlarfg( m-i, alpha, a( min( i+2, m ), i ), 1,
312  $ tauq( i ) )
313  e( i ) = alpha
314  a( i+1, i ) = one
315 *
316 * Apply H(i)**H to A(i+1:m,i+1:n) from the left
317 *
318  CALL zlarf( 'Left', m-i, n-i, a( i+1, i ), 1,
319  $ dconjg( tauq( i ) ), a( i+1, i+1 ), lda,
320  $ work )
321  a( i+1, i ) = e( i )
322  ELSE
323  tauq( i ) = zero
324  END IF
325  20 CONTINUE
326  END IF
327  RETURN
328 *
329 * End of ZGEBD2
330 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgebrd ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
complex*16, dimension( * )  TAUQ,
complex*16, dimension( * )  TAUP,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGEBRD

Download ZGEBRD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEBRD reduces a general complex M-by-N matrix A to upper or lower
 bidiagonal form B by a unitary transformation: Q**H * A * P = B.

 If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
Parameters
[in]M
          M is INTEGER
          The number of rows in the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns in the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N general matrix to be reduced.
          On exit,
          if m >= n, the diagonal and the first superdiagonal are
            overwritten with the upper bidiagonal matrix B; the
            elements below the diagonal, with the array TAUQ, represent
            the unitary matrix Q as a product of elementary
            reflectors, and the elements above the first superdiagonal,
            with the array TAUP, represent the unitary matrix P as
            a product of elementary reflectors;
          if m < n, the diagonal and the first subdiagonal are
            overwritten with the lower bidiagonal matrix B; the
            elements below the first subdiagonal, with the array TAUQ,
            represent the unitary matrix Q as a product of
            elementary reflectors, and the elements above the diagonal,
            with the array TAUP, represent the unitary matrix P as
            a product of elementary reflectors.
          See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]D
          D is DOUBLE PRECISION array, dimension (min(M,N))
          The diagonal elements of the bidiagonal matrix B:
          D(i) = A(i,i).
[out]E
          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
          The off-diagonal elements of the bidiagonal matrix B:
          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
[out]TAUQ
          TAUQ is COMPLEX*16 array dimension (min(M,N))
          The scalar factors of the elementary reflectors which
          represent the unitary matrix Q. See Further Details.
[out]TAUP
          TAUP is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors which
          represent the unitary matrix P. See Further Details.
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK.  LWORK >= max(1,M,N).
          For optimum performance LWORK >= (M+N)*NB, where NB
          is the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  The matrices Q and P are represented as products of elementary
  reflectors:

  If m >= n,

     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)

  Each H(i) and G(i) has the form:

     H(i) = I - tauq * v * v**H  and G(i) = I - taup * u * u**H

  where tauq and taup are complex scalars, and v and u are complex
  vectors; v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in
  A(i+1:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in
  A(i,i+2:n); tauq is stored in TAUQ(i) and taup in TAUP(i).

  If m < n,

     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)

  Each H(i) and G(i) has the form:

     H(i) = I - tauq * v * v**H  and G(i) = I - taup * u * u**H

  where tauq and taup are complex scalars, and v and u are complex
  vectors; v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in
  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in
  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).

  The contents of A on exit are illustrated by the following examples:

  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):

    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
    (  v1  v2  v3  v4  v5 )

  where d and e denote diagonal and off-diagonal elements of B, vi
  denotes an element of the vector defining H(i), and ui an element of
  the vector defining G(i).

Definition at line 207 of file zgebrd.f.

207 *
208 * -- LAPACK computational routine (version 3.4.0) --
209 * -- LAPACK is a software package provided by Univ. of Tennessee, --
210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211 * November 2011
212 *
213 * .. Scalar Arguments ..
214  INTEGER info, lda, lwork, m, n
215 * ..
216 * .. Array Arguments ..
217  DOUBLE PRECISION d( * ), e( * )
218  COMPLEX*16 a( lda, * ), taup( * ), tauq( * ), work( * )
219 * ..
220 *
221 * =====================================================================
222 *
223 * .. Parameters ..
224  COMPLEX*16 one
225  parameter( one = ( 1.0d+0, 0.0d+0 ) )
226 * ..
227 * .. Local Scalars ..
228  LOGICAL lquery
229  INTEGER i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb,
230  $ nbmin, nx
231  DOUBLE PRECISION ws
232 * ..
233 * .. External Subroutines ..
234  EXTERNAL xerbla, zgebd2, zgemm, zlabrd
235 * ..
236 * .. Intrinsic Functions ..
237  INTRINSIC dble, max, min
238 * ..
239 * .. External Functions ..
240  INTEGER ilaenv
241  EXTERNAL ilaenv
242 * ..
243 * .. Executable Statements ..
244 *
245 * Test the input parameters
246 *
247  info = 0
248  nb = max( 1, ilaenv( 1, 'ZGEBRD', ' ', m, n, -1, -1 ) )
249  lwkopt = ( m+n )*nb
250  work( 1 ) = dble( lwkopt )
251  lquery = ( lwork.EQ.-1 )
252  IF( m.LT.0 ) THEN
253  info = -1
254  ELSE IF( n.LT.0 ) THEN
255  info = -2
256  ELSE IF( lda.LT.max( 1, m ) ) THEN
257  info = -4
258  ELSE IF( lwork.LT.max( 1, m, n ) .AND. .NOT.lquery ) THEN
259  info = -10
260  END IF
261  IF( info.LT.0 ) THEN
262  CALL xerbla( 'ZGEBRD', -info )
263  RETURN
264  ELSE IF( lquery ) THEN
265  RETURN
266  END IF
267 *
268 * Quick return if possible
269 *
270  minmn = min( m, n )
271  IF( minmn.EQ.0 ) THEN
272  work( 1 ) = 1
273  RETURN
274  END IF
275 *
276  ws = max( m, n )
277  ldwrkx = m
278  ldwrky = n
279 *
280  IF( nb.GT.1 .AND. nb.LT.minmn ) THEN
281 *
282 * Set the crossover point NX.
283 *
284  nx = max( nb, ilaenv( 3, 'ZGEBRD', ' ', m, n, -1, -1 ) )
285 *
286 * Determine when to switch from blocked to unblocked code.
287 *
288  IF( nx.LT.minmn ) THEN
289  ws = ( m+n )*nb
290  IF( lwork.LT.ws ) THEN
291 *
292 * Not enough work space for the optimal NB, consider using
293 * a smaller block size.
294 *
295  nbmin = ilaenv( 2, 'ZGEBRD', ' ', m, n, -1, -1 )
296  IF( lwork.GE.( m+n )*nbmin ) THEN
297  nb = lwork / ( m+n )
298  ELSE
299  nb = 1
300  nx = minmn
301  END IF
302  END IF
303  END IF
304  ELSE
305  nx = minmn
306  END IF
307 *
308  DO 30 i = 1, minmn - nx, nb
309 *
310 * Reduce rows and columns i:i+ib-1 to bidiagonal form and return
311 * the matrices X and Y which are needed to update the unreduced
312 * part of the matrix
313 *
314  CALL zlabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),
315  $ tauq( i ), taup( i ), work, ldwrkx,
316  $ work( ldwrkx*nb+1 ), ldwrky )
317 *
318 * Update the trailing submatrix A(i+ib:m,i+ib:n), using
319 * an update of the form A := A - V*Y**H - X*U**H
320 *
321  CALL zgemm( 'No transpose', 'Conjugate transpose', m-i-nb+1,
322  $ n-i-nb+1, nb, -one, a( i+nb, i ), lda,
323  $ work( ldwrkx*nb+nb+1 ), ldwrky, one,
324  $ a( i+nb, i+nb ), lda )
325  CALL zgemm( 'No transpose', 'No transpose', m-i-nb+1, n-i-nb+1,
326  $ nb, -one, work( nb+1 ), ldwrkx, a( i, i+nb ), lda,
327  $ one, a( i+nb, i+nb ), lda )
328 *
329 * Copy diagonal and off-diagonal elements of B back into A
330 *
331  IF( m.GE.n ) THEN
332  DO 10 j = i, i + nb - 1
333  a( j, j ) = d( j )
334  a( j, j+1 ) = e( j )
335  10 CONTINUE
336  ELSE
337  DO 20 j = i, i + nb - 1
338  a( j, j ) = d( j )
339  a( j+1, j ) = e( j )
340  20 CONTINUE
341  END IF
342  30 CONTINUE
343 *
344 * Use unblocked code to reduce the remainder of the matrix
345 *
346  CALL zgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),
347  $ tauq( i ), taup( i ), work, iinfo )
348  work( 1 ) = ws
349  RETURN
350 *
351 * End of ZGEBRD
352 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlabrd(M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY)
ZLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form.
Definition: zlabrd.f:214
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine zgebd2(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO)
ZGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm.
Definition: zgebd2.f:191
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgecon ( character  NORM,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
double precision  ANORM,
double precision  RCOND,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGECON

Download ZGECON + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGECON estimates the reciprocal of the condition number of a general
 complex matrix A, in either the 1-norm or the infinity-norm, using
 the LU factorization computed by ZGETRF.

 An estimate is obtained for norm(inv(A)), and the reciprocal of the
 condition number is computed as
    RCOND = 1 / ( norm(A) * norm(inv(A)) ).
Parameters
[in]NORM
          NORM is CHARACTER*1
          Specifies whether the 1-norm condition number or the
          infinity-norm condition number is required:
          = '1' or 'O':  1-norm;
          = 'I':         Infinity-norm.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The factors L and U from the factorization A = P*L*U
          as computed by ZGETRF.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]ANORM
          ANORM is DOUBLE PRECISION
          If NORM = '1' or 'O', the 1-norm of the original matrix A.
          If NORM = 'I', the infinity-norm of the original matrix A.
[out]RCOND
          RCOND is DOUBLE PRECISION
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(norm(A) * norm(inv(A))).
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (2*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 126 of file zgecon.f.

126 *
127 * -- LAPACK computational routine (version 3.4.0) --
128 * -- LAPACK is a software package provided by Univ. of Tennessee, --
129 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
130 * November 2011
131 *
132 * .. Scalar Arguments ..
133  CHARACTER norm
134  INTEGER info, lda, n
135  DOUBLE PRECISION anorm, rcond
136 * ..
137 * .. Array Arguments ..
138  DOUBLE PRECISION rwork( * )
139  COMPLEX*16 a( lda, * ), work( * )
140 * ..
141 *
142 * =====================================================================
143 *
144 * .. Parameters ..
145  DOUBLE PRECISION one, zero
146  parameter( one = 1.0d+0, zero = 0.0d+0 )
147 * ..
148 * .. Local Scalars ..
149  LOGICAL onenrm
150  CHARACTER normin
151  INTEGER ix, kase, kase1
152  DOUBLE PRECISION ainvnm, scale, sl, smlnum, su
153  COMPLEX*16 zdum
154 * ..
155 * .. Local Arrays ..
156  INTEGER isave( 3 )
157 * ..
158 * .. External Functions ..
159  LOGICAL lsame
160  INTEGER izamax
161  DOUBLE PRECISION dlamch
162  EXTERNAL lsame, izamax, dlamch
163 * ..
164 * .. External Subroutines ..
165  EXTERNAL xerbla, zdrscl, zlacn2, zlatrs
166 * ..
167 * .. Intrinsic Functions ..
168  INTRINSIC abs, dble, dimag, max
169 * ..
170 * .. Statement Functions ..
171  DOUBLE PRECISION cabs1
172 * ..
173 * .. Statement Function definitions ..
174  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
175 * ..
176 * .. Executable Statements ..
177 *
178 * Test the input parameters.
179 *
180  info = 0
181  onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )
182  IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN
183  info = -1
184  ELSE IF( n.LT.0 ) THEN
185  info = -2
186  ELSE IF( lda.LT.max( 1, n ) ) THEN
187  info = -4
188  ELSE IF( anorm.LT.zero ) THEN
189  info = -5
190  END IF
191  IF( info.NE.0 ) THEN
192  CALL xerbla( 'ZGECON', -info )
193  RETURN
194  END IF
195 *
196 * Quick return if possible
197 *
198  rcond = zero
199  IF( n.EQ.0 ) THEN
200  rcond = one
201  RETURN
202  ELSE IF( anorm.EQ.zero ) THEN
203  RETURN
204  END IF
205 *
206  smlnum = dlamch( 'Safe minimum' )
207 *
208 * Estimate the norm of inv(A).
209 *
210  ainvnm = zero
211  normin = 'N'
212  IF( onenrm ) THEN
213  kase1 = 1
214  ELSE
215  kase1 = 2
216  END IF
217  kase = 0
218  10 CONTINUE
219  CALL zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
220  IF( kase.NE.0 ) THEN
221  IF( kase.EQ.kase1 ) THEN
222 *
223 * Multiply by inv(L).
224 *
225  CALL zlatrs( 'Lower', 'No transpose', 'Unit', normin, n, a,
226  $ lda, work, sl, rwork, info )
227 *
228 * Multiply by inv(U).
229 *
230  CALL zlatrs( 'Upper', 'No transpose', 'Non-unit', normin, n,
231  $ a, lda, work, su, rwork( n+1 ), info )
232  ELSE
233 *
234 * Multiply by inv(U**H).
235 *
236  CALL zlatrs( 'Upper', 'Conjugate transpose', 'Non-unit',
237  $ normin, n, a, lda, work, su, rwork( n+1 ),
238  $ info )
239 *
240 * Multiply by inv(L**H).
241 *
242  CALL zlatrs( 'Lower', 'Conjugate transpose', 'Unit', normin,
243  $ n, a, lda, work, sl, rwork, info )
244  END IF
245 *
246 * Divide X by 1/(SL*SU) if doing so will not cause overflow.
247 *
248  scale = sl*su
249  normin = 'Y'
250  IF( scale.NE.one ) THEN
251  ix = izamax( n, work, 1 )
252  IF( scale.LT.cabs1( work( ix ) )*smlnum .OR. scale.EQ.zero )
253  $ GO TO 20
254  CALL zdrscl( n, scale, work, 1 )
255  END IF
256  GO TO 10
257  END IF
258 *
259 * Compute the estimate of the reciprocal condition number.
260 *
261  IF( ainvnm.NE.zero )
262  $ rcond = ( one / ainvnm ) / anorm
263 *
264  20 CONTINUE
265  RETURN
266 *
267 * End of ZGECON
268 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
ZLATRS solves a triangular system of equations with the scale factor set to prevent overflow...
Definition: zlatrs.f:241
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: zdrscl.f:86

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeequ ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  R,
double precision, dimension( * )  C,
double precision  ROWCND,
double precision  COLCND,
double precision  AMAX,
integer  INFO 
)

ZGEEQU

Download ZGEEQU + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEEQU computes row and column scalings intended to equilibrate an
 M-by-N matrix A and reduce its condition number.  R returns the row
 scale factors and C the column scale factors, chosen to try to make
 the largest element in each row and column of the matrix B with
 elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.

 R(i) and C(j) are restricted to be between SMLNUM = smallest safe
 number and BIGNUM = largest safe number.  Use of these scaling
 factors is not guaranteed to reduce the condition number of A but
 works well in practice.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The M-by-N matrix whose equilibration factors are
          to be computed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]R
          R is DOUBLE PRECISION array, dimension (M)
          If INFO = 0 or INFO > M, R contains the row scale factors
          for A.
[out]C
          C is DOUBLE PRECISION array, dimension (N)
          If INFO = 0,  C contains the column scale factors for A.
[out]ROWCND
          ROWCND is DOUBLE PRECISION
          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
          AMAX is neither too large nor too small, it is not worth
          scaling by R.
[out]COLCND
          COLCND is DOUBLE PRECISION
          If INFO = 0, COLCND contains the ratio of the smallest
          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
          worth scaling by C.
[out]AMAX
          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i,  and i is
                <= M:  the i-th row of A is exactly zero
                >  M:  the (i-M)-th column of A is exactly zero
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 142 of file zgeequ.f.

142 *
143 * -- LAPACK computational routine (version 3.4.0) --
144 * -- LAPACK is a software package provided by Univ. of Tennessee, --
145 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146 * November 2011
147 *
148 * .. Scalar Arguments ..
149  INTEGER info, lda, m, n
150  DOUBLE PRECISION amax, colcnd, rowcnd
151 * ..
152 * .. Array Arguments ..
153  DOUBLE PRECISION c( * ), r( * )
154  COMPLEX*16 a( lda, * )
155 * ..
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160  DOUBLE PRECISION one, zero
161  parameter( one = 1.0d+0, zero = 0.0d+0 )
162 * ..
163 * .. Local Scalars ..
164  INTEGER i, j
165  DOUBLE PRECISION bignum, rcmax, rcmin, smlnum
166  COMPLEX*16 zdum
167 * ..
168 * .. External Functions ..
169  DOUBLE PRECISION dlamch
170  EXTERNAL dlamch
171 * ..
172 * .. External Subroutines ..
173  EXTERNAL xerbla
174 * ..
175 * .. Intrinsic Functions ..
176  INTRINSIC abs, dble, dimag, max, min
177 * ..
178 * .. Statement Functions ..
179  DOUBLE PRECISION cabs1
180 * ..
181 * .. Statement Function definitions ..
182  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
183 * ..
184 * .. Executable Statements ..
185 *
186 * Test the input parameters.
187 *
188  info = 0
189  IF( m.LT.0 ) THEN
190  info = -1
191  ELSE IF( n.LT.0 ) THEN
192  info = -2
193  ELSE IF( lda.LT.max( 1, m ) ) THEN
194  info = -4
195  END IF
196  IF( info.NE.0 ) THEN
197  CALL xerbla( 'ZGEEQU', -info )
198  RETURN
199  END IF
200 *
201 * Quick return if possible
202 *
203  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
204  rowcnd = one
205  colcnd = one
206  amax = zero
207  RETURN
208  END IF
209 *
210 * Get machine constants.
211 *
212  smlnum = dlamch( 'S' )
213  bignum = one / smlnum
214 *
215 * Compute row scale factors.
216 *
217  DO 10 i = 1, m
218  r( i ) = zero
219  10 CONTINUE
220 *
221 * Find the maximum element in each row.
222 *
223  DO 30 j = 1, n
224  DO 20 i = 1, m
225  r( i ) = max( r( i ), cabs1( a( i, j ) ) )
226  20 CONTINUE
227  30 CONTINUE
228 *
229 * Find the maximum and minimum scale factors.
230 *
231  rcmin = bignum
232  rcmax = zero
233  DO 40 i = 1, m
234  rcmax = max( rcmax, r( i ) )
235  rcmin = min( rcmin, r( i ) )
236  40 CONTINUE
237  amax = rcmax
238 *
239  IF( rcmin.EQ.zero ) THEN
240 *
241 * Find the first zero scale factor and return an error code.
242 *
243  DO 50 i = 1, m
244  IF( r( i ).EQ.zero ) THEN
245  info = i
246  RETURN
247  END IF
248  50 CONTINUE
249  ELSE
250 *
251 * Invert the scale factors.
252 *
253  DO 60 i = 1, m
254  r( i ) = one / min( max( r( i ), smlnum ), bignum )
255  60 CONTINUE
256 *
257 * Compute ROWCND = min(R(I)) / max(R(I))
258 *
259  rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
260  END IF
261 *
262 * Compute column scale factors
263 *
264  DO 70 j = 1, n
265  c( j ) = zero
266  70 CONTINUE
267 *
268 * Find the maximum element in each column,
269 * assuming the row scaling computed above.
270 *
271  DO 90 j = 1, n
272  DO 80 i = 1, m
273  c( j ) = max( c( j ), cabs1( a( i, j ) )*r( i ) )
274  80 CONTINUE
275  90 CONTINUE
276 *
277 * Find the maximum and minimum scale factors.
278 *
279  rcmin = bignum
280  rcmax = zero
281  DO 100 j = 1, n
282  rcmin = min( rcmin, c( j ) )
283  rcmax = max( rcmax, c( j ) )
284  100 CONTINUE
285 *
286  IF( rcmin.EQ.zero ) THEN
287 *
288 * Find the first zero scale factor and return an error code.
289 *
290  DO 110 j = 1, n
291  IF( c( j ).EQ.zero ) THEN
292  info = m + j
293  RETURN
294  END IF
295  110 CONTINUE
296  ELSE
297 *
298 * Invert the scale factors.
299 *
300  DO 120 j = 1, n
301  c( j ) = one / min( max( c( j ), smlnum ), bignum )
302  120 CONTINUE
303 *
304 * Compute COLCND = min(C(J)) / max(C(J))
305 *
306  colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
307  END IF
308 *
309  RETURN
310 *
311 * End of ZGEEQU
312 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeequb ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  R,
double precision, dimension( * )  C,
double precision  ROWCND,
double precision  COLCND,
double precision  AMAX,
integer  INFO 
)

ZGEEQUB

Download ZGEEQUB + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEEQUB computes row and column scalings intended to equilibrate an
 M-by-N matrix A and reduce its condition number.  R returns the row
 scale factors and C the column scale factors, chosen to try to make
 the largest element in each row and column of the matrix B with
 elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most
 the radix.

 R(i) and C(j) are restricted to be a power of the radix between
 SMLNUM = smallest safe number and BIGNUM = largest safe number.  Use
 of these scaling factors is not guaranteed to reduce the condition
 number of A but works well in practice.

 This routine differs from ZGEEQU by restricting the scaling factors
 to a power of the radix.  Baring over- and underflow, scaling by
 these factors introduces no additional rounding errors.  However, the
 scaled entries' magnitured are no longer approximately 1 but lie
 between sqrt(radix) and 1/sqrt(radix).
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The M-by-N matrix whose equilibration factors are
          to be computed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]R
          R is DOUBLE PRECISION array, dimension (M)
          If INFO = 0 or INFO > M, R contains the row scale factors
          for A.
[out]C
          C is DOUBLE PRECISION array, dimension (N)
          If INFO = 0,  C contains the column scale factors for A.
[out]ROWCND
          ROWCND is DOUBLE PRECISION
          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
          AMAX is neither too large nor too small, it is not worth
          scaling by R.
[out]COLCND
          COLCND is DOUBLE PRECISION
          If INFO = 0, COLCND contains the ratio of the smallest
          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
          worth scaling by C.
[out]AMAX
          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i,  and i is
                <= M:  the i-th row of A is exactly zero
                >  M:  the (i-M)-th column of A is exactly zero
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 149 of file zgeequb.f.

149 *
150 * -- LAPACK computational routine (version 3.4.0) --
151 * -- LAPACK is a software package provided by Univ. of Tennessee, --
152 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
153 * November 2011
154 *
155 * .. Scalar Arguments ..
156  INTEGER info, lda, m, n
157  DOUBLE PRECISION amax, colcnd, rowcnd
158 * ..
159 * .. Array Arguments ..
160  DOUBLE PRECISION c( * ), r( * )
161  COMPLEX*16 a( lda, * )
162 * ..
163 *
164 * =====================================================================
165 *
166 * .. Parameters ..
167  DOUBLE PRECISION one, zero
168  parameter( one = 1.0d+0, zero = 0.0d+0 )
169 * ..
170 * .. Local Scalars ..
171  INTEGER i, j
172  DOUBLE PRECISION bignum, rcmax, rcmin, smlnum, radix, logrdx
173  COMPLEX*16 zdum
174 * ..
175 * .. External Functions ..
176  DOUBLE PRECISION dlamch
177  EXTERNAL dlamch
178 * ..
179 * .. External Subroutines ..
180  EXTERNAL xerbla
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC abs, max, min, log, dble, dimag
184 * ..
185 * .. Statement Functions ..
186  DOUBLE PRECISION cabs1
187 * ..
188 * .. Statement Function definitions ..
189  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
190 * ..
191 * .. Executable Statements ..
192 *
193 * Test the input parameters.
194 *
195  info = 0
196  IF( m.LT.0 ) THEN
197  info = -1
198  ELSE IF( n.LT.0 ) THEN
199  info = -2
200  ELSE IF( lda.LT.max( 1, m ) ) THEN
201  info = -4
202  END IF
203  IF( info.NE.0 ) THEN
204  CALL xerbla( 'ZGEEQUB', -info )
205  RETURN
206  END IF
207 *
208 * Quick return if possible.
209 *
210  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
211  rowcnd = one
212  colcnd = one
213  amax = zero
214  RETURN
215  END IF
216 *
217 * Get machine constants. Assume SMLNUM is a power of the radix.
218 *
219  smlnum = dlamch( 'S' )
220  bignum = one / smlnum
221  radix = dlamch( 'B' )
222  logrdx = log( radix )
223 *
224 * Compute row scale factors.
225 *
226  DO 10 i = 1, m
227  r( i ) = zero
228  10 CONTINUE
229 *
230 * Find the maximum element in each row.
231 *
232  DO 30 j = 1, n
233  DO 20 i = 1, m
234  r( i ) = max( r( i ), cabs1( a( i, j ) ) )
235  20 CONTINUE
236  30 CONTINUE
237  DO i = 1, m
238  IF( r( i ).GT.zero ) THEN
239  r( i ) = radix**int( log(r( i ) ) / logrdx )
240  END IF
241  END DO
242 *
243 * Find the maximum and minimum scale factors.
244 *
245  rcmin = bignum
246  rcmax = zero
247  DO 40 i = 1, m
248  rcmax = max( rcmax, r( i ) )
249  rcmin = min( rcmin, r( i ) )
250  40 CONTINUE
251  amax = rcmax
252 *
253  IF( rcmin.EQ.zero ) THEN
254 *
255 * Find the first zero scale factor and return an error code.
256 *
257  DO 50 i = 1, m
258  IF( r( i ).EQ.zero ) THEN
259  info = i
260  RETURN
261  END IF
262  50 CONTINUE
263  ELSE
264 *
265 * Invert the scale factors.
266 *
267  DO 60 i = 1, m
268  r( i ) = one / min( max( r( i ), smlnum ), bignum )
269  60 CONTINUE
270 *
271 * Compute ROWCND = min(R(I)) / max(R(I)).
272 *
273  rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
274  END IF
275 *
276 * Compute column scale factors.
277 *
278  DO 70 j = 1, n
279  c( j ) = zero
280  70 CONTINUE
281 *
282 * Find the maximum element in each column,
283 * assuming the row scaling computed above.
284 *
285  DO 90 j = 1, n
286  DO 80 i = 1, m
287  c( j ) = max( c( j ), cabs1( a( i, j ) )*r( i ) )
288  80 CONTINUE
289  IF( c( j ).GT.zero ) THEN
290  c( j ) = radix**int( log( c( j ) ) / logrdx )
291  END IF
292  90 CONTINUE
293 *
294 * Find the maximum and minimum scale factors.
295 *
296  rcmin = bignum
297  rcmax = zero
298  DO 100 j = 1, n
299  rcmin = min( rcmin, c( j ) )
300  rcmax = max( rcmax, c( j ) )
301  100 CONTINUE
302 *
303  IF( rcmin.EQ.zero ) THEN
304 *
305 * Find the first zero scale factor and return an error code.
306 *
307  DO 110 j = 1, n
308  IF( c( j ).EQ.zero ) THEN
309  info = m + j
310  RETURN
311  END IF
312  110 CONTINUE
313  ELSE
314 *
315 * Invert the scale factors.
316 *
317  DO 120 j = 1, n
318  c( j ) = one / min( max( c( j ), smlnum ), bignum )
319  120 CONTINUE
320 *
321 * Compute COLCND = min(C(J)) / max(C(J)).
322 *
323  colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
324  END IF
325 *
326  RETURN
327 *
328 * End of ZGEEQUB
329 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgehd2 ( integer  N,
integer  ILO,
integer  IHI,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm.

Download ZGEHD2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEHD2 reduces a complex general matrix A to upper Hessenberg form H
 by a unitary similarity transformation:  Q**H * A * Q = H .
Parameters
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER

          It is assumed that A is already upper triangular in rows
          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
          set by a previous call to ZGEBAL; otherwise they should be
          set to 1 and N respectively. See Further Details.
          1 <= ILO <= IHI <= max(1,N).
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the n by n general matrix to be reduced.
          On exit, the upper triangle and the first subdiagonal of A
          are overwritten with the upper Hessenberg matrix H, and the
          elements below the first subdiagonal, with the array TAU,
          represent the unitary matrix Q as a product of elementary
          reflectors. See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]TAU
          TAU is COMPLEX*16 array, dimension (N-1)
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  The matrix Q is represented as a product of (ihi-ilo) elementary
  reflectors

     Q = H(ilo) H(ilo+1) . . . H(ihi-1).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
  exit in A(i+2:ihi,i), and tau in TAU(i).

  The contents of A are illustrated by the following example, with
  n = 7, ilo = 2 and ihi = 6:

  on entry,                        on exit,

  ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
  (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
  (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
  (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
  (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
  (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
  (                         a )    (                          a )

  where a denotes an element of the original matrix A, h denotes a
  modified element of the upper Hessenberg matrix H, and vi denotes an
  element of the vector defining H(i).

Definition at line 151 of file zgehd2.f.

151 *
152 * -- LAPACK computational routine (version 3.4.2) --
153 * -- LAPACK is a software package provided by Univ. of Tennessee, --
154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155 * September 2012
156 *
157 * .. Scalar Arguments ..
158  INTEGER ihi, ilo, info, lda, n
159 * ..
160 * .. Array Arguments ..
161  COMPLEX*16 a( lda, * ), tau( * ), work( * )
162 * ..
163 *
164 * =====================================================================
165 *
166 * .. Parameters ..
167  COMPLEX*16 one
168  parameter( one = ( 1.0d+0, 0.0d+0 ) )
169 * ..
170 * .. Local Scalars ..
171  INTEGER i
172  COMPLEX*16 alpha
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL xerbla, zlarf, zlarfg
176 * ..
177 * .. Intrinsic Functions ..
178  INTRINSIC dconjg, max, min
179 * ..
180 * .. Executable Statements ..
181 *
182 * Test the input parameters
183 *
184  info = 0
185  IF( n.LT.0 ) THEN
186  info = -1
187  ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
188  info = -2
189  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
190  info = -3
191  ELSE IF( lda.LT.max( 1, n ) ) THEN
192  info = -5
193  END IF
194  IF( info.NE.0 ) THEN
195  CALL xerbla( 'ZGEHD2', -info )
196  RETURN
197  END IF
198 *
199  DO 10 i = ilo, ihi - 1
200 *
201 * Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
202 *
203  alpha = a( i+1, i )
204  CALL zlarfg( ihi-i, alpha, a( min( i+2, n ), i ), 1, tau( i ) )
205  a( i+1, i ) = one
206 *
207 * Apply H(i) to A(1:ihi,i+1:ihi) from the right
208 *
209  CALL zlarf( 'Right', ihi, ihi-i, a( i+1, i ), 1, tau( i ),
210  $ a( 1, i+1 ), lda, work )
211 *
212 * Apply H(i)**H to A(i+1:ihi,i+1:n) from the left
213 *
214  CALL zlarf( 'Left', ihi-i, n-i, a( i+1, i ), 1,
215  $ dconjg( tau( i ) ), a( i+1, i+1 ), lda, work )
216 *
217  a( i+1, i ) = alpha
218  10 CONTINUE
219 *
220  RETURN
221 *
222 * End of ZGEHD2
223 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgehrd ( integer  N,
integer  ILO,
integer  IHI,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGEHRD

Download ZGEHRD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEHRD reduces a complex general matrix A to upper Hessenberg form H by
 an unitary similarity transformation:  Q**H * A * Q = H .
Parameters
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER

          It is assumed that A is already upper triangular in rows
          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
          set by a previous call to ZGEBAL; otherwise they should be
          set to 1 and N respectively. See Further Details.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the N-by-N general matrix to be reduced.
          On exit, the upper triangle and the first subdiagonal of A
          are overwritten with the upper Hessenberg matrix H, and the
          elements below the first subdiagonal, with the array TAU,
          represent the unitary matrix Q as a product of elementary
          reflectors. See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]TAU
          TAU is COMPLEX*16 array, dimension (N-1)
          The scalar factors of the elementary reflectors (see Further
          Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
          zero.
[out]WORK
          WORK is COMPLEX*16 array, dimension (LWORK)
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK.  LWORK >= max(1,N).
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Further Details:
  The matrix Q is represented as a product of (ihi-ilo) elementary
  reflectors

     Q = H(ilo) H(ilo+1) . . . H(ihi-1).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
  exit in A(i+2:ihi,i), and tau in TAU(i).

  The contents of A are illustrated by the following example, with
  n = 7, ilo = 2 and ihi = 6:

  on entry,                        on exit,

  ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
  (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
  (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
  (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
  (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
  (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
  (                         a )    (                          a )

  where a denotes an element of the original matrix A, h denotes a
  modified element of the upper Hessenberg matrix H, and vi denotes an
  element of the vector defining H(i).

  This file is a slight modification of LAPACK-3.0's DGEHRD
  subroutine incorporating improvements proposed by Quintana-Orti and
  Van de Geijn (2006). (See DLAHR2.)

Definition at line 169 of file zgehrd.f.

169 *
170 * -- LAPACK computational routine (version 3.6.0) --
171 * -- LAPACK is a software package provided by Univ. of Tennessee, --
172 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
173 * November 2015
174 *
175 * .. Scalar Arguments ..
176  INTEGER ihi, ilo, info, lda, lwork, n
177 * ..
178 * .. Array Arguments ..
179  COMPLEX*16 a( lda, * ), tau( * ), work( * )
180 * ..
181 *
182 * =====================================================================
183 *
184 * .. Parameters ..
185  INTEGER nbmax, ldt, tsize
186  parameter( nbmax = 64, ldt = nbmax+1,
187  $ tsize = ldt*nbmax )
188  COMPLEX*16 zero, one
189  parameter( zero = ( 0.0d+0, 0.0d+0 ),
190  $ one = ( 1.0d+0, 0.0d+0 ) )
191 * ..
192 * .. Local Scalars ..
193  LOGICAL lquery
194  INTEGER i, ib, iinfo, iwt, j, ldwork, lwkopt, nb,
195  $ nbmin, nh, nx
196  COMPLEX*16 ei
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL zaxpy, zgehd2, zgemm, zlahr2, zlarfb, ztrmm,
200  $ xerbla
201 * ..
202 * .. Intrinsic Functions ..
203  INTRINSIC max, min
204 * ..
205 * .. External Functions ..
206  INTEGER ilaenv
207  EXTERNAL ilaenv
208 * ..
209 * .. Executable Statements ..
210 *
211 * Test the input parameters
212 *
213  info = 0
214  lquery = ( lwork.EQ.-1 )
215  IF( n.LT.0 ) THEN
216  info = -1
217  ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
218  info = -2
219  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
220  info = -3
221  ELSE IF( lda.LT.max( 1, n ) ) THEN
222  info = -5
223  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
224  info = -8
225  END IF
226 *
227  IF( info.EQ.0 ) THEN
228 *
229 * Compute the workspace requirements
230 *
231  nb = min( nbmax, ilaenv( 1, 'ZGEHRD', ' ', n, ilo, ihi, -1 ) )
232  lwkopt = n*nb + tsize
233  work( 1 ) = lwkopt
234  ENDIF
235 *
236  IF( info.NE.0 ) THEN
237  CALL xerbla( 'ZGEHRD', -info )
238  RETURN
239  ELSE IF( lquery ) THEN
240  RETURN
241  END IF
242 *
243 * Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
244 *
245  DO 10 i = 1, ilo - 1
246  tau( i ) = zero
247  10 CONTINUE
248  DO 20 i = max( 1, ihi ), n - 1
249  tau( i ) = zero
250  20 CONTINUE
251 *
252 * Quick return if possible
253 *
254  nh = ihi - ilo + 1
255  IF( nh.LE.1 ) THEN
256  work( 1 ) = 1
257  RETURN
258  END IF
259 *
260 * Determine the block size
261 *
262  nb = min( nbmax, ilaenv( 1, 'ZGEHRD', ' ', n, ilo, ihi, -1 ) )
263  nbmin = 2
264  IF( nb.GT.1 .AND. nb.LT.nh ) THEN
265 *
266 * Determine when to cross over from blocked to unblocked code
267 * (last block is always handled by unblocked code)
268 *
269  nx = max( nb, ilaenv( 3, 'ZGEHRD', ' ', n, ilo, ihi, -1 ) )
270  IF( nx.LT.nh ) THEN
271 *
272 * Determine if workspace is large enough for blocked code
273 *
274  IF( lwork.LT.n*nb+tsize ) THEN
275 *
276 * Not enough workspace to use optimal NB: determine the
277 * minimum value of NB, and reduce NB or force use of
278 * unblocked code
279 *
280  nbmin = max( 2, ilaenv( 2, 'ZGEHRD', ' ', n, ilo, ihi,
281  $ -1 ) )
282  IF( lwork.GE.(n*nbmin + tsize) ) THEN
283  nb = (lwork-tsize) / n
284  ELSE
285  nb = 1
286  END IF
287  END IF
288  END IF
289  END IF
290  ldwork = n
291 *
292  IF( nb.LT.nbmin .OR. nb.GE.nh ) THEN
293 *
294 * Use unblocked code below
295 *
296  i = ilo
297 *
298  ELSE
299 *
300 * Use blocked code
301 *
302  iwt = 1 + n*nb
303  DO 40 i = ilo, ihi - 1 - nx, nb
304  ib = min( nb, ihi-i )
305 *
306 * Reduce columns i:i+ib-1 to Hessenberg form, returning the
307 * matrices V and T of the block reflector H = I - V*T*V**H
308 * which performs the reduction, and also the matrix Y = A*V*T
309 *
310  CALL zlahr2( ihi, i, ib, a( 1, i ), lda, tau( i ),
311  $ work( iwt ), ldt, work, ldwork )
312 *
313 * Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
314 * right, computing A := A - Y * V**H. V(i+ib,ib-1) must be set
315 * to 1
316 *
317  ei = a( i+ib, i+ib-1 )
318  a( i+ib, i+ib-1 ) = one
319  CALL zgemm( 'No transpose', 'Conjugate transpose',
320  $ ihi, ihi-i-ib+1,
321  $ ib, -one, work, ldwork, a( i+ib, i ), lda, one,
322  $ a( 1, i+ib ), lda )
323  a( i+ib, i+ib-1 ) = ei
324 *
325 * Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
326 * right
327 *
328  CALL ztrmm( 'Right', 'Lower', 'Conjugate transpose',
329  $ 'Unit', i, ib-1,
330  $ one, a( i+1, i ), lda, work, ldwork )
331  DO 30 j = 0, ib-2
332  CALL zaxpy( i, -one, work( ldwork*j+1 ), 1,
333  $ a( 1, i+j+1 ), 1 )
334  30 CONTINUE
335 *
336 * Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
337 * left
338 *
339  CALL zlarfb( 'Left', 'Conjugate transpose', 'Forward',
340  $ 'Columnwise',
341  $ ihi-i, n-i-ib+1, ib, a( i+1, i ), lda,
342  $ work( iwt ), ldt, a( i+1, i+ib ), lda,
343  $ work, ldwork )
344  40 CONTINUE
345  END IF
346 *
347 * Use unblocked code to reduce the rest of the matrix
348 *
349  CALL zgehd2( n, i, ihi, a, lda, tau, work, iinfo )
350  work( 1 ) = lwkopt
351 *
352  RETURN
353 *
354 * End of ZGEHRD
355 *
subroutine zlahr2(N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)
ZLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elemen...
Definition: zlahr2.f:183
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:179
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: zlarfb.f:197
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53
subroutine zgehd2(N, ILO, IHI, A, LDA, TAU, WORK, INFO)
ZGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm...
Definition: zgehd2.f:151
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgelq2 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm.

Download ZGELQ2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGELQ2 computes an LQ factorization of a complex m by n matrix A:
 A = L * Q.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the m by n matrix A.
          On exit, the elements on and below the diagonal of the array
          contain the m by min(m,n) lower trapezoidal matrix L (L is
          lower triangular if m <= n); the elements above the diagonal,
          with the array TAU, represent the unitary matrix Q as a
          product of elementary reflectors (see Further Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (M)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(k)**H . . . H(2)**H H(1)**H, where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i-1) = 0 and v(i) = 1; conjg(v(i+1:n)) is stored on exit in
  A(i,i+1:n), and tau in TAU(i).

Definition at line 123 of file zgelq2.f.

123 *
124 * -- LAPACK computational routine (version 3.4.2) --
125 * -- LAPACK is a software package provided by Univ. of Tennessee, --
126 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
127 * September 2012
128 *
129 * .. Scalar Arguments ..
130  INTEGER info, lda, m, n
131 * ..
132 * .. Array Arguments ..
133  COMPLEX*16 a( lda, * ), tau( * ), work( * )
134 * ..
135 *
136 * =====================================================================
137 *
138 * .. Parameters ..
139  COMPLEX*16 one
140  parameter( one = ( 1.0d+0, 0.0d+0 ) )
141 * ..
142 * .. Local Scalars ..
143  INTEGER i, k
144  COMPLEX*16 alpha
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL xerbla, zlacgv, zlarf, zlarfg
148 * ..
149 * .. Intrinsic Functions ..
150  INTRINSIC max, min
151 * ..
152 * .. Executable Statements ..
153 *
154 * Test the input arguments
155 *
156  info = 0
157  IF( m.LT.0 ) THEN
158  info = -1
159  ELSE IF( n.LT.0 ) THEN
160  info = -2
161  ELSE IF( lda.LT.max( 1, m ) ) THEN
162  info = -4
163  END IF
164  IF( info.NE.0 ) THEN
165  CALL xerbla( 'ZGELQ2', -info )
166  RETURN
167  END IF
168 *
169  k = min( m, n )
170 *
171  DO 10 i = 1, k
172 *
173 * Generate elementary reflector H(i) to annihilate A(i,i+1:n)
174 *
175  CALL zlacgv( n-i+1, a( i, i ), lda )
176  alpha = a( i, i )
177  CALL zlarfg( n-i+1, alpha, a( i, min( i+1, n ) ), lda,
178  $ tau( i ) )
179  IF( i.LT.m ) THEN
180 *
181 * Apply H(i) to A(i+1:m,i:n) from the right
182 *
183  a( i, i ) = one
184  CALL zlarf( 'Right', m-i, n-i+1, a( i, i ), lda, tau( i ),
185  $ a( i+1, i ), lda, work )
186  END IF
187  a( i, i ) = alpha
188  CALL zlacgv( n-i+1, a( i, i ), lda )
189  10 CONTINUE
190  RETURN
191 *
192 * End of ZGELQ2
193 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgelqf ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGELQF

Download ZGELQF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGELQF computes an LQ factorization of a complex M-by-N matrix A:
 A = L * Q.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the elements on and below the diagonal of the array
          contain the m-by-min(m,n) lower trapezoidal matrix L (L is
          lower triangular if m <= n); the elements above the diagonal,
          with the array TAU, represent the unitary matrix Q as a
          product of elementary reflectors (see Further Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,M).
          For optimum performance LWORK >= M*NB, where NB is the
          optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(k)**H . . . H(2)**H H(1)**H, where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i-1) = 0 and v(i) = 1; conjg(v(i+1:n)) is stored on exit in
  A(i,i+1:n), and tau in TAU(i).

Definition at line 137 of file zgelqf.f.

137 *
138 * -- LAPACK computational routine (version 3.4.0) --
139 * -- LAPACK is a software package provided by Univ. of Tennessee, --
140 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
141 * November 2011
142 *
143 * .. Scalar Arguments ..
144  INTEGER info, lda, lwork, m, n
145 * ..
146 * .. Array Arguments ..
147  COMPLEX*16 a( lda, * ), tau( * ), work( * )
148 * ..
149 *
150 * =====================================================================
151 *
152 * .. Local Scalars ..
153  LOGICAL lquery
154  INTEGER i, ib, iinfo, iws, k, ldwork, lwkopt, nb,
155  $ nbmin, nx
156 * ..
157 * .. External Subroutines ..
158  EXTERNAL xerbla, zgelq2, zlarfb, zlarft
159 * ..
160 * .. Intrinsic Functions ..
161  INTRINSIC max, min
162 * ..
163 * .. External Functions ..
164  INTEGER ilaenv
165  EXTERNAL ilaenv
166 * ..
167 * .. Executable Statements ..
168 *
169 * Test the input arguments
170 *
171  info = 0
172  nb = ilaenv( 1, 'ZGELQF', ' ', m, n, -1, -1 )
173  lwkopt = m*nb
174  work( 1 ) = lwkopt
175  lquery = ( lwork.EQ.-1 )
176  IF( m.LT.0 ) THEN
177  info = -1
178  ELSE IF( n.LT.0 ) THEN
179  info = -2
180  ELSE IF( lda.LT.max( 1, m ) ) THEN
181  info = -4
182  ELSE IF( lwork.LT.max( 1, m ) .AND. .NOT.lquery ) THEN
183  info = -7
184  END IF
185  IF( info.NE.0 ) THEN
186  CALL xerbla( 'ZGELQF', -info )
187  RETURN
188  ELSE IF( lquery ) THEN
189  RETURN
190  END IF
191 *
192 * Quick return if possible
193 *
194  k = min( m, n )
195  IF( k.EQ.0 ) THEN
196  work( 1 ) = 1
197  RETURN
198  END IF
199 *
200  nbmin = 2
201  nx = 0
202  iws = m
203  IF( nb.GT.1 .AND. nb.LT.k ) THEN
204 *
205 * Determine when to cross over from blocked to unblocked code.
206 *
207  nx = max( 0, ilaenv( 3, 'ZGELQF', ' ', m, n, -1, -1 ) )
208  IF( nx.LT.k ) THEN
209 *
210 * Determine if workspace is large enough for blocked code.
211 *
212  ldwork = m
213  iws = ldwork*nb
214  IF( lwork.LT.iws ) THEN
215 *
216 * Not enough workspace to use optimal NB: reduce NB and
217 * determine the minimum value of NB.
218 *
219  nb = lwork / ldwork
220  nbmin = max( 2, ilaenv( 2, 'ZGELQF', ' ', m, n, -1,
221  $ -1 ) )
222  END IF
223  END IF
224  END IF
225 *
226  IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
227 *
228 * Use blocked code initially
229 *
230  DO 10 i = 1, k - nx, nb
231  ib = min( k-i+1, nb )
232 *
233 * Compute the LQ factorization of the current block
234 * A(i:i+ib-1,i:n)
235 *
236  CALL zgelq2( ib, n-i+1, a( i, i ), lda, tau( i ), work,
237  $ iinfo )
238  IF( i+ib.LE.m ) THEN
239 *
240 * Form the triangular factor of the block reflector
241 * H = H(i) H(i+1) . . . H(i+ib-1)
242 *
243  CALL zlarft( 'Forward', 'Rowwise', n-i+1, ib, a( i, i ),
244  $ lda, tau( i ), work, ldwork )
245 *
246 * Apply H to A(i+ib:m,i:n) from the right
247 *
248  CALL zlarfb( 'Right', 'No transpose', 'Forward',
249  $ 'Rowwise', m-i-ib+1, n-i+1, ib, a( i, i ),
250  $ lda, work, ldwork, a( i+ib, i ), lda,
251  $ work( ib+1 ), ldwork )
252  END IF
253  10 CONTINUE
254  ELSE
255  i = 1
256  END IF
257 *
258 * Use unblocked code to factor the last or only block.
259 *
260  IF( i.LE.k )
261  $ CALL zgelq2( m-i+1, n-i+1, a( i, i ), lda, tau( i ), work,
262  $ iinfo )
263 *
264  work( 1 ) = iws
265  RETURN
266 *
267 * End of ZGELQF
268 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: zlarfb.f:197
subroutine zlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
ZLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: zlarft.f:165
subroutine zgelq2(M, N, A, LDA, TAU, WORK, INFO)
ZGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: zgelq2.f:123
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgemqrt ( character  SIDE,
character  TRANS,
integer  M,
integer  N,
integer  K,
integer  NB,
complex*16, dimension( ldv, * )  V,
integer  LDV,
complex*16, dimension( ldt, * )  T,
integer  LDT,
complex*16, dimension( ldc, * )  C,
integer  LDC,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZGEMQRT

Download ZGEMQRT + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEMQRT overwrites the general complex M-by-N matrix C with

                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q C            C Q
 TRANS = 'C':    Q**H C            C Q**H

 where Q is a complex orthogonal matrix defined as the product of K
 elementary reflectors:

       Q = H(1) H(2) . . . H(K) = I - V T V**H

 generated using the compact WY representation as returned by ZGEQRT. 

 Q is of order M if SIDE = 'L' and of order N  if SIDE = 'R'.
Parameters
[in]SIDE
          SIDE is CHARACTER*1
          = 'L': apply Q or Q**H from the Left;
          = 'R': apply Q or Q**H from the Right.
[in]TRANS
          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'C':  Transpose, apply Q**H.
[in]M
          M is INTEGER
          The number of rows of the matrix C. M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix C. N >= 0.
[in]K
          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.
[in]NB
          NB is INTEGER
          The block size used for the storage of T.  K >= NB >= 1.
          This must be the same value of NB used to generate T
          in CGEQRT.
[in]V
          V is COMPLEX*16 array, dimension (LDV,K)
          The i-th column must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          CGEQRT in the first K columns of its array argument A.
[in]LDV
          LDV is INTEGER
          The leading dimension of the array V.
          If SIDE = 'L', LDA >= max(1,M);
          if SIDE = 'R', LDA >= max(1,N).
[in]T
          T is COMPLEX*16 array, dimension (LDT,K)
          The upper triangular factors of the block reflectors
          as returned by CGEQRT, stored as a NB-by-N matrix.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T.  LDT >= NB.
[in,out]C
          C is COMPLEX*16 array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q C, Q**H C, C Q**H or C Q.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).
[out]WORK
          WORK is COMPLEX*16 array. The dimension of WORK is
           N*NB if SIDE = 'L', or  M*NB if SIDE = 'R'.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013

Definition at line 170 of file zgemqrt.f.

170 *
171 * -- LAPACK computational routine (version 3.5.0) --
172 * -- LAPACK is a software package provided by Univ. of Tennessee, --
173 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
174 * November 2013
175 *
176 * .. Scalar Arguments ..
177  CHARACTER side, trans
178  INTEGER info, k, ldv, ldc, m, n, nb, ldt
179 * ..
180 * .. Array Arguments ..
181  COMPLEX*16 v( ldv, * ), c( ldc, * ), t( ldt, * ), work( * )
182 * ..
183 *
184 * =====================================================================
185 *
186 * ..
187 * .. Local Scalars ..
188  LOGICAL left, right, tran, notran
189  INTEGER i, ib, ldwork, kf, q
190 * ..
191 * .. External Functions ..
192  LOGICAL lsame
193  EXTERNAL lsame
194 * ..
195 * .. External Subroutines ..
196  EXTERNAL xerbla, zlarfb
197 * ..
198 * .. Intrinsic Functions ..
199  INTRINSIC max, min
200 * ..
201 * .. Executable Statements ..
202 *
203 * .. Test the input arguments ..
204 *
205  info = 0
206  left = lsame( side, 'L' )
207  right = lsame( side, 'R' )
208  tran = lsame( trans, 'C' )
209  notran = lsame( trans, 'N' )
210 *
211  IF( left ) THEN
212  ldwork = max( 1, n )
213  q = m
214  ELSE IF ( right ) THEN
215  ldwork = max( 1, m )
216  q = n
217  END IF
218  IF( .NOT.left .AND. .NOT.right ) THEN
219  info = -1
220  ELSE IF( .NOT.tran .AND. .NOT.notran ) THEN
221  info = -2
222  ELSE IF( m.LT.0 ) THEN
223  info = -3
224  ELSE IF( n.LT.0 ) THEN
225  info = -4
226  ELSE IF( k.LT.0 .OR. k.GT.q ) THEN
227  info = -5
228  ELSE IF( nb.LT.1 .OR. (nb.GT.k .AND. k.GT.0)) THEN
229  info = -6
230  ELSE IF( ldv.LT.max( 1, q ) ) THEN
231  info = -8
232  ELSE IF( ldt.LT.nb ) THEN
233  info = -10
234  ELSE IF( ldc.LT.max( 1, m ) ) THEN
235  info = -12
236  END IF
237 *
238  IF( info.NE.0 ) THEN
239  CALL xerbla( 'ZGEMQRT', -info )
240  RETURN
241  END IF
242 *
243 * .. Quick return if possible ..
244 *
245  IF( m.EQ.0 .OR. n.EQ.0 .OR. k.EQ.0 ) RETURN
246 *
247  IF( left .AND. tran ) THEN
248 *
249  DO i = 1, k, nb
250  ib = min( nb, k-i+1 )
251  CALL zlarfb( 'L', 'C', 'F', 'C', m-i+1, n, ib,
252  $ v( i, i ), ldv, t( 1, i ), ldt,
253  $ c( i, 1 ), ldc, work, ldwork )
254  END DO
255 *
256  ELSE IF( right .AND. notran ) THEN
257 *
258  DO i = 1, k, nb
259  ib = min( nb, k-i+1 )
260  CALL zlarfb( 'R', 'N', 'F', 'C', m, n-i+1, ib,
261  $ v( i, i ), ldv, t( 1, i ), ldt,
262  $ c( 1, i ), ldc, work, ldwork )
263  END DO
264 *
265  ELSE IF( left .AND. notran ) THEN
266 *
267  kf = ((k-1)/nb)*nb+1
268  DO i = kf, 1, -nb
269  ib = min( nb, k-i+1 )
270  CALL zlarfb( 'L', 'N', 'F', 'C', m-i+1, n, ib,
271  $ v( i, i ), ldv, t( 1, i ), ldt,
272  $ c( i, 1 ), ldc, work, ldwork )
273  END DO
274 *
275  ELSE IF( right .AND. tran ) THEN
276 *
277  kf = ((k-1)/nb)*nb+1
278  DO i = kf, 1, -nb
279  ib = min( nb, k-i+1 )
280  CALL zlarfb( 'R', 'C', 'F', 'C', m, n-i+1, ib,
281  $ v( i, i ), ldv, t( 1, i ), ldt,
282  $ c( 1, i ), ldc, work, ldwork )
283  END DO
284 *
285  END IF
286 *
287  RETURN
288 *
289 * End of ZGEMQRT
290 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: zlarfb.f:197
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeql2 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZGEQL2 computes the QL factorization of a general rectangular matrix using an unblocked algorithm.

Download ZGEQL2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQL2 computes a QL factorization of a complex m by n matrix A:
 A = Q * L.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the m by n matrix A.
          On exit, if m >= n, the lower triangle of the subarray
          A(m-n+1:m,1:n) contains the n by n lower triangular matrix L;
          if m <= n, the elements on and below the (n-m)-th
          superdiagonal contain the m by n lower trapezoidal matrix L;
          the remaining elements, with the array TAU, represent the
          unitary matrix Q as a product of elementary reflectors
          (see Further Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(k) . . . H(2) H(1), where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(m-k+i+1:m) = 0 and v(m-k+i) = 1; v(1:m-k+i-1) is stored on exit in
  A(1:m-k+i-1,n-k+i), and tau in TAU(i).

Definition at line 125 of file zgeql2.f.

125 *
126 * -- LAPACK computational routine (version 3.4.2) --
127 * -- LAPACK is a software package provided by Univ. of Tennessee, --
128 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
129 * September 2012
130 *
131 * .. Scalar Arguments ..
132  INTEGER info, lda, m, n
133 * ..
134 * .. Array Arguments ..
135  COMPLEX*16 a( lda, * ), tau( * ), work( * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  COMPLEX*16 one
142  parameter( one = ( 1.0d+0, 0.0d+0 ) )
143 * ..
144 * .. Local Scalars ..
145  INTEGER i, k
146  COMPLEX*16 alpha
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL xerbla, zlarf, zlarfg
150 * ..
151 * .. Intrinsic Functions ..
152  INTRINSIC dconjg, max, min
153 * ..
154 * .. Executable Statements ..
155 *
156 * Test the input arguments
157 *
158  info = 0
159  IF( m.LT.0 ) THEN
160  info = -1
161  ELSE IF( n.LT.0 ) THEN
162  info = -2
163  ELSE IF( lda.LT.max( 1, m ) ) THEN
164  info = -4
165  END IF
166  IF( info.NE.0 ) THEN
167  CALL xerbla( 'ZGEQL2', -info )
168  RETURN
169  END IF
170 *
171  k = min( m, n )
172 *
173  DO 10 i = k, 1, -1
174 *
175 * Generate elementary reflector H(i) to annihilate
176 * A(1:m-k+i-1,n-k+i)
177 *
178  alpha = a( m-k+i, n-k+i )
179  CALL zlarfg( m-k+i, alpha, a( 1, n-k+i ), 1, tau( i ) )
180 *
181 * Apply H(i)**H to A(1:m-k+i,1:n-k+i-1) from the left
182 *
183  a( m-k+i, n-k+i ) = one
184  CALL zlarf( 'Left', m-k+i, n-k+i-1, a( 1, n-k+i ), 1,
185  $ dconjg( tau( i ) ), a, lda, work )
186  a( m-k+i, n-k+i ) = alpha
187  10 CONTINUE
188  RETURN
189 *
190 * End of ZGEQL2
191 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeqlf ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGEQLF

Download ZGEQLF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQLF computes a QL factorization of a complex M-by-N matrix A:
 A = Q * L.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit,
          if m >= n, the lower triangle of the subarray
          A(m-n+1:m,1:n) contains the N-by-N lower triangular matrix L;
          if m <= n, the elements on and below the (n-m)-th
          superdiagonal contain the M-by-N lower trapezoidal matrix L;
          the remaining elements, with the array TAU, represent the
          unitary matrix Q as a product of elementary reflectors
          (see Further Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,N).
          For optimum performance LWORK >= N*NB, where NB is
          the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(k) . . . H(2) H(1), where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(m-k+i+1:m) = 0 and v(m-k+i) = 1; v(1:m-k+i-1) is stored on exit in
  A(1:m-k+i-1,n-k+i), and tau in TAU(i).

Definition at line 140 of file zgeqlf.f.

140 *
141 * -- LAPACK computational routine (version 3.4.0) --
142 * -- LAPACK is a software package provided by Univ. of Tennessee, --
143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144 * November 2011
145 *
146 * .. Scalar Arguments ..
147  INTEGER info, lda, lwork, m, n
148 * ..
149 * .. Array Arguments ..
150  COMPLEX*16 a( lda, * ), tau( * ), work( * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Local Scalars ..
156  LOGICAL lquery
157  INTEGER i, ib, iinfo, iws, k, ki, kk, ldwork, lwkopt,
158  $ mu, nb, nbmin, nu, nx
159 * ..
160 * .. External Subroutines ..
161  EXTERNAL xerbla, zgeql2, zlarfb, zlarft
162 * ..
163 * .. Intrinsic Functions ..
164  INTRINSIC max, min
165 * ..
166 * .. External Functions ..
167  INTEGER ilaenv
168  EXTERNAL ilaenv
169 * ..
170 * .. Executable Statements ..
171 *
172 * Test the input arguments
173 *
174  info = 0
175  lquery = ( lwork.EQ.-1 )
176  IF( m.LT.0 ) THEN
177  info = -1
178  ELSE IF( n.LT.0 ) THEN
179  info = -2
180  ELSE IF( lda.LT.max( 1, m ) ) THEN
181  info = -4
182  END IF
183 *
184  IF( info.EQ.0 ) THEN
185  k = min( m, n )
186  IF( k.EQ.0 ) THEN
187  lwkopt = 1
188  ELSE
189  nb = ilaenv( 1, 'ZGEQLF', ' ', m, n, -1, -1 )
190  lwkopt = n*nb
191  END IF
192  work( 1 ) = lwkopt
193 *
194  IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
195  info = -7
196  END IF
197  END IF
198 *
199  IF( info.NE.0 ) THEN
200  CALL xerbla( 'ZGEQLF', -info )
201  RETURN
202  ELSE IF( lquery ) THEN
203  RETURN
204  END IF
205 *
206 * Quick return if possible
207 *
208  IF( k.EQ.0 ) THEN
209  RETURN
210  END IF
211 *
212  nbmin = 2
213  nx = 1
214  iws = n
215  IF( nb.GT.1 .AND. nb.LT.k ) THEN
216 *
217 * Determine when to cross over from blocked to unblocked code.
218 *
219  nx = max( 0, ilaenv( 3, 'ZGEQLF', ' ', m, n, -1, -1 ) )
220  IF( nx.LT.k ) THEN
221 *
222 * Determine if workspace is large enough for blocked code.
223 *
224  ldwork = n
225  iws = ldwork*nb
226  IF( lwork.LT.iws ) THEN
227 *
228 * Not enough workspace to use optimal NB: reduce NB and
229 * determine the minimum value of NB.
230 *
231  nb = lwork / ldwork
232  nbmin = max( 2, ilaenv( 2, 'ZGEQLF', ' ', m, n, -1,
233  $ -1 ) )
234  END IF
235  END IF
236  END IF
237 *
238  IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
239 *
240 * Use blocked code initially.
241 * The last kk columns are handled by the block method.
242 *
243  ki = ( ( k-nx-1 ) / nb )*nb
244  kk = min( k, ki+nb )
245 *
246  DO 10 i = k - kk + ki + 1, k - kk + 1, -nb
247  ib = min( k-i+1, nb )
248 *
249 * Compute the QL factorization of the current block
250 * A(1:m-k+i+ib-1,n-k+i:n-k+i+ib-1)
251 *
252  CALL zgeql2( m-k+i+ib-1, ib, a( 1, n-k+i ), lda, tau( i ),
253  $ work, iinfo )
254  IF( n-k+i.GT.1 ) THEN
255 *
256 * Form the triangular factor of the block reflector
257 * H = H(i+ib-1) . . . H(i+1) H(i)
258 *
259  CALL zlarft( 'Backward', 'Columnwise', m-k+i+ib-1, ib,
260  $ a( 1, n-k+i ), lda, tau( i ), work, ldwork )
261 *
262 * Apply H**H to A(1:m-k+i+ib-1,1:n-k+i-1) from the left
263 *
264  CALL zlarfb( 'Left', 'Conjugate transpose', 'Backward',
265  $ 'Columnwise', m-k+i+ib-1, n-k+i-1, ib,
266  $ a( 1, n-k+i ), lda, work, ldwork, a, lda,
267  $ work( ib+1 ), ldwork )
268  END IF
269  10 CONTINUE
270  mu = m - k + i + nb - 1
271  nu = n - k + i + nb - 1
272  ELSE
273  mu = m
274  nu = n
275  END IF
276 *
277 * Use unblocked code to factor the last or only block
278 *
279  IF( mu.GT.0 .AND. nu.GT.0 )
280  $ CALL zgeql2( mu, nu, a, lda, tau, work, iinfo )
281 *
282  work( 1 ) = iws
283  RETURN
284 *
285 * End of ZGEQLF
286 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: zlarfb.f:197
subroutine zlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
ZLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: zlarft.f:165
subroutine zgeql2(M, N, A, LDA, TAU, WORK, INFO)
ZGEQL2 computes the QL factorization of a general rectangular matrix using an unblocked algorithm...
Definition: zgeql2.f:125
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeqp3 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGEQP3

Download ZGEQP3 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQP3 computes a QR factorization with column pivoting of a
 matrix A:  A*P = Q*R  using Level 3 BLAS.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A. M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the upper triangle of the array contains the
          min(M,N)-by-N upper trapezoidal matrix R; the elements below
          the diagonal, together with the array TAU, represent the
          unitary matrix Q as a product of min(M,N) elementary
          reflectors.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[in,out]JPVT
          JPVT is INTEGER array, dimension (N)
          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
          to the front of A*P (a leading column); if JPVT(J)=0,
          the J-th column of A is a free column.
          On exit, if JPVT(J)=K, then the J-th column of A*P was the
          the K-th column of A.
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= N+1.
          For optimal performance LWORK >= ( N+1 )*NB, where NB
          is the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (2*N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit.
          < 0: if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(1) H(2) . . . H(k), where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a real/complex vector
  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
  A(i+1:m,i), and tau in TAU(i).
Contributors:
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain X. Sun, Computer Science Dept., Duke University, USA

Definition at line 161 of file zgeqp3.f.

161 *
162 * -- LAPACK computational routine (version 3.6.0) --
163 * -- LAPACK is a software package provided by Univ. of Tennessee, --
164 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165 * November 2015
166 *
167 * .. Scalar Arguments ..
168  INTEGER info, lda, lwork, m, n
169 * ..
170 * .. Array Arguments ..
171  INTEGER jpvt( * )
172  DOUBLE PRECISION rwork( * )
173  COMPLEX*16 a( lda, * ), tau( * ), work( * )
174 * ..
175 *
176 * =====================================================================
177 *
178 * .. Parameters ..
179  INTEGER inb, inbmin, ixover
180  parameter( inb = 1, inbmin = 2, ixover = 3 )
181 * ..
182 * .. Local Scalars ..
183  LOGICAL lquery
184  INTEGER fjb, iws, j, jb, lwkopt, minmn, minws, na, nb,
185  $ nbmin, nfxd, nx, sm, sminmn, sn, topbmn
186 * ..
187 * .. External Subroutines ..
188  EXTERNAL xerbla, zgeqrf, zlaqp2, zlaqps, zswap, zunmqr
189 * ..
190 * .. External Functions ..
191  INTEGER ilaenv
192  DOUBLE PRECISION dznrm2
193  EXTERNAL ilaenv, dznrm2
194 * ..
195 * .. Intrinsic Functions ..
196  INTRINSIC int, max, min
197 * ..
198 * .. Executable Statements ..
199 *
200 * Test input arguments
201 * ====================
202 *
203  info = 0
204  lquery = ( lwork.EQ.-1 )
205  IF( m.LT.0 ) THEN
206  info = -1
207  ELSE IF( n.LT.0 ) THEN
208  info = -2
209  ELSE IF( lda.LT.max( 1, m ) ) THEN
210  info = -4
211  END IF
212 *
213  IF( info.EQ.0 ) THEN
214  minmn = min( m, n )
215  IF( minmn.EQ.0 ) THEN
216  iws = 1
217  lwkopt = 1
218  ELSE
219  iws = n + 1
220  nb = ilaenv( inb, 'ZGEQRF', ' ', m, n, -1, -1 )
221  lwkopt = ( n + 1 )*nb
222  END IF
223  work( 1 ) = dcmplx( lwkopt )
224 *
225  IF( ( lwork.LT.iws ) .AND. .NOT.lquery ) THEN
226  info = -8
227  END IF
228  END IF
229 *
230  IF( info.NE.0 ) THEN
231  CALL xerbla( 'ZGEQP3', -info )
232  RETURN
233  ELSE IF( lquery ) THEN
234  RETURN
235  END IF
236 *
237 * Move initial columns up front.
238 *
239  nfxd = 1
240  DO 10 j = 1, n
241  IF( jpvt( j ).NE.0 ) THEN
242  IF( j.NE.nfxd ) THEN
243  CALL zswap( m, a( 1, j ), 1, a( 1, nfxd ), 1 )
244  jpvt( j ) = jpvt( nfxd )
245  jpvt( nfxd ) = j
246  ELSE
247  jpvt( j ) = j
248  END IF
249  nfxd = nfxd + 1
250  ELSE
251  jpvt( j ) = j
252  END IF
253  10 CONTINUE
254  nfxd = nfxd - 1
255 *
256 * Factorize fixed columns
257 * =======================
258 *
259 * Compute the QR factorization of fixed columns and update
260 * remaining columns.
261 *
262  IF( nfxd.GT.0 ) THEN
263  na = min( m, nfxd )
264 *CC CALL ZGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
265  CALL zgeqrf( m, na, a, lda, tau, work, lwork, info )
266  iws = max( iws, int( work( 1 ) ) )
267  IF( na.LT.n ) THEN
268 *CC CALL ZUNM2R( 'Left', 'Conjugate Transpose', M, N-NA,
269 *CC $ NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK,
270 *CC $ INFO )
271  CALL zunmqr( 'Left', 'Conjugate Transpose', m, n-na, na, a,
272  $ lda, tau, a( 1, na+1 ), lda, work, lwork,
273  $ info )
274  iws = max( iws, int( work( 1 ) ) )
275  END IF
276  END IF
277 *
278 * Factorize free columns
279 * ======================
280 *
281  IF( nfxd.LT.minmn ) THEN
282 *
283  sm = m - nfxd
284  sn = n - nfxd
285  sminmn = minmn - nfxd
286 *
287 * Determine the block size.
288 *
289  nb = ilaenv( inb, 'ZGEQRF', ' ', sm, sn, -1, -1 )
290  nbmin = 2
291  nx = 0
292 *
293  IF( ( nb.GT.1 ) .AND. ( nb.LT.sminmn ) ) THEN
294 *
295 * Determine when to cross over from blocked to unblocked code.
296 *
297  nx = max( 0, ilaenv( ixover, 'ZGEQRF', ' ', sm, sn, -1,
298  $ -1 ) )
299 *
300 *
301  IF( nx.LT.sminmn ) THEN
302 *
303 * Determine if workspace is large enough for blocked code.
304 *
305  minws = ( sn+1 )*nb
306  iws = max( iws, minws )
307  IF( lwork.LT.minws ) THEN
308 *
309 * Not enough workspace to use optimal NB: Reduce NB and
310 * determine the minimum value of NB.
311 *
312  nb = lwork / ( sn+1 )
313  nbmin = max( 2, ilaenv( inbmin, 'ZGEQRF', ' ', sm, sn,
314  $ -1, -1 ) )
315 *
316 *
317  END IF
318  END IF
319  END IF
320 *
321 * Initialize partial column norms. The first N elements of work
322 * store the exact column norms.
323 *
324  DO 20 j = nfxd + 1, n
325  rwork( j ) = dznrm2( sm, a( nfxd+1, j ), 1 )
326  rwork( n+j ) = rwork( j )
327  20 CONTINUE
328 *
329  IF( ( nb.GE.nbmin ) .AND. ( nb.LT.sminmn ) .AND.
330  $ ( nx.LT.sminmn ) ) THEN
331 *
332 * Use blocked code initially.
333 *
334  j = nfxd + 1
335 *
336 * Compute factorization: while loop.
337 *
338 *
339  topbmn = minmn - nx
340  30 CONTINUE
341  IF( j.LE.topbmn ) THEN
342  jb = min( nb, topbmn-j+1 )
343 *
344 * Factorize JB columns among columns J:N.
345 *
346  CALL zlaqps( m, n-j+1, j-1, jb, fjb, a( 1, j ), lda,
347  $ jpvt( j ), tau( j ), rwork( j ),
348  $ rwork( n+j ), work( 1 ), work( jb+1 ),
349  $ n-j+1 )
350 *
351  j = j + fjb
352  GO TO 30
353  END IF
354  ELSE
355  j = nfxd + 1
356  END IF
357 *
358 * Use unblocked code to factor the last or only block.
359 *
360 *
361  IF( j.LE.minmn )
362  $ CALL zlaqp2( m, n-j+1, j-1, a( 1, j ), lda, jpvt( j ),
363  $ tau( j ), rwork( j ), rwork( n+j ), work( 1 ) )
364 *
365  END IF
366 *
367  work( 1 ) = dcmplx( lwkopt )
368  RETURN
369 *
370 * End of ZGEQP3
371 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dznrm2(N, X, INCX)
DZNRM2
Definition: dznrm2.f:56
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zlaqp2(M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)
ZLAQP2 computes a QR factorization with column pivoting of the matrix block.
Definition: zlaqp2.f:151
subroutine zlaqps(M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)
ZLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BL...
Definition: zlaqps.f:179

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeqpf ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGEQPF

Download ZGEQPF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 This routine is deprecated and has been replaced by routine ZGEQP3.

 ZGEQPF computes a QR factorization with column pivoting of a
 complex M-by-N matrix A: A*P = Q*R.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A. M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A. N >= 0
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the upper triangle of the array contains the
          min(M,N)-by-N upper triangular matrix R; the elements
          below the diagonal, together with the array TAU,
          represent the unitary matrix Q as a product of
          min(m,n) elementary reflectors.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[in,out]JPVT
          JPVT is INTEGER array, dimension (N)
          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
          to the front of A*P (a leading column); if JPVT(i) = 0,
          the i-th column of A is a free column.
          On exit, if JPVT(i) = k, then the i-th column of A*P
          was the k-th column of A.
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[out]WORK
          WORK is COMPLEX*16 array, dimension (N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (2*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(1) H(2) . . . H(n)

  Each H(i) has the form

     H = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).

  The matrix P is represented in jpvt as follows: If
     jpvt(j) = i
  then the jth column of P is the ith canonical unit vector.

  Partial column norm updating strategy modified by
    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
    University of Zagreb, Croatia.
  -- April 2011                                                      --
  For more details see LAPACK Working Note 176.

Definition at line 150 of file zgeqpf.f.

150 *
151 * -- LAPACK computational routine (version 3.4.0) --
152 * -- LAPACK is a software package provided by Univ. of Tennessee, --
153 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
154 * November 2011
155 *
156 * .. Scalar Arguments ..
157  INTEGER info, lda, m, n
158 * ..
159 * .. Array Arguments ..
160  INTEGER jpvt( * )
161  DOUBLE PRECISION rwork( * )
162  COMPLEX*16 a( lda, * ), tau( * ), work( * )
163 * ..
164 *
165 * =====================================================================
166 *
167 * .. Parameters ..
168  DOUBLE PRECISION zero, one
169  parameter( zero = 0.0d+0, one = 1.0d+0 )
170 * ..
171 * .. Local Scalars ..
172  INTEGER i, itemp, j, ma, mn, pvt
173  DOUBLE PRECISION temp, temp2, tol3z
174  COMPLEX*16 aii
175 * ..
176 * .. External Subroutines ..
177  EXTERNAL xerbla, zgeqr2, zlarf, zlarfg, zswap, zunm2r
178 * ..
179 * .. Intrinsic Functions ..
180  INTRINSIC abs, dcmplx, dconjg, max, min, sqrt
181 * ..
182 * .. External Functions ..
183  INTEGER idamax
184  DOUBLE PRECISION dlamch, dznrm2
185  EXTERNAL idamax, dlamch, dznrm2
186 * ..
187 * .. Executable Statements ..
188 *
189 * Test the input arguments
190 *
191  info = 0
192  IF( m.LT.0 ) THEN
193  info = -1
194  ELSE IF( n.LT.0 ) THEN
195  info = -2
196  ELSE IF( lda.LT.max( 1, m ) ) THEN
197  info = -4
198  END IF
199  IF( info.NE.0 ) THEN
200  CALL xerbla( 'ZGEQPF', -info )
201  RETURN
202  END IF
203 *
204  mn = min( m, n )
205  tol3z = sqrt(dlamch('Epsilon'))
206 *
207 * Move initial columns up front
208 *
209  itemp = 1
210  DO 10 i = 1, n
211  IF( jpvt( i ).NE.0 ) THEN
212  IF( i.NE.itemp ) THEN
213  CALL zswap( m, a( 1, i ), 1, a( 1, itemp ), 1 )
214  jpvt( i ) = jpvt( itemp )
215  jpvt( itemp ) = i
216  ELSE
217  jpvt( i ) = i
218  END IF
219  itemp = itemp + 1
220  ELSE
221  jpvt( i ) = i
222  END IF
223  10 CONTINUE
224  itemp = itemp - 1
225 *
226 * Compute the QR factorization and update remaining columns
227 *
228  IF( itemp.GT.0 ) THEN
229  ma = min( itemp, m )
230  CALL zgeqr2( m, ma, a, lda, tau, work, info )
231  IF( ma.LT.n ) THEN
232  CALL zunm2r( 'Left', 'Conjugate transpose', m, n-ma, ma, a,
233  $ lda, tau, a( 1, ma+1 ), lda, work, info )
234  END IF
235  END IF
236 *
237  IF( itemp.LT.mn ) THEN
238 *
239 * Initialize partial column norms. The first n elements of
240 * work store the exact column norms.
241 *
242  DO 20 i = itemp + 1, n
243  rwork( i ) = dznrm2( m-itemp, a( itemp+1, i ), 1 )
244  rwork( n+i ) = rwork( i )
245  20 CONTINUE
246 *
247 * Compute factorization
248 *
249  DO 40 i = itemp + 1, mn
250 *
251 * Determine ith pivot column and swap if necessary
252 *
253  pvt = ( i-1 ) + idamax( n-i+1, rwork( i ), 1 )
254 *
255  IF( pvt.NE.i ) THEN
256  CALL zswap( m, a( 1, pvt ), 1, a( 1, i ), 1 )
257  itemp = jpvt( pvt )
258  jpvt( pvt ) = jpvt( i )
259  jpvt( i ) = itemp
260  rwork( pvt ) = rwork( i )
261  rwork( n+pvt ) = rwork( n+i )
262  END IF
263 *
264 * Generate elementary reflector H(i)
265 *
266  aii = a( i, i )
267  CALL zlarfg( m-i+1, aii, a( min( i+1, m ), i ), 1,
268  $ tau( i ) )
269  a( i, i ) = aii
270 *
271  IF( i.LT.n ) THEN
272 *
273 * Apply H(i) to A(i:m,i+1:n) from the left
274 *
275  aii = a( i, i )
276  a( i, i ) = dcmplx( one )
277  CALL zlarf( 'Left', m-i+1, n-i, a( i, i ), 1,
278  $ dconjg( tau( i ) ), a( i, i+1 ), lda, work )
279  a( i, i ) = aii
280  END IF
281 *
282 * Update partial column norms
283 *
284  DO 30 j = i + 1, n
285  IF( rwork( j ).NE.zero ) THEN
286 *
287 * NOTE: The following 4 lines follow from the analysis in
288 * Lapack Working Note 176.
289 *
290  temp = abs( a( i, j ) ) / rwork( j )
291  temp = max( zero, ( one+temp )*( one-temp ) )
292  temp2 = temp*( rwork( j ) / rwork( n+j ) )**2
293  IF( temp2 .LE. tol3z ) THEN
294  IF( m-i.GT.0 ) THEN
295  rwork( j ) = dznrm2( m-i, a( i+1, j ), 1 )
296  rwork( n+j ) = rwork( j )
297  ELSE
298  rwork( j ) = zero
299  rwork( n+j ) = zero
300  END IF
301  ELSE
302  rwork( j ) = rwork( j )*sqrt( temp )
303  END IF
304  END IF
305  30 CONTINUE
306 *
307  40 CONTINUE
308  END IF
309  RETURN
310 *
311 * End of ZGEQPF
312 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgeqr2(M, N, A, LDA, TAU, WORK, INFO)
ZGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: zgeqr2.f:123
subroutine zunm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
ZUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition: zunm2r.f:161
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
double precision function dznrm2(N, X, INCX)
DZNRM2
Definition: dznrm2.f:56
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeqr2 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.

Download ZGEQR2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQR2 computes a QR factorization of a complex m by n matrix A:
 A = Q * R.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the m by n matrix A.
          On exit, the elements on and above the diagonal of the array
          contain the min(m,n) by n upper trapezoidal matrix R (R is
          upper triangular if m >= n); the elements below the diagonal,
          with the array TAU, represent the unitary matrix Q as a
          product of elementary reflectors (see Further Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(1) H(2) . . . H(k), where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
  and tau in TAU(i).

Definition at line 123 of file zgeqr2.f.

123 *
124 * -- LAPACK computational routine (version 3.4.2) --
125 * -- LAPACK is a software package provided by Univ. of Tennessee, --
126 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
127 * September 2012
128 *
129 * .. Scalar Arguments ..
130  INTEGER info, lda, m, n
131 * ..
132 * .. Array Arguments ..
133  COMPLEX*16 a( lda, * ), tau( * ), work( * )
134 * ..
135 *
136 * =====================================================================
137 *
138 * .. Parameters ..
139  COMPLEX*16 one
140  parameter( one = ( 1.0d+0, 0.0d+0 ) )
141 * ..
142 * .. Local Scalars ..
143  INTEGER i, k
144  COMPLEX*16 alpha
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL xerbla, zlarf, zlarfg
148 * ..
149 * .. Intrinsic Functions ..
150  INTRINSIC dconjg, max, min
151 * ..
152 * .. Executable Statements ..
153 *
154 * Test the input arguments
155 *
156  info = 0
157  IF( m.LT.0 ) THEN
158  info = -1
159  ELSE IF( n.LT.0 ) THEN
160  info = -2
161  ELSE IF( lda.LT.max( 1, m ) ) THEN
162  info = -4
163  END IF
164  IF( info.NE.0 ) THEN
165  CALL xerbla( 'ZGEQR2', -info )
166  RETURN
167  END IF
168 *
169  k = min( m, n )
170 *
171  DO 10 i = 1, k
172 *
173 * Generate elementary reflector H(i) to annihilate A(i+1:m,i)
174 *
175  CALL zlarfg( m-i+1, a( i, i ), a( min( i+1, m ), i ), 1,
176  $ tau( i ) )
177  IF( i.LT.n ) THEN
178 *
179 * Apply H(i)**H to A(i:m,i+1:n) from the left
180 *
181  alpha = a( i, i )
182  a( i, i ) = one
183  CALL zlarf( 'Left', m-i+1, n-i, a( i, i ), 1,
184  $ dconjg( tau( i ) ), a( i, i+1 ), lda, work )
185  a( i, i ) = alpha
186  END IF
187  10 CONTINUE
188  RETURN
189 *
190 * End of ZGEQR2
191 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeqr2p ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZGEQR2P computes the QR factorization of a general rectangular matrix with non-negative diagonal elements using an unblocked algorithm.

Download ZGEQR2P + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQR2P computes a QR factorization of a complex m by n matrix A:
 A = Q * R. The diagonal entries of R are real and nonnegative.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the m by n matrix A.
          On exit, the elements on and above the diagonal of the array
          contain the min(m,n) by n upper trapezoidal matrix R (R is
          upper triangular if m >= n). The diagonal entries of R 
          are real and nonnegative; the elements below the diagonal,
          with the array TAU, represent the unitary matrix Q as a
          product of elementary reflectors (see Further Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(1) H(2) . . . H(k), where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
  and tau in TAU(i).

 See Lapack Working Note 203 for details

Definition at line 126 of file zgeqr2p.f.

126 *
127 * -- LAPACK computational routine (version 3.6.0) --
128 * -- LAPACK is a software package provided by Univ. of Tennessee, --
129 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
130 * November 2015
131 *
132 * .. Scalar Arguments ..
133  INTEGER info, lda, m, n
134 * ..
135 * .. Array Arguments ..
136  COMPLEX*16 a( lda, * ), tau( * ), work( * )
137 * ..
138 *
139 * =====================================================================
140 *
141 * .. Parameters ..
142  COMPLEX*16 one
143  parameter( one = ( 1.0d+0, 0.0d+0 ) )
144 * ..
145 * .. Local Scalars ..
146  INTEGER i, k
147  COMPLEX*16 alpha
148 * ..
149 * .. External Subroutines ..
150  EXTERNAL xerbla, zlarf, zlarfgp
151 * ..
152 * .. Intrinsic Functions ..
153  INTRINSIC dconjg, max, min
154 * ..
155 * .. Executable Statements ..
156 *
157 * Test the input arguments
158 *
159  info = 0
160  IF( m.LT.0 ) THEN
161  info = -1
162  ELSE IF( n.LT.0 ) THEN
163  info = -2
164  ELSE IF( lda.LT.max( 1, m ) ) THEN
165  info = -4
166  END IF
167  IF( info.NE.0 ) THEN
168  CALL xerbla( 'ZGEQR2P', -info )
169  RETURN
170  END IF
171 *
172  k = min( m, n )
173 *
174  DO 10 i = 1, k
175 *
176 * Generate elementary reflector H(i) to annihilate A(i+1:m,i)
177 *
178  CALL zlarfgp( m-i+1, a( i, i ), a( min( i+1, m ), i ), 1,
179  $ tau( i ) )
180  IF( i.LT.n ) THEN
181 *
182 * Apply H(i)**H to A(i:m,i+1:n) from the left
183 *
184  alpha = a( i, i )
185  a( i, i ) = one
186  CALL zlarf( 'Left', m-i+1, n-i, a( i, i ), 1,
187  $ dconjg( tau( i ) ), a( i, i+1 ), lda, work )
188  a( i, i ) = alpha
189  END IF
190  10 CONTINUE
191  RETURN
192 *
193 * End of ZGEQR2P
194 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarfgp(N, ALPHA, X, INCX, TAU)
ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: zlarfgp.f:106
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeqrf ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGEQRF

Download ZGEQRF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQRF computes a QR factorization of a complex M-by-N matrix A:
 A = Q * R.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the elements on and above the diagonal of the array
          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
          upper triangular if m >= n); the elements below the diagonal,
          with the array TAU, represent the unitary matrix Q as a
          product of min(m,n) elementary reflectors (see Further
          Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,N).
          For optimum performance LWORK >= N*NB, where NB is
          the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(1) H(2) . . . H(k), where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
  and tau in TAU(i).

Definition at line 138 of file zgeqrf.f.

138 *
139 * -- LAPACK computational routine (version 3.4.0) --
140 * -- LAPACK is a software package provided by Univ. of Tennessee, --
141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142 * November 2011
143 *
144 * .. Scalar Arguments ..
145  INTEGER info, lda, lwork, m, n
146 * ..
147 * .. Array Arguments ..
148  COMPLEX*16 a( lda, * ), tau( * ), work( * )
149 * ..
150 *
151 * =====================================================================
152 *
153 * .. Local Scalars ..
154  LOGICAL lquery
155  INTEGER i, ib, iinfo, iws, k, ldwork, lwkopt, nb,
156  $ nbmin, nx
157 * ..
158 * .. External Subroutines ..
159  EXTERNAL xerbla, zgeqr2, zlarfb, zlarft
160 * ..
161 * .. Intrinsic Functions ..
162  INTRINSIC max, min
163 * ..
164 * .. External Functions ..
165  INTEGER ilaenv
166  EXTERNAL ilaenv
167 * ..
168 * .. Executable Statements ..
169 *
170 * Test the input arguments
171 *
172  info = 0
173  nb = ilaenv( 1, 'ZGEQRF', ' ', m, n, -1, -1 )
174  lwkopt = n*nb
175  work( 1 ) = lwkopt
176  lquery = ( lwork.EQ.-1 )
177  IF( m.LT.0 ) THEN
178  info = -1
179  ELSE IF( n.LT.0 ) THEN
180  info = -2
181  ELSE IF( lda.LT.max( 1, m ) ) THEN
182  info = -4
183  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
184  info = -7
185  END IF
186  IF( info.NE.0 ) THEN
187  CALL xerbla( 'ZGEQRF', -info )
188  RETURN
189  ELSE IF( lquery ) THEN
190  RETURN
191  END IF
192 *
193 * Quick return if possible
194 *
195  k = min( m, n )
196  IF( k.EQ.0 ) THEN
197  work( 1 ) = 1
198  RETURN
199  END IF
200 *
201  nbmin = 2
202  nx = 0
203  iws = n
204  IF( nb.GT.1 .AND. nb.LT.k ) THEN
205 *
206 * Determine when to cross over from blocked to unblocked code.
207 *
208  nx = max( 0, ilaenv( 3, 'ZGEQRF', ' ', m, n, -1, -1 ) )
209  IF( nx.LT.k ) THEN
210 *
211 * Determine if workspace is large enough for blocked code.
212 *
213  ldwork = n
214  iws = ldwork*nb
215  IF( lwork.LT.iws ) THEN
216 *
217 * Not enough workspace to use optimal NB: reduce NB and
218 * determine the minimum value of NB.
219 *
220  nb = lwork / ldwork
221  nbmin = max( 2, ilaenv( 2, 'ZGEQRF', ' ', m, n, -1,
222  $ -1 ) )
223  END IF
224  END IF
225  END IF
226 *
227  IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
228 *
229 * Use blocked code initially
230 *
231  DO 10 i = 1, k - nx, nb
232  ib = min( k-i+1, nb )
233 *
234 * Compute the QR factorization of the current block
235 * A(i:m,i:i+ib-1)
236 *
237  CALL zgeqr2( m-i+1, ib, a( i, i ), lda, tau( i ), work,
238  $ iinfo )
239  IF( i+ib.LE.n ) THEN
240 *
241 * Form the triangular factor of the block reflector
242 * H = H(i) H(i+1) . . . H(i+ib-1)
243 *
244  CALL zlarft( 'Forward', 'Columnwise', m-i+1, ib,
245  $ a( i, i ), lda, tau( i ), work, ldwork )
246 *
247 * Apply H**H to A(i:m,i+ib:n) from the left
248 *
249  CALL zlarfb( 'Left', 'Conjugate transpose', 'Forward',
250  $ 'Columnwise', m-i+1, n-i-ib+1, ib,
251  $ a( i, i ), lda, work, ldwork, a( i, i+ib ),
252  $ lda, work( ib+1 ), ldwork )
253  END IF
254  10 CONTINUE
255  ELSE
256  i = 1
257  END IF
258 *
259 * Use unblocked code to factor the last or only block.
260 *
261  IF( i.LE.k )
262  $ CALL zgeqr2( m-i+1, n-i+1, a( i, i ), lda, tau( i ), work,
263  $ iinfo )
264 *
265  work( 1 ) = iws
266  RETURN
267 *
268 * End of ZGEQRF
269 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: zlarfb.f:197
subroutine zgeqr2(M, N, A, LDA, TAU, WORK, INFO)
ZGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: zgeqr2.f:123
subroutine zlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
ZLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: zlarft.f:165
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

subroutine zgeqrfp ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGEQRFP

Download ZGEQRFP + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQRFP computes a QR factorization of a complex M-by-N matrix A:
 A = Q * R. The diagonal entries of R are real and nonnegative.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the elements on and above the diagonal of the array
          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
          upper triangular if m >= n). The diagonal entries of R
          are real and nonnegative; The elements below the diagonal,
          with the array TAU, represent the unitary matrix Q as a
          product of min(m,n) elementary reflectors (see Further
          Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,N).
          For optimum performance LWORK >= N*NB, where NB is
          the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(1) H(2) . . . H(k), where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
  and tau in TAU(i).

 See Lapack Working Note 203 for details

Definition at line 141 of file zgeqrfp.f.

141 *
142 * -- LAPACK computational routine (version 3.6.0) --
143 * -- LAPACK is a software package provided by Univ. of Tennessee, --
144 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
145 * November 2015
146 *
147 * .. Scalar Arguments ..
148  INTEGER info, lda, lwork, m, n
149 * ..
150 * .. Array Arguments ..
151  COMPLEX*16 a( lda, * ), tau( * ), work( * )
152 * ..
153 *
154 * =====================================================================
155 *
156 * .. Local Scalars ..
157  LOGICAL lquery
158  INTEGER i, ib, iinfo, iws, k, ldwork, lwkopt, nb,
159  $ nbmin, nx
160 * ..
161 * .. External Subroutines ..
162  EXTERNAL xerbla, zgeqr2p, zlarfb, zlarft
163 * ..
164 * .. Intrinsic Functions ..
165  INTRINSIC max, min
166 * ..
167 * .. External Functions ..
168  INTEGER ilaenv
169  EXTERNAL ilaenv
170 * ..
171 * .. Executable Statements ..
172 *
173 * Test the input arguments
174 *
175  info = 0
176  nb = ilaenv( 1, 'ZGEQRF', ' ', m, n, -1, -1 )
177  lwkopt = n*nb
178  work( 1 ) = lwkopt
179  lquery = ( lwork.EQ.-1 )
180  IF( m.LT.0 ) THEN
181  info = -1
182  ELSE IF( n.LT.0 ) THEN
183  info = -2
184  ELSE IF( lda.LT.max( 1, m ) ) THEN
185  info = -4
186  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
187  info = -7
188  END IF
189  IF( info.NE.0 ) THEN
190  CALL xerbla( 'ZGEQRFP', -info )
191  RETURN
192  ELSE IF( lquery ) THEN
193  RETURN
194  END IF
195 *
196 * Quick return if possible
197 *
198  k = min( m, n )
199  IF( k.EQ.0 ) THEN
200  work( 1 ) = 1
201  RETURN
202  END IF
203 *
204  nbmin = 2
205  nx = 0
206  iws = n
207  IF( nb.GT.1 .AND. nb.LT.k ) THEN
208 *
209 * Determine when to cross over from blocked to unblocked code.
210 *
211  nx = max( 0, ilaenv( 3, 'ZGEQRF', ' ', m, n, -1, -1 ) )
212  IF( nx.LT.k ) THEN
213 *
214 * Determine if workspace is large enough for blocked code.
215 *
216  ldwork = n
217  iws = ldwork*nb
218  IF( lwork.LT.iws ) THEN
219 *
220 * Not enough workspace to use optimal NB: reduce NB and
221 * determine the minimum value of NB.
222 *
223  nb = lwork / ldwork
224  nbmin = max( 2, ilaenv( 2, 'ZGEQRF', ' ', m, n, -1,
225  $ -1 ) )
226  END IF
227  END IF
228  END IF
229 *
230  IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
231 *
232 * Use blocked code initially
233 *
234  DO 10 i = 1, k - nx, nb
235  ib = min( k-i+1, nb )
236 *
237 * Compute the QR factorization of the current block
238 * A(i:m,i:i+ib-1)
239 *
240  CALL zgeqr2p( m-i+1, ib, a( i, i ), lda, tau( i ), work,
241  $ iinfo )
242  IF( i+ib.LE.n ) THEN
243 *
244 * Form the triangular factor of the block reflector
245 * H = H(i) H(i+1) . . . H(i+ib-1)
246 *
247  CALL zlarft( 'Forward', 'Columnwise', m-i+1, ib,
248  $ a( i, i ), lda, tau( i ), work, ldwork )
249 *
250 * Apply H**H to A(i:m,i+ib:n) from the left
251 *
252  CALL zlarfb( 'Left', 'Conjugate transpose', 'Forward',
253  $ 'Columnwise', m-i+1, n-i-ib+1, ib,
254  $ a( i, i ), lda, work, ldwork, a( i, i+ib ),
255  $ lda, work( ib+1 ), ldwork )
256  END IF
257  10 CONTINUE
258  ELSE
259  i = 1
260  END IF
261 *
262 * Use unblocked code to factor the last or only block.
263 *
264  IF( i.LE.k )
265  $ CALL zgeqr2p( m-i+1, n-i+1, a( i, i ), lda, tau( i ), work,
266  $ iinfo )
267 *
268  work( 1 ) = iws
269  RETURN
270 *
271 * End of ZGEQRFP
272 *
subroutine zgeqr2p(M, N, A, LDA, TAU, WORK, INFO)
ZGEQR2P computes the QR factorization of a general rectangular matrix with non-negative diagonal elem...
Definition: zgeqr2p.f:126
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: zlarfb.f:197
subroutine zlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
ZLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: zlarft.f:165
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeqrt ( integer  M,
integer  N,
integer  NB,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldt, * )  T,
integer  LDT,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZGEQRT

Download ZGEQRT + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQRT computes a blocked QR factorization of a complex M-by-N matrix A
 using the compact WY representation of Q.  
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]NB
          NB is INTEGER
          The block size to be used in the blocked QR.  MIN(M,N) >= NB >= 1.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the elements on and above the diagonal of the array
          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
          upper triangular if M >= N); the elements below the diagonal
          are the columns of V.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]T
          T is COMPLEX*16 array, dimension (LDT,MIN(M,N))
          The upper triangular block reflectors stored in compact form
          as a sequence of upper triangular blocks.  See below
          for further details.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T.  LDT >= NB.
[out]WORK
          WORK is COMPLEX*16 array, dimension (NB*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Further Details:
  The matrix V stores the elementary reflectors H(i) in the i-th column
  below the diagonal. For example, if M=5 and N=3, the matrix V is

               V = (  1       )
                   ( v1  1    )
                   ( v1 v2  1 )
                   ( v1 v2 v3 )
                   ( v1 v2 v3 )

  where the vi's represent the vectors which define H(i), which are returned
  in the matrix A.  The 1's along the diagonal of V are not stored in A.

  Let K=MIN(M,N).  The number of blocks is B = ceiling(K/NB), where each
  block is of order NB except for the last block, which is of order 
  IB = K - (B-1)*NB.  For each of the B blocks, a upper triangular block
  reflector factor is computed: T1, T2, ..., TB.  The NB-by-NB (and IB-by-IB 
  for the last block) T's are stored in the NB-by-N matrix T as

               T = (T1 T2 ... TB).

Definition at line 143 of file zgeqrt.f.

143 *
144 * -- LAPACK computational routine (version 3.5.0) --
145 * -- LAPACK is a software package provided by Univ. of Tennessee, --
146 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
147 * November 2013
148 *
149 * .. Scalar Arguments ..
150  INTEGER info, lda, ldt, m, n, nb
151 * ..
152 * .. Array Arguments ..
153  COMPLEX*16 a( lda, * ), t( ldt, * ), work( * )
154 * ..
155 *
156 * =====================================================================
157 *
158 * ..
159 * .. Local Scalars ..
160  INTEGER i, ib, iinfo, k
161  LOGICAL use_recursive_qr
162  parameter( use_recursive_qr=.true. )
163 * ..
164 * .. External Subroutines ..
165  EXTERNAL zgeqrt2, zgeqrt3, zlarfb, xerbla
166 * ..
167 * .. Executable Statements ..
168 *
169 * Test the input arguments
170 *
171  info = 0
172  IF( m.LT.0 ) THEN
173  info = -1
174  ELSE IF( n.LT.0 ) THEN
175  info = -2
176  ELSE IF( nb.LT.1 .OR. ( nb.GT.min(m,n) .AND. min(m,n).GT.0 ) )THEN
177  info = -3
178  ELSE IF( lda.LT.max( 1, m ) ) THEN
179  info = -5
180  ELSE IF( ldt.LT.nb ) THEN
181  info = -7
182  END IF
183  IF( info.NE.0 ) THEN
184  CALL xerbla( 'ZGEQRT', -info )
185  RETURN
186  END IF
187 *
188 * Quick return if possible
189 *
190  k = min( m, n )
191  IF( k.EQ.0 ) RETURN
192 *
193 * Blocked loop of length K
194 *
195  DO i = 1, k, nb
196  ib = min( k-i+1, nb )
197 *
198 * Compute the QR factorization of the current block A(I:M,I:I+IB-1)
199 *
200  IF( use_recursive_qr ) THEN
201  CALL zgeqrt3( m-i+1, ib, a(i,i), lda, t(1,i), ldt, iinfo )
202  ELSE
203  CALL zgeqrt2( m-i+1, ib, a(i,i), lda, t(1,i), ldt, iinfo )
204  END IF
205  IF( i+ib.LE.n ) THEN
206 *
207 * Update by applying H**H to A(I:M,I+IB:N) from the left
208 *
209  CALL zlarfb( 'L', 'C', 'F', 'C', m-i+1, n-i-ib+1, ib,
210  $ a( i, i ), lda, t( 1, i ), ldt,
211  $ a( i, i+ib ), lda, work , n-i-ib+1 )
212  END IF
213  END DO
214  RETURN
215 *
216 * End of ZGEQRT
217 *
subroutine zgeqrt2(M, N, A, LDA, T, LDT, INFO)
ZGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY represen...
Definition: zgeqrt2.f:129
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: zlarfb.f:197
recursive subroutine zgeqrt3(M, N, A, LDA, T, LDT, INFO)
ZGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact...
Definition: zgeqrt3.f:134

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgeqrt2 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldt, * )  T,
integer  LDT,
integer  INFO 
)

ZGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY representation of Q.

Download ZGEQRT2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQRT2 computes a QR factorization of a complex M-by-N matrix A, 
 using the compact WY representation of Q. 
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= N.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the complex M-by-N matrix A.  On exit, the elements on and
          above the diagonal contain the N-by-N upper triangular matrix R; the
          elements below the diagonal are the columns of V.  See below for
          further details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]T
          T is COMPLEX*16 array, dimension (LDT,N)
          The N-by-N upper triangular factor of the block reflector.
          The elements on and above the diagonal contain the block
          reflector T; the elements below the diagonal are not used.
          See below for further details.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T.  LDT >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  The matrix V stores the elementary reflectors H(i) in the i-th column
  below the diagonal. For example, if M=5 and N=3, the matrix V is

               V = (  1       )
                   ( v1  1    )
                   ( v1 v2  1 )
                   ( v1 v2 v3 )
                   ( v1 v2 v3 )

  where the vi's represent the vectors which define H(i), which are returned
  in the matrix A.  The 1's along the diagonal of V are not stored in A.  The
  block reflector H is then given by

               H = I - V * T * V**H

  where V**H is the conjugate transpose of V.

Definition at line 129 of file zgeqrt2.f.

129 *
130 * -- LAPACK computational routine (version 3.4.2) --
131 * -- LAPACK is a software package provided by Univ. of Tennessee, --
132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133 * September 2012
134 *
135 * .. Scalar Arguments ..
136  INTEGER info, lda, ldt, m, n
137 * ..
138 * .. Array Arguments ..
139  COMPLEX*16 a( lda, * ), t( ldt, * )
140 * ..
141 *
142 * =====================================================================
143 *
144 * .. Parameters ..
145  COMPLEX*16 one, zero
146  parameter( one = (1.0d+00,0.0d+00), zero = (0.0d+00,0.0d+00) )
147 * ..
148 * .. Local Scalars ..
149  INTEGER i, k
150  COMPLEX*16 aii, alpha
151 * ..
152 * .. External Subroutines ..
153  EXTERNAL zlarfg, zgemv, zgerc, ztrmv, xerbla
154 * ..
155 * .. Executable Statements ..
156 *
157 * Test the input arguments
158 *
159  info = 0
160  IF( m.LT.0 ) THEN
161  info = -1
162  ELSE IF( n.LT.0 ) THEN
163  info = -2
164  ELSE IF( lda.LT.max( 1, m ) ) THEN
165  info = -4
166  ELSE IF( ldt.LT.max( 1, n ) ) THEN
167  info = -6
168  END IF
169  IF( info.NE.0 ) THEN
170  CALL xerbla( 'ZGEQRT2', -info )
171  RETURN
172  END IF
173 *
174  k = min( m, n )
175 *
176  DO i = 1, k
177 *
178 * Generate elem. refl. H(i) to annihilate A(i+1:m,i), tau(I) -> T(I,1)
179 *
180  CALL zlarfg( m-i+1, a( i, i ), a( min( i+1, m ), i ), 1,
181  $ t( i, 1 ) )
182  IF( i.LT.n ) THEN
183 *
184 * Apply H(i) to A(I:M,I+1:N) from the left
185 *
186  aii = a( i, i )
187  a( i, i ) = one
188 *
189 * W(1:N-I) := A(I:M,I+1:N)^H * A(I:M,I) [W = T(:,N)]
190 *
191  CALL zgemv( 'C',m-i+1, n-i, one, a( i, i+1 ), lda,
192  $ a( i, i ), 1, zero, t( 1, n ), 1 )
193 *
194 * A(I:M,I+1:N) = A(I:m,I+1:N) + alpha*A(I:M,I)*W(1:N-1)^H
195 *
196  alpha = -conjg(t( i, 1 ))
197  CALL zgerc( m-i+1, n-i, alpha, a( i, i ), 1,
198  $ t( 1, n ), 1, a( i, i+1 ), lda )
199  a( i, i ) = aii
200  END IF
201  END DO
202 *
203  DO i = 2, n
204  aii = a( i, i )
205  a( i, i ) = one
206 *
207 * T(1:I-1,I) := alpha * A(I:M,1:I-1)**H * A(I:M,I)
208 *
209  alpha = -t( i, 1 )
210  CALL zgemv( 'C', m-i+1, i-1, alpha, a( i, 1 ), lda,
211  $ a( i, i ), 1, zero, t( 1, i ), 1 )
212  a( i, i ) = aii
213 *
214 * T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
215 *
216  CALL ztrmv( 'U', 'N', 'N', i-1, t, ldt, t( 1, i ), 1 )
217 *
218 * T(I,I) = tau(I)
219 *
220  t( i, i ) = t( i, 1 )
221  t( i, 1) = zero
222  END DO
223 
224 *
225 * End of ZGEQRT2
226 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine ztrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
ZTRMV
Definition: ztrmv.f:149
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108
subroutine zgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERC
Definition: zgerc.f:132

Here is the call graph for this function:

Here is the caller graph for this function:

recursive subroutine zgeqrt3 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldt, * )  T,
integer  LDT,
integer  INFO 
)

ZGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact WY representation of Q.

Download ZGEQRT3 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGEQRT3 recursively computes a QR factorization of a complex M-by-N 
 matrix A, using the compact WY representation of Q. 

 Based on the algorithm of Elmroth and Gustavson, 
 IBM J. Res. Develop. Vol 44 No. 4 July 2000.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= N.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the complex M-by-N matrix A.  On exit, the elements on 
          and above the diagonal contain the N-by-N upper triangular matrix R;
          the elements below the diagonal are the columns of V.  See below for
          further details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]T
          T is COMPLEX*16 array, dimension (LDT,N)
          The N-by-N upper triangular factor of the block reflector.
          The elements on and above the diagonal contain the block
          reflector T; the elements below the diagonal are not used.
          See below for further details.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T.  LDT >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  The matrix V stores the elementary reflectors H(i) in the i-th column
  below the diagonal. For example, if M=5 and N=3, the matrix V is

               V = (  1       )
                   ( v1  1    )
                   ( v1 v2  1 )
                   ( v1 v2 v3 )
                   ( v1 v2 v3 )

  where the vi's represent the vectors which define H(i), which are returned
  in the matrix A.  The 1's along the diagonal of V are not stored in A.  The
  block reflector H is then given by

               H = I - V * T * V**H

  where V**H is the conjugate transpose of V.

  For details of the algorithm, see Elmroth and Gustavson (cited above).

Definition at line 134 of file zgeqrt3.f.

134 *
135 * -- LAPACK computational routine (version 3.4.2) --
136 * -- LAPACK is a software package provided by Univ. of Tennessee, --
137 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138 * September 2012
139 *
140 * .. Scalar Arguments ..
141  INTEGER info, lda, m, n, ldt
142 * ..
143 * .. Array Arguments ..
144  COMPLEX*16 a( lda, * ), t( ldt, * )
145 * ..
146 *
147 * =====================================================================
148 *
149 * .. Parameters ..
150  COMPLEX*16 one
151  parameter( one = (1.0d+00,0.0d+00) )
152 * ..
153 * .. Local Scalars ..
154  INTEGER i, i1, j, j1, n1, n2, iinfo
155 * ..
156 * .. External Subroutines ..
157  EXTERNAL zlarfg, ztrmm, zgemm, xerbla
158 * ..
159 * .. Executable Statements ..
160 *
161  info = 0
162  IF( n .LT. 0 ) THEN
163  info = -2
164  ELSE IF( m .LT. n ) THEN
165  info = -1
166  ELSE IF( lda .LT. max( 1, m ) ) THEN
167  info = -4
168  ELSE IF( ldt .LT. max( 1, n ) ) THEN
169  info = -6
170  END IF
171  IF( info.NE.0 ) THEN
172  CALL xerbla( 'ZGEQRT3', -info )
173  RETURN
174  END IF
175 *
176  IF( n.EQ.1 ) THEN
177 *
178 * Compute Householder transform when N=1
179 *
180  CALL zlarfg( m, a, a( min( 2, m ), 1 ), 1, t )
181 *
182  ELSE
183 *
184 * Otherwise, split A into blocks...
185 *
186  n1 = n/2
187  n2 = n-n1
188  j1 = min( n1+1, n )
189  i1 = min( n+1, m )
190 *
191 * Compute A(1:M,1:N1) <- (Y1,R1,T1), where Q1 = I - Y1 T1 Y1^H
192 *
193  CALL zgeqrt3( m, n1, a, lda, t, ldt, iinfo )
194 *
195 * Compute A(1:M,J1:N) = Q1^H A(1:M,J1:N) [workspace: T(1:N1,J1:N)]
196 *
197  DO j=1,n2
198  DO i=1,n1
199  t( i, j+n1 ) = a( i, j+n1 )
200  END DO
201  END DO
202  CALL ztrmm( 'L', 'L', 'C', 'U', n1, n2, one,
203  & a, lda, t( 1, j1 ), ldt )
204 *
205  CALL zgemm( 'C', 'N', n1, n2, m-n1, one, a( j1, 1 ), lda,
206  & a( j1, j1 ), lda, one, t( 1, j1 ), ldt)
207 *
208  CALL ztrmm( 'L', 'U', 'C', 'N', n1, n2, one,
209  & t, ldt, t( 1, j1 ), ldt )
210 *
211  CALL zgemm( 'N', 'N', m-n1, n2, n1, -one, a( j1, 1 ), lda,
212  & t( 1, j1 ), ldt, one, a( j1, j1 ), lda )
213 *
214  CALL ztrmm( 'L', 'L', 'N', 'U', n1, n2, one,
215  & a, lda, t( 1, j1 ), ldt )
216 *
217  DO j=1,n2
218  DO i=1,n1
219  a( i, j+n1 ) = a( i, j+n1 ) - t( i, j+n1 )
220  END DO
221  END DO
222 *
223 * Compute A(J1:M,J1:N) <- (Y2,R2,T2) where Q2 = I - Y2 T2 Y2^H
224 *
225  CALL zgeqrt3( m-n1, n2, a( j1, j1 ), lda,
226  & t( j1, j1 ), ldt, iinfo )
227 *
228 * Compute T3 = T(1:N1,J1:N) = -T1 Y1^H Y2 T2
229 *
230  DO i=1,n1
231  DO j=1,n2
232  t( i, j+n1 ) = conjg(a( j+n1, i ))
233  END DO
234  END DO
235 *
236  CALL ztrmm( 'R', 'L', 'N', 'U', n1, n2, one,
237  & a( j1, j1 ), lda, t( 1, j1 ), ldt )
238 *
239  CALL zgemm( 'C', 'N', n1, n2, m-n, one, a( i1, 1 ), lda,
240  & a( i1, j1 ), lda, one, t( 1, j1 ), ldt )
241 *
242  CALL ztrmm( 'L', 'U', 'N', 'N', n1, n2, -one, t, ldt,
243  & t( 1, j1 ), ldt )
244 *
245  CALL ztrmm( 'R', 'U', 'N', 'N', n1, n2, one,
246  & t( j1, j1 ), ldt, t( 1, j1 ), ldt )
247 *
248 * Y = (Y1,Y2); R = [ R1 A(1:N1,J1:N) ]; T = [T1 T3]
249 * [ 0 R2 ] [ 0 T2]
250 *
251  END IF
252 *
253  RETURN
254 *
255 * End of ZGEQRT3
256 *
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:179
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
recursive subroutine zgeqrt3(M, N, A, LDA, T, LDT, INFO)
ZGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact...
Definition: zgeqrt3.f:134
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgerfs ( character  TRANS,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGERFS

Download ZGERFS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGERFS improves the computed solution to a system of linear
 equations and provides error bounds and backward error estimates for
 the solution.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the form of the system of equations:
          = 'N':  A * X = B     (No transpose)
          = 'T':  A**T * X = B  (Transpose)
          = 'C':  A**H * X = B  (Conjugate transpose)
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X.  NRHS >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The original N-by-N matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
          The factors L and U from the factorization A = P*L*U
          as computed by ZGETRF.
[in]LDAF
          LDAF is INTEGER
          The leading dimension of the array AF.  LDAF >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices from ZGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[in]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
          The right hand side matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[in,out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by ZGETRS.
          On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]FERR
          FERR is DOUBLE PRECISION array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Internal Parameters:
  ITMAX is the maximum number of steps of iterative refinement.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 188 of file zgerfs.f.

188 *
189 * -- LAPACK computational routine (version 3.4.0) --
190 * -- LAPACK is a software package provided by Univ. of Tennessee, --
191 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
192 * November 2011
193 *
194 * .. Scalar Arguments ..
195  CHARACTER trans
196  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs
197 * ..
198 * .. Array Arguments ..
199  INTEGER ipiv( * )
200  DOUBLE PRECISION berr( * ), ferr( * ), rwork( * )
201  COMPLEX*16 a( lda, * ), af( ldaf, * ), b( ldb, * ),
202  $ work( * ), x( ldx, * )
203 * ..
204 *
205 * =====================================================================
206 *
207 * .. Parameters ..
208  INTEGER itmax
209  parameter( itmax = 5 )
210  DOUBLE PRECISION zero
211  parameter( zero = 0.0d+0 )
212  COMPLEX*16 one
213  parameter( one = ( 1.0d+0, 0.0d+0 ) )
214  DOUBLE PRECISION two
215  parameter( two = 2.0d+0 )
216  DOUBLE PRECISION three
217  parameter( three = 3.0d+0 )
218 * ..
219 * .. Local Scalars ..
220  LOGICAL notran
221  CHARACTER transn, transt
222  INTEGER count, i, j, k, kase, nz
223  DOUBLE PRECISION eps, lstres, s, safe1, safe2, safmin, xk
224  COMPLEX*16 zdum
225 * ..
226 * .. Local Arrays ..
227  INTEGER isave( 3 )
228 * ..
229 * .. External Functions ..
230  LOGICAL lsame
231  DOUBLE PRECISION dlamch
232  EXTERNAL lsame, dlamch
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL xerbla, zaxpy, zcopy, zgemv, zgetrs, zlacn2
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, dble, dimag, max
239 * ..
240 * .. Statement Functions ..
241  DOUBLE PRECISION cabs1
242 * ..
243 * .. Statement Function definitions ..
244  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
245 * ..
246 * .. Executable Statements ..
247 *
248 * Test the input parameters.
249 *
250  info = 0
251  notran = lsame( trans, 'N' )
252  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
253  $ lsame( trans, 'C' ) ) THEN
254  info = -1
255  ELSE IF( n.LT.0 ) THEN
256  info = -2
257  ELSE IF( nrhs.LT.0 ) THEN
258  info = -3
259  ELSE IF( lda.LT.max( 1, n ) ) THEN
260  info = -5
261  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
262  info = -7
263  ELSE IF( ldb.LT.max( 1, n ) ) THEN
264  info = -10
265  ELSE IF( ldx.LT.max( 1, n ) ) THEN
266  info = -12
267  END IF
268  IF( info.NE.0 ) THEN
269  CALL xerbla( 'ZGERFS', -info )
270  RETURN
271  END IF
272 *
273 * Quick return if possible
274 *
275  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
276  DO 10 j = 1, nrhs
277  ferr( j ) = zero
278  berr( j ) = zero
279  10 CONTINUE
280  RETURN
281  END IF
282 *
283  IF( notran ) THEN
284  transn = 'N'
285  transt = 'C'
286  ELSE
287  transn = 'C'
288  transt = 'N'
289  END IF
290 *
291 * NZ = maximum number of nonzero elements in each row of A, plus 1
292 *
293  nz = n + 1
294  eps = dlamch( 'Epsilon' )
295  safmin = dlamch( 'Safe minimum' )
296  safe1 = nz*safmin
297  safe2 = safe1 / eps
298 *
299 * Do for each right hand side
300 *
301  DO 140 j = 1, nrhs
302 *
303  count = 1
304  lstres = three
305  20 CONTINUE
306 *
307 * Loop until stopping criterion is satisfied.
308 *
309 * Compute residual R = B - op(A) * X,
310 * where op(A) = A, A**T, or A**H, depending on TRANS.
311 *
312  CALL zcopy( n, b( 1, j ), 1, work, 1 )
313  CALL zgemv( trans, n, n, -one, a, lda, x( 1, j ), 1, one, work,
314  $ 1 )
315 *
316 * Compute componentwise relative backward error from formula
317 *
318 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
319 *
320 * where abs(Z) is the componentwise absolute value of the matrix
321 * or vector Z. If the i-th component of the denominator is less
322 * than SAFE2, then SAFE1 is added to the i-th components of the
323 * numerator and denominator before dividing.
324 *
325  DO 30 i = 1, n
326  rwork( i ) = cabs1( b( i, j ) )
327  30 CONTINUE
328 *
329 * Compute abs(op(A))*abs(X) + abs(B).
330 *
331  IF( notran ) THEN
332  DO 50 k = 1, n
333  xk = cabs1( x( k, j ) )
334  DO 40 i = 1, n
335  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
336  40 CONTINUE
337  50 CONTINUE
338  ELSE
339  DO 70 k = 1, n
340  s = zero
341  DO 60 i = 1, n
342  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
343  60 CONTINUE
344  rwork( k ) = rwork( k ) + s
345  70 CONTINUE
346  END IF
347  s = zero
348  DO 80 i = 1, n
349  IF( rwork( i ).GT.safe2 ) THEN
350  s = max( s, cabs1( work( i ) ) / rwork( i ) )
351  ELSE
352  s = max( s, ( cabs1( work( i ) )+safe1 ) /
353  $ ( rwork( i )+safe1 ) )
354  END IF
355  80 CONTINUE
356  berr( j ) = s
357 *
358 * Test stopping criterion. Continue iterating if
359 * 1) The residual BERR(J) is larger than machine epsilon, and
360 * 2) BERR(J) decreased by at least a factor of 2 during the
361 * last iteration, and
362 * 3) At most ITMAX iterations tried.
363 *
364  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
365  $ count.LE.itmax ) THEN
366 *
367 * Update solution and try again.
368 *
369  CALL zgetrs( trans, n, 1, af, ldaf, ipiv, work, n, info )
370  CALL zaxpy( n, one, work, 1, x( 1, j ), 1 )
371  lstres = berr( j )
372  count = count + 1
373  GO TO 20
374  END IF
375 *
376 * Bound error from formula
377 *
378 * norm(X - XTRUE) / norm(X) .le. FERR =
379 * norm( abs(inv(op(A)))*
380 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
381 *
382 * where
383 * norm(Z) is the magnitude of the largest component of Z
384 * inv(op(A)) is the inverse of op(A)
385 * abs(Z) is the componentwise absolute value of the matrix or
386 * vector Z
387 * NZ is the maximum number of nonzeros in any row of A, plus 1
388 * EPS is machine epsilon
389 *
390 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
391 * is incremented by SAFE1 if the i-th component of
392 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
393 *
394 * Use ZLACN2 to estimate the infinity-norm of the matrix
395 * inv(op(A)) * diag(W),
396 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
397 *
398  DO 90 i = 1, n
399  IF( rwork( i ).GT.safe2 ) THEN
400  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
401  ELSE
402  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
403  $ safe1
404  END IF
405  90 CONTINUE
406 *
407  kase = 0
408  100 CONTINUE
409  CALL zlacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
410  IF( kase.NE.0 ) THEN
411  IF( kase.EQ.1 ) THEN
412 *
413 * Multiply by diag(W)*inv(op(A)**H).
414 *
415  CALL zgetrs( transt, n, 1, af, ldaf, ipiv, work, n,
416  $ info )
417  DO 110 i = 1, n
418  work( i ) = rwork( i )*work( i )
419  110 CONTINUE
420  ELSE
421 *
422 * Multiply by inv(op(A))*diag(W).
423 *
424  DO 120 i = 1, n
425  work( i ) = rwork( i )*work( i )
426  120 CONTINUE
427  CALL zgetrs( transn, n, 1, af, ldaf, ipiv, work, n,
428  $ info )
429  END IF
430  GO TO 100
431  END IF
432 *
433 * Normalize error.
434 *
435  lstres = zero
436  DO 130 i = 1, n
437  lstres = max( lstres, cabs1( x( i, j ) ) )
438  130 CONTINUE
439  IF( lstres.NE.zero )
440  $ ferr( j ) = ferr( j ) / lstres
441 *
442  140 CONTINUE
443 *
444  RETURN
445 *
446 * End of ZGERFS
447 *
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
Definition: zgetrs.f:123

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgerfsx ( character  TRANS,
character  EQUED,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldaf, * )  AF,
integer  LDAF,
integer, dimension( * )  IPIV,
double precision, dimension( * )  R,
double precision, dimension( * )  C,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx , * )  X,
integer  LDX,
double precision  RCOND,
double precision, dimension( * )  BERR,
integer  N_ERR_BNDS,
double precision, dimension( nrhs, * )  ERR_BNDS_NORM,
double precision, dimension( nrhs, * )  ERR_BNDS_COMP,
integer  NPARAMS,
double precision, dimension( * )  PARAMS,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGERFSX

Download ZGERFSX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
    ZGERFSX improves the computed solution to a system of linear
    equations and provides error bounds and backward error estimates
    for the solution.  In addition to normwise error bound, the code
    provides maximum componentwise error bound if possible.  See
    comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
    error bounds.

    The original system of linear equations may have been equilibrated
    before calling this routine, as described by arguments EQUED, R
    and C below. In this case, the solution and error bounds returned
    are for the original unequilibrated system.
     Some optional parameters are bundled in the PARAMS array.  These
     settings determine how refinement is performed, but often the
     defaults are acceptable.  If the defaults are acceptable, users
     can pass NPARAMS = 0 which prevents the source code from accessing
     the PARAMS argument.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
     Specifies the form of the system of equations:
       = 'N':  A * X = B     (No transpose)
       = 'T':  A**T * X = B  (Transpose)
       = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
[in]EQUED
          EQUED is CHARACTER*1
     Specifies the form of equilibration that was done to A
     before calling this routine. This is needed to compute
     the solution and error bounds correctly.
       = 'N':  No equilibration
       = 'R':  Row equilibration, i.e., A has been premultiplied by
               diag(R).
       = 'C':  Column equilibration, i.e., A has been postmultiplied
               by diag(C).
       = 'B':  Both row and column equilibration, i.e., A has been
               replaced by diag(R) * A * diag(C).
               The right hand side B has been changed accordingly.
[in]N
          N is INTEGER
     The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
     The number of right hand sides, i.e., the number of columns
     of the matrices B and X.  NRHS >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
     The original N-by-N matrix A.
[in]LDA
          LDA is INTEGER
     The leading dimension of the array A.  LDA >= max(1,N).
[in]AF
          AF is COMPLEX*16 array, dimension (LDAF,N)
     The factors L and U from the factorization A = P*L*U
     as computed by ZGETRF.
[in]LDAF
          LDAF is INTEGER
     The leading dimension of the array AF.  LDAF >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
     The pivot indices from ZGETRF; for 1<=i<=N, row i of the
     matrix was interchanged with row IPIV(i).
[in]R
          R is DOUBLE PRECISION array, dimension (N)
     The row scale factors for A.  If EQUED = 'R' or 'B', A is
     multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
     is not accessed.  
     If R is accessed, each element of R should be a power of the radix
     to ensure a reliable solution and error estimates. Scaling by
     powers of the radix does not cause rounding errors unless the
     result underflows or overflows. Rounding errors during scaling
     lead to refining with a matrix that is not equivalent to the
     input matrix, producing error estimates that may not be
     reliable.
[in]C
          C is DOUBLE PRECISION array, dimension (N)
     The column scale factors for A.  If EQUED = 'C' or 'B', A is
     multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
     is not accessed.
     If C is accessed, each element of C should be a power of the radix
     to ensure a reliable solution and error estimates. Scaling by
     powers of the radix does not cause rounding errors unless the
     result underflows or overflows. Rounding errors during scaling
     lead to refining with a matrix that is not equivalent to the
     input matrix, producing error estimates that may not be
     reliable.
[in]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
     The right hand side matrix B.
[in]LDB
          LDB is INTEGER
     The leading dimension of the array B.  LDB >= max(1,N).
[in,out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
     On entry, the solution matrix X, as computed by ZGETRS.
     On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
     The leading dimension of the array X.  LDX >= max(1,N).
[out]RCOND
          RCOND is DOUBLE PRECISION
     Reciprocal scaled condition number.  This is an estimate of the
     reciprocal Skeel condition number of the matrix A after
     equilibration (if done).  If this is less than the machine
     precision (in particular, if it is zero), the matrix is singular
     to working precision.  Note that the error may still be small even
     if this number is very small and the matrix appears ill-
     conditioned.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
     Componentwise relative backward error.  This is the
     componentwise relative backward error of each solution vector X(j)
     (i.e., the smallest relative change in any element of A or B that
     makes X(j) an exact solution).
[in]N_ERR_BNDS
          N_ERR_BNDS is INTEGER
     Number of error bounds to return for each right hand side
     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
     ERR_BNDS_COMP below.
[out]ERR_BNDS_NORM
          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     normwise relative error, which is defined as follows:

     Normwise relative error in the ith solution vector:
             max_j (abs(XTRUE(j,i) - X(j,i)))
            ------------------------------
                  max_j abs(X(j,i))

     The array is indexed by the type of error information as described
     below. There currently are up to three pieces of information
     returned.

     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_NORM(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated normwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*A, where S scales each row by a power of the
              radix so all absolute row sums of Z are approximately 1.

     See Lapack Working Note 165 for further details and extra
     cautions.
[out]ERR_BNDS_COMP
          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     componentwise relative error, which is defined as follows:

     Componentwise relative error in the ith solution vector:
                    abs(XTRUE(j,i) - X(j,i))
             max_j ----------------------
                         abs(X(j,i))

     The array is indexed by the right-hand side i (on which the
     componentwise relative error depends), and the type of error
     information as described below. There currently are up to three
     pieces of information returned for each right-hand side. If
     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most
     the first (:,N_ERR_BNDS) entries are returned.

     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_COMP(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated componentwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*(A*diag(x)), where x is the solution for the
              current right-hand side and S scales each row of
              A*diag(x) by a power of the radix so all absolute row
              sums of Z are approximately 1.

     See Lapack Working Note 165 for further details and extra
     cautions.
[in]NPARAMS
          NPARAMS is INTEGER
     Specifies the number of parameters set in PARAMS.  If .LE. 0, the
     PARAMS array is never referenced and default values are used.
[in,out]PARAMS
          PARAMS is DOUBLE PRECISION array, dimension NPARAMS
     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
     that entry will be filled with default value used for that
     parameter.  Only positions up to NPARAMS are accessed; defaults
     are used for higher-numbered parameters.

       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
            refinement or not.
         Default: 1.0D+0
            = 0.0 : No refinement is performed, and no error bounds are
                    computed.
            = 1.0 : Use the double-precision refinement algorithm,
                    possibly with doubled-single computations if the
                    compilation environment does not support DOUBLE
                    PRECISION.
              (other values are reserved for future use)

       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
            computations allowed for refinement.
         Default: 10
         Aggressive: Set to 100 to permit convergence using approximate
                     factorizations or factorizations other than LU. If
                     the factorization uses a technique other than
                     Gaussian elimination, the guarantees in
                     err_bnds_norm and err_bnds_comp may no longer be
                     trustworthy.

       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
            will attempt to find a solution with small componentwise
            relative error in the double-precision algorithm.  Positive
            is true, 0.0 is false.
         Default: 1.0 (attempt componentwise convergence)
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (2*N)
[out]INFO
          INFO is INTEGER
       = 0:  Successful exit. The solution to every right-hand side is
         guaranteed.
       < 0:  If INFO = -i, the i-th argument had an illegal value
       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
         has been completed, but the factor U is exactly singular, so
         the solution and error bounds could not be computed. RCOND = 0
         is returned.
       = N+J: The solution corresponding to the Jth right-hand side is
         not guaranteed. The solutions corresponding to other right-
         hand sides K with K > J may not be guaranteed as well, but
         only the first such right-hand side is reported. If a small
         componentwise error is not requested (PARAMS(3) = 0.0) then
         the Jth right-hand side is the first with a normwise error
         bound that is not guaranteed (the smallest J such
         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
         the Jth right-hand side is the first with either a normwise or
         componentwise error bound that is not guaranteed (the smallest
         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
         about all of the right-hand sides check ERR_BNDS_NORM or
         ERR_BNDS_COMP.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 416 of file zgerfsx.f.

416 *
417 * -- LAPACK computational routine (version 3.4.0) --
418 * -- LAPACK is a software package provided by Univ. of Tennessee, --
419 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
420 * November 2011
421 *
422 * .. Scalar Arguments ..
423  CHARACTER trans, equed
424  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs, nparams,
425  $ n_err_bnds
426  DOUBLE PRECISION rcond
427 * ..
428 * .. Array Arguments ..
429  INTEGER ipiv( * )
430  COMPLEX*16 a( lda, * ), af( ldaf, * ), b( ldb, * ),
431  $ x( ldx , * ), work( * )
432  DOUBLE PRECISION r( * ), c( * ), params( * ), berr( * ),
433  $ err_bnds_norm( nrhs, * ),
434  $ err_bnds_comp( nrhs, * ), rwork( * )
435 * ..
436 *
437 * ==================================================================
438 *
439 * .. Parameters ..
440  DOUBLE PRECISION zero, one
441  parameter( zero = 0.0d+0, one = 1.0d+0 )
442  DOUBLE PRECISION itref_default, ithresh_default
443  DOUBLE PRECISION componentwise_default, rthresh_default
444  DOUBLE PRECISION dzthresh_default
445  parameter( itref_default = 1.0d+0 )
446  parameter( ithresh_default = 10.0d+0 )
447  parameter( componentwise_default = 1.0d+0 )
448  parameter( rthresh_default = 0.5d+0 )
449  parameter( dzthresh_default = 0.25d+0 )
450  INTEGER la_linrx_itref_i, la_linrx_ithresh_i,
451  $ la_linrx_cwise_i
452  parameter( la_linrx_itref_i = 1,
453  $ la_linrx_ithresh_i = 2 )
454  parameter( la_linrx_cwise_i = 3 )
455  INTEGER la_linrx_trust_i, la_linrx_err_i,
456  $ la_linrx_rcond_i
457  parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
458  parameter( la_linrx_rcond_i = 3 )
459 * ..
460 * .. Local Scalars ..
461  CHARACTER(1) norm
462  LOGICAL rowequ, colequ, notran
463  INTEGER j, trans_type, prec_type, ref_type
464  INTEGER n_norms
465  DOUBLE PRECISION anorm, rcond_tmp
466  DOUBLE PRECISION illrcond_thresh, err_lbnd, cwise_wrong
467  LOGICAL ignore_cwise
468  INTEGER ithresh
469  DOUBLE PRECISION rthresh, unstable_thresh
470 * ..
471 * .. External Subroutines ..
473 * ..
474 * .. Intrinsic Functions ..
475  INTRINSIC max, sqrt, transfer
476 * ..
477 * .. External Functions ..
478  EXTERNAL lsame, blas_fpinfo_x, ilatrans, ilaprec
480  DOUBLE PRECISION dlamch, zlange, zla_gercond_x, zla_gercond_c
481  LOGICAL lsame
482  INTEGER blas_fpinfo_x
483  INTEGER ilatrans, ilaprec
484 * ..
485 * .. Executable Statements ..
486 *
487 * Check the input parameters.
488 *
489  info = 0
490  trans_type = ilatrans( trans )
491  ref_type = int( itref_default )
492  IF ( nparams .GE. la_linrx_itref_i ) THEN
493  IF ( params( la_linrx_itref_i ) .LT. 0.0d+0 ) THEN
494  params( la_linrx_itref_i ) = itref_default
495  ELSE
496  ref_type = params( la_linrx_itref_i )
497  END IF
498  END IF
499 *
500 * Set default parameters.
501 *
502  illrcond_thresh = dble( n ) * dlamch( 'Epsilon' )
503  ithresh = int( ithresh_default )
504  rthresh = rthresh_default
505  unstable_thresh = dzthresh_default
506  ignore_cwise = componentwise_default .EQ. 0.0d+0
507 *
508  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
509  IF ( params( la_linrx_ithresh_i ).LT.0.0d+0 ) THEN
510  params(la_linrx_ithresh_i) = ithresh
511  ELSE
512  ithresh = int( params( la_linrx_ithresh_i ) )
513  END IF
514  END IF
515  IF ( nparams.GE.la_linrx_cwise_i ) THEN
516  IF ( params( la_linrx_cwise_i ).LT.0.0d+0 ) THEN
517  IF ( ignore_cwise ) THEN
518  params( la_linrx_cwise_i ) = 0.0d+0
519  ELSE
520  params( la_linrx_cwise_i ) = 1.0d+0
521  END IF
522  ELSE
523  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0d+0
524  END IF
525  END IF
526  IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
527  n_norms = 0
528  ELSE IF ( ignore_cwise ) THEN
529  n_norms = 1
530  ELSE
531  n_norms = 2
532  END IF
533 *
534  notran = lsame( trans, 'N' )
535  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
536  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
537 *
538 * Test input parameters.
539 *
540  IF( trans_type.EQ.-1 ) THEN
541  info = -1
542  ELSE IF( .NOT.rowequ .AND. .NOT.colequ .AND.
543  $ .NOT.lsame( equed, 'N' ) ) THEN
544  info = -2
545  ELSE IF( n.LT.0 ) THEN
546  info = -3
547  ELSE IF( nrhs.LT.0 ) THEN
548  info = -4
549  ELSE IF( lda.LT.max( 1, n ) ) THEN
550  info = -6
551  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
552  info = -8
553  ELSE IF( ldb.LT.max( 1, n ) ) THEN
554  info = -13
555  ELSE IF( ldx.LT.max( 1, n ) ) THEN
556  info = -15
557  END IF
558  IF( info.NE.0 ) THEN
559  CALL xerbla( 'ZGERFSX', -info )
560  RETURN
561  END IF
562 *
563 * Quick return if possible.
564 *
565  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
566  rcond = 1.0d+0
567  DO j = 1, nrhs
568  berr( j ) = 0.0d+0
569  IF ( n_err_bnds .GE. 1 ) THEN
570  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
571  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
572  END IF
573  IF ( n_err_bnds .GE. 2 ) THEN
574  err_bnds_norm( j, la_linrx_err_i ) = 0.0d+0
575  err_bnds_comp( j, la_linrx_err_i ) = 0.0d+0
576  END IF
577  IF ( n_err_bnds .GE. 3 ) THEN
578  err_bnds_norm( j, la_linrx_rcond_i ) = 1.0d+0
579  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0d+0
580  END IF
581  END DO
582  RETURN
583  END IF
584 *
585 * Default to failure.
586 *
587  rcond = 0.0d+0
588  DO j = 1, nrhs
589  berr( j ) = 1.0d+0
590  IF ( n_err_bnds .GE. 1 ) THEN
591  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
592  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
593  END IF
594  IF ( n_err_bnds .GE. 2 ) THEN
595  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
596  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
597  END IF
598  IF ( n_err_bnds .GE. 3 ) THEN
599  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0d+0
600  err_bnds_comp( j, la_linrx_rcond_i ) = 0.0d+0
601  END IF
602  END DO
603 *
604 * Compute the norm of A and the reciprocal of the condition
605 * number of A.
606 *
607  IF( notran ) THEN
608  norm = 'I'
609  ELSE
610  norm = '1'
611  END IF
612  anorm = zlange( norm, n, n, a, lda, rwork )
613  CALL zgecon( norm, n, af, ldaf, anorm, rcond, work, rwork, info )
614 *
615 * Perform refinement on each right-hand side
616 *
617  IF ( ref_type .NE. 0 ) THEN
618 
619  prec_type = ilaprec( 'E' )
620 
621  IF ( notran ) THEN
622  CALL zla_gerfsx_extended( prec_type, trans_type, n,
623  $ nrhs, a, lda, af, ldaf, ipiv, colequ, c, b,
624  $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
625  $ err_bnds_comp, work, rwork, work(n+1),
626  $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
627  $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
628  $ info )
629  ELSE
630  CALL zla_gerfsx_extended( prec_type, trans_type, n,
631  $ nrhs, a, lda, af, ldaf, ipiv, rowequ, r, b,
632  $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
633  $ err_bnds_comp, work, rwork, work(n+1),
634  $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
635  $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
636  $ info )
637  END IF
638  END IF
639 
640  err_lbnd = max( 10.0d+0, sqrt( dble( n ) ) ) * dlamch( 'Epsilon' )
641  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
642 *
643 * Compute scaled normwise condition number cond(A*C).
644 *
645  IF ( colequ .AND. notran ) THEN
646  rcond_tmp = zla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
647  $ c, .true., info, work, rwork )
648  ELSE IF ( rowequ .AND. .NOT. notran ) THEN
649  rcond_tmp = zla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
650  $ r, .true., info, work, rwork )
651  ELSE
652  rcond_tmp = zla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
653  $ c, .false., info, work, rwork )
654  END IF
655  DO j = 1, nrhs
656 *
657 * Cap the error at 1.0.
658 *
659  IF ( n_err_bnds .GE. la_linrx_err_i
660  $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0d+0 )
661  $ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
662 *
663 * Threshold the error (see LAWN).
664 *
665  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
666  err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
667  err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+0
668  IF ( info .LE. n ) info = n + j
669  ELSE IF (err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd)
670  $ THEN
671  err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
672  err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
673  END IF
674 *
675 * Save the condition number.
676 *
677  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
678  err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
679  END IF
680  END DO
681  END IF
682 
683  IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 2 ) THEN
684 *
685 * Compute componentwise condition number cond(A*diag(Y(:,J))) for
686 * each right-hand side using the current solution as an estimate of
687 * the true solution. If the componentwise error estimate is too
688 * large, then the solution is a lousy estimate of truth and the
689 * estimated RCOND may be too optimistic. To avoid misleading users,
690 * the inverse condition number is set to 0.0 when the estimated
691 * cwise error is at least CWISE_WRONG.
692 *
693  cwise_wrong = sqrt( dlamch( 'Epsilon' ) )
694  DO j = 1, nrhs
695  IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
696  $ THEN
697  rcond_tmp = zla_gercond_x( trans, n, a, lda, af, ldaf,
698  $ ipiv, x(1,j), info, work, rwork )
699  ELSE
700  rcond_tmp = 0.0d+0
701  END IF
702 *
703 * Cap the error at 1.0.
704 *
705  IF ( n_err_bnds .GE. la_linrx_err_i
706  $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0d+0 )
707  $ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
708 *
709 * Threshold the error (see LAWN).
710 *
711  IF ( rcond_tmp .LT. illrcond_thresh ) THEN
712  err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
713  err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
714  IF ( params( la_linrx_cwise_i ) .EQ. 1.0d+0
715  $ .AND. info.LT.n + j ) info = n + j
716  ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
717  $ .LT. err_lbnd ) THEN
718  err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
719  err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
720  END IF
721 *
722 * Save the condition number.
723 *
724  IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
725  err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
726  END IF
727 
728  END DO
729  END IF
730 *
731  RETURN
732 *
733 * End of ZGERFSX
734 *
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:60
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function zla_gercond_x(TRANS, N, A, LDA, AF, LDAF, IPIV, X, INFO, WORK, RWORK)
ZLA_GERCOND_X computes the infinity norm condition number of op(A)*diag(x) for general matrices...
double precision function zla_gercond_c(TRANS, N, A, LDA, AF, LDAF, IPIV, C, CAPPLY, INFO, WORK, RWORK)
ZLA_GERCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general matrices...
subroutine zgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
ZGECON
Definition: zgecon.f:126
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:60
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:117
subroutine zla_gerfsx_extended(PREC_TYPE, TRANS_TYPE, N, NRHS, A, LDA, AF, LDAF, IPIV, COLEQU, C, B, LDB, Y, LDY, BERR_OUT, N_NORMS, ERRS_N, ERRS_C, RES, AYB, DY, Y_TAIL, RCOND, ITHRESH, RTHRESH, DZ_UB, IGNORE_CWISE, INFO)
ZLA_GERFSX_EXTENDED

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgerq2 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  INFO 
)

ZGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm.

Download ZGERQ2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGERQ2 computes an RQ factorization of a complex m by n matrix A:
 A = R * Q.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the m by n matrix A.
          On exit, if m <= n, the upper triangle of the subarray
          A(1:m,n-m+1:n) contains the m by m upper triangular matrix R;
          if m >= n, the elements on and above the (m-n)-th subdiagonal
          contain the m by n upper trapezoidal matrix R; the remaining
          elements, with the array TAU, represent the unitary matrix
          Q as a product of elementary reflectors (see Further
          Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (M)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(1)**H H(2)**H . . . H(k)**H, where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(n-k+i+1:n) = 0 and v(n-k+i) = 1; conjg(v(1:n-k+i-1)) is stored on
  exit in A(m-k+i,1:n-k+i-1), and tau in TAU(i).

Definition at line 125 of file zgerq2.f.

125 *
126 * -- LAPACK computational routine (version 3.4.2) --
127 * -- LAPACK is a software package provided by Univ. of Tennessee, --
128 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
129 * September 2012
130 *
131 * .. Scalar Arguments ..
132  INTEGER info, lda, m, n
133 * ..
134 * .. Array Arguments ..
135  COMPLEX*16 a( lda, * ), tau( * ), work( * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  COMPLEX*16 one
142  parameter( one = ( 1.0d+0, 0.0d+0 ) )
143 * ..
144 * .. Local Scalars ..
145  INTEGER i, k
146  COMPLEX*16 alpha
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL xerbla, zlacgv, zlarf, zlarfg
150 * ..
151 * .. Intrinsic Functions ..
152  INTRINSIC max, min
153 * ..
154 * .. Executable Statements ..
155 *
156 * Test the input arguments
157 *
158  info = 0
159  IF( m.LT.0 ) THEN
160  info = -1
161  ELSE IF( n.LT.0 ) THEN
162  info = -2
163  ELSE IF( lda.LT.max( 1, m ) ) THEN
164  info = -4
165  END IF
166  IF( info.NE.0 ) THEN
167  CALL xerbla( 'ZGERQ2', -info )
168  RETURN
169  END IF
170 *
171  k = min( m, n )
172 *
173  DO 10 i = k, 1, -1
174 *
175 * Generate elementary reflector H(i) to annihilate
176 * A(m-k+i,1:n-k+i-1)
177 *
178  CALL zlacgv( n-k+i, a( m-k+i, 1 ), lda )
179  alpha = a( m-k+i, n-k+i )
180  CALL zlarfg( n-k+i, alpha, a( m-k+i, 1 ), lda, tau( i ) )
181 *
182 * Apply H(i) to A(1:m-k+i-1,1:n-k+i) from the right
183 *
184  a( m-k+i, n-k+i ) = one
185  CALL zlarf( 'Right', m-k+i-1, n-k+i, a( m-k+i, 1 ), lda,
186  $ tau( i ), a, lda, work )
187  a( m-k+i, n-k+i ) = alpha
188  CALL zlacgv( n-k+i-1, a( m-k+i, 1 ), lda )
189  10 CONTINUE
190  RETURN
191 *
192 * End of ZGERQ2
193 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgerqf ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGERQF

Download ZGERQF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGERQF computes an RQ factorization of a complex M-by-N matrix A:
 A = R * Q.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit,
          if m <= n, the upper triangle of the subarray
          A(1:m,n-m+1:n) contains the M-by-M upper triangular matrix R;
          if m >= n, the elements on and above the (m-n)-th subdiagonal
          contain the M-by-N upper trapezoidal matrix R;
          the remaining elements, with the array TAU, represent the
          unitary matrix Q as a product of min(m,n) elementary
          reflectors (see Further Details).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]TAU
          TAU is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,M).
          For optimum performance LWORK >= M*NB, where NB is
          the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  The matrix Q is represented as a product of elementary reflectors

     Q = H(1)**H H(2)**H . . . H(k)**H, where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**H

  where tau is a complex scalar, and v is a complex vector with
  v(n-k+i+1:n) = 0 and v(n-k+i) = 1; conjg(v(1:n-k+i-1)) is stored on
  exit in A(m-k+i,1:n-k+i-1), and tau in TAU(i).

Definition at line 140 of file zgerqf.f.

140 *
141 * -- LAPACK computational routine (version 3.4.0) --
142 * -- LAPACK is a software package provided by Univ. of Tennessee, --
143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144 * November 2011
145 *
146 * .. Scalar Arguments ..
147  INTEGER info, lda, lwork, m, n
148 * ..
149 * .. Array Arguments ..
150  COMPLEX*16 a( lda, * ), tau( * ), work( * )
151 * ..
152 *
153 * =====================================================================
154 *
155 * .. Local Scalars ..
156  LOGICAL lquery
157  INTEGER i, ib, iinfo, iws, k, ki, kk, ldwork, lwkopt,
158  $ mu, nb, nbmin, nu, nx
159 * ..
160 * .. External Subroutines ..
161  EXTERNAL xerbla, zgerq2, zlarfb, zlarft
162 * ..
163 * .. Intrinsic Functions ..
164  INTRINSIC max, min
165 * ..
166 * .. External Functions ..
167  INTEGER ilaenv
168  EXTERNAL ilaenv
169 * ..
170 * .. Executable Statements ..
171 *
172 * Test the input arguments
173 *
174  info = 0
175  lquery = ( lwork.EQ.-1 )
176  IF( m.LT.0 ) THEN
177  info = -1
178  ELSE IF( n.LT.0 ) THEN
179  info = -2
180  ELSE IF( lda.LT.max( 1, m ) ) THEN
181  info = -4
182  END IF
183 *
184  IF( info.EQ.0 ) THEN
185  k = min( m, n )
186  IF( k.EQ.0 ) THEN
187  lwkopt = 1
188  ELSE
189  nb = ilaenv( 1, 'ZGERQF', ' ', m, n, -1, -1 )
190  lwkopt = m*nb
191  END IF
192  work( 1 ) = lwkopt
193 *
194  IF( lwork.LT.max( 1, m ) .AND. .NOT.lquery ) THEN
195  info = -7
196  END IF
197  END IF
198 *
199  IF( info.NE.0 ) THEN
200  CALL xerbla( 'ZGERQF', -info )
201  RETURN
202  ELSE IF( lquery ) THEN
203  RETURN
204  END IF
205 *
206 * Quick return if possible
207 *
208  IF( k.EQ.0 ) THEN
209  RETURN
210  END IF
211 *
212  nbmin = 2
213  nx = 1
214  iws = m
215  IF( nb.GT.1 .AND. nb.LT.k ) THEN
216 *
217 * Determine when to cross over from blocked to unblocked code.
218 *
219  nx = max( 0, ilaenv( 3, 'ZGERQF', ' ', m, n, -1, -1 ) )
220  IF( nx.LT.k ) THEN
221 *
222 * Determine if workspace is large enough for blocked code.
223 *
224  ldwork = m
225  iws = ldwork*nb
226  IF( lwork.LT.iws ) THEN
227 *
228 * Not enough workspace to use optimal NB: reduce NB and
229 * determine the minimum value of NB.
230 *
231  nb = lwork / ldwork
232  nbmin = max( 2, ilaenv( 2, 'ZGERQF', ' ', m, n, -1,
233  $ -1 ) )
234  END IF
235  END IF
236  END IF
237 *
238  IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
239 *
240 * Use blocked code initially.
241 * The last kk rows are handled by the block method.
242 *
243  ki = ( ( k-nx-1 ) / nb )*nb
244  kk = min( k, ki+nb )
245 *
246  DO 10 i = k - kk + ki + 1, k - kk + 1, -nb
247  ib = min( k-i+1, nb )
248 *
249 * Compute the RQ factorization of the current block
250 * A(m-k+i:m-k+i+ib-1,1:n-k+i+ib-1)
251 *
252  CALL zgerq2( ib, n-k+i+ib-1, a( m-k+i, 1 ), lda, tau( i ),
253  $ work, iinfo )
254  IF( m-k+i.GT.1 ) THEN
255 *
256 * Form the triangular factor of the block reflector
257 * H = H(i+ib-1) . . . H(i+1) H(i)
258 *
259  CALL zlarft( 'Backward', 'Rowwise', n-k+i+ib-1, ib,
260  $ a( m-k+i, 1 ), lda, tau( i ), work, ldwork )
261 *
262 * Apply H to A(1:m-k+i-1,1:n-k+i+ib-1) from the right
263 *
264  CALL zlarfb( 'Right', 'No transpose', 'Backward',
265  $ 'Rowwise', m-k+i-1, n-k+i+ib-1, ib,
266  $ a( m-k+i, 1 ), lda, work, ldwork, a, lda,
267  $ work( ib+1 ), ldwork )
268  END IF
269  10 CONTINUE
270  mu = m - k + i + nb - 1
271  nu = n - k + i + nb - 1
272  ELSE
273  mu = m
274  nu = n
275  END IF
276 *
277 * Use unblocked code to factor the last or only block
278 *
279  IF( mu.GT.0 .AND. nu.GT.0 )
280  $ CALL zgerq2( mu, nu, a, lda, tau, work, iinfo )
281 *
282  work( 1 ) = iws
283  RETURN
284 *
285 * End of ZGERQF
286 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: zlarfb.f:197
subroutine zlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
ZLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: zlarft.f:165
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zgerq2(M, N, A, LDA, TAU, WORK, INFO)
ZGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: zgerq2.f:125

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgetf2 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

ZGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row interchanges (unblocked algorithm).

Download ZGETF2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGETF2 computes an LU factorization of a general m-by-n matrix A
 using partial pivoting with row interchanges.

 The factorization has the form
    A = P * L * U
 where P is a permutation matrix, L is lower triangular with unit
 diagonal elements (lower trapezoidal if m > n), and U is upper
 triangular (upper trapezoidal if m < n).

 This is the right-looking Level 2 BLAS version of the algorithm.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the m by n matrix to be factored.
          On exit, the factors L and U from the factorization
          A = P*L*U; the unit diagonal elements of L are not stored.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]IPIV
          IPIV is INTEGER array, dimension (min(M,N))
          The pivot indices; for 1 <= i <= min(M,N), row i of the
          matrix was interchanged with row IPIV(i).
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -k, the k-th argument had an illegal value
          > 0: if INFO = k, U(k,k) is exactly zero. The factorization
               has been completed, but the factor U is exactly
               singular, and division by zero will occur if it is used
               to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 110 of file zgetf2.f.

110 *
111 * -- LAPACK computational routine (version 3.4.2) --
112 * -- LAPACK is a software package provided by Univ. of Tennessee, --
113 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
114 * September 2012
115 *
116 * .. Scalar Arguments ..
117  INTEGER info, lda, m, n
118 * ..
119 * .. Array Arguments ..
120  INTEGER ipiv( * )
121  COMPLEX*16 a( lda, * )
122 * ..
123 *
124 * =====================================================================
125 *
126 * .. Parameters ..
127  COMPLEX*16 one, zero
128  parameter( one = ( 1.0d+0, 0.0d+0 ),
129  $ zero = ( 0.0d+0, 0.0d+0 ) )
130 * ..
131 * .. Local Scalars ..
132  DOUBLE PRECISION sfmin
133  INTEGER i, j, jp
134 * ..
135 * .. External Functions ..
136  DOUBLE PRECISION dlamch
137  INTEGER izamax
138  EXTERNAL dlamch, izamax
139 * ..
140 * .. External Subroutines ..
141  EXTERNAL xerbla, zgeru, zscal, zswap
142 * ..
143 * .. Intrinsic Functions ..
144  INTRINSIC max, min
145 * ..
146 * .. Executable Statements ..
147 *
148 * Test the input parameters.
149 *
150  info = 0
151  IF( m.LT.0 ) THEN
152  info = -1
153  ELSE IF( n.LT.0 ) THEN
154  info = -2
155  ELSE IF( lda.LT.max( 1, m ) ) THEN
156  info = -4
157  END IF
158  IF( info.NE.0 ) THEN
159  CALL xerbla( 'ZGETF2', -info )
160  RETURN
161  END IF
162 *
163 * Quick return if possible
164 *
165  IF( m.EQ.0 .OR. n.EQ.0 )
166  $ RETURN
167 *
168 * Compute machine safe minimum
169 *
170  sfmin = dlamch('S')
171 *
172  DO 10 j = 1, min( m, n )
173 *
174 * Find pivot and test for singularity.
175 *
176  jp = j - 1 + izamax( m-j+1, a( j, j ), 1 )
177  ipiv( j ) = jp
178  IF( a( jp, j ).NE.zero ) THEN
179 *
180 * Apply the interchange to columns 1:N.
181 *
182  IF( jp.NE.j )
183  $ CALL zswap( n, a( j, 1 ), lda, a( jp, 1 ), lda )
184 *
185 * Compute elements J+1:M of J-th column.
186 *
187  IF( j.LT.m ) THEN
188  IF( abs(a( j, j )) .GE. sfmin ) THEN
189  CALL zscal( m-j, one / a( j, j ), a( j+1, j ), 1 )
190  ELSE
191  DO 20 i = 1, m-j
192  a( j+i, j ) = a( j+i, j ) / a( j, j )
193  20 CONTINUE
194  END IF
195  END IF
196 *
197  ELSE IF( info.EQ.0 ) THEN
198 *
199  info = j
200  END IF
201 *
202  IF( j.LT.min( m, n ) ) THEN
203 *
204 * Update trailing submatrix.
205 *
206  CALL zgeru( m-j, n-j, -one, a( j+1, j ), 1, a( j, j+1 ),
207  $ lda, a( j+1, j+1 ), lda )
208  END IF
209  10 CONTINUE
210  RETURN
211 *
212 * End of ZGETF2
213 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
subroutine zgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERU
Definition: zgeru.f:132
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgetrf ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

ZGETRF

Download ZGETRF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGETRF computes an LU factorization of a general M-by-N matrix A
 using partial pivoting with row interchanges.

 The factorization has the form
    A = P * L * U
 where P is a permutation matrix, L is lower triangular with unit
 diagonal elements (lower trapezoidal if m > n), and U is upper
 triangular (upper trapezoidal if m < n).

 This is the right-looking Level 3 BLAS version of the algorithm.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix to be factored.
          On exit, the factors L and U from the factorization
          A = P*L*U; the unit diagonal elements of L are not stored.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]IPIV
          IPIV is INTEGER array, dimension (min(M,N))
          The pivot indices; for 1 <= i <= min(M,N), row i of the
          matrix was interchanged with row IPIV(i).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
                has been completed, but the factor U is exactly
                singular, and division by zero will occur if it is used
                to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 110 of file zgetrf.f.

110 *
111 * -- LAPACK computational routine (version 3.6.0) --
112 * -- LAPACK is a software package provided by Univ. of Tennessee, --
113 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
114 * November 2015
115 *
116 * .. Scalar Arguments ..
117  INTEGER info, lda, m, n
118 * ..
119 * .. Array Arguments ..
120  INTEGER ipiv( * )
121  COMPLEX*16 a( lda, * )
122 * ..
123 *
124 * =====================================================================
125 *
126 * .. Parameters ..
127  COMPLEX*16 one
128  parameter( one = ( 1.0d+0, 0.0d+0 ) )
129 * ..
130 * .. Local Scalars ..
131  INTEGER i, iinfo, j, jb, nb
132 * ..
133 * .. External Subroutines ..
134  EXTERNAL xerbla, zgemm, zgetrf2, zlaswp, ztrsm
135 * ..
136 * .. External Functions ..
137  INTEGER ilaenv
138  EXTERNAL ilaenv
139 * ..
140 * .. Intrinsic Functions ..
141  INTRINSIC max, min
142 * ..
143 * .. Executable Statements ..
144 *
145 * Test the input parameters.
146 *
147  info = 0
148  IF( m.LT.0 ) THEN
149  info = -1
150  ELSE IF( n.LT.0 ) THEN
151  info = -2
152  ELSE IF( lda.LT.max( 1, m ) ) THEN
153  info = -4
154  END IF
155  IF( info.NE.0 ) THEN
156  CALL xerbla( 'ZGETRF', -info )
157  RETURN
158  END IF
159 *
160 * Quick return if possible
161 *
162  IF( m.EQ.0 .OR. n.EQ.0 )
163  $ RETURN
164 *
165 * Determine the block size for this environment.
166 *
167  nb = ilaenv( 1, 'ZGETRF', ' ', m, n, -1, -1 )
168  IF( nb.LE.1 .OR. nb.GE.min( m, n ) ) THEN
169 *
170 * Use unblocked code.
171 *
172  CALL zgetrf2( m, n, a, lda, ipiv, info )
173  ELSE
174 *
175 * Use blocked code.
176 *
177  DO 20 j = 1, min( m, n ), nb
178  jb = min( min( m, n )-j+1, nb )
179 *
180 * Factor diagonal and subdiagonal blocks and test for exact
181 * singularity.
182 *
183  CALL zgetrf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo )
184 *
185 * Adjust INFO and the pivot indices.
186 *
187  IF( info.EQ.0 .AND. iinfo.GT.0 )
188  $ info = iinfo + j - 1
189  DO 10 i = j, min( m, j+jb-1 )
190  ipiv( i ) = j - 1 + ipiv( i )
191  10 CONTINUE
192 *
193 * Apply interchanges to columns 1:J-1.
194 *
195  CALL zlaswp( j-1, a, lda, j, j+jb-1, ipiv, 1 )
196 *
197  IF( j+jb.LE.n ) THEN
198 *
199 * Apply interchanges to columns J+JB:N.
200 *
201  CALL zlaswp( n-j-jb+1, a( 1, j+jb ), lda, j, j+jb-1,
202  $ ipiv, 1 )
203 *
204 * Compute block row of U.
205 *
206  CALL ztrsm( 'Left', 'Lower', 'No transpose', 'Unit', jb,
207  $ n-j-jb+1, one, a( j, j ), lda, a( j, j+jb ),
208  $ lda )
209  IF( j+jb.LE.m ) THEN
210 *
211 * Update trailing submatrix.
212 *
213  CALL zgemm( 'No transpose', 'No transpose', m-j-jb+1,
214  $ n-j-jb+1, jb, -one, a( j+jb, j ), lda,
215  $ a( j, j+jb ), lda, one, a( j+jb, j+jb ),
216  $ lda )
217  END IF
218  END IF
219  20 CONTINUE
220  END IF
221  RETURN
222 *
223 * End of ZGETRF
224 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
subroutine zlaswp(N, A, LDA, K1, K2, IPIV, INCX)
ZLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: zlaswp.f:116
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
recursive subroutine zgetrf2(M, N, A, LDA, IPIV, INFO)
ZGETRF2
Definition: zgetrf2.f:115

Here is the call graph for this function:

recursive subroutine zgetrf2 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

ZGETRF2

Purpose:
 ZGETRF2 computes an LU factorization of a general M-by-N matrix A
 using partial pivoting with row interchanges.

 The factorization has the form
    A = P * L * U
 where P is a permutation matrix, L is lower triangular with unit
 diagonal elements (lower trapezoidal if m > n), and U is upper
 triangular (upper trapezoidal if m < n).

 This is the recursive version of the algorithm. It divides
 the matrix into four submatrices:
            
        [  A11 | A12  ]  where A11 is n1 by n1 and A22 is n2 by n2
    A = [ -----|----- ]  with n1 = min(m,n)
        [  A21 | A22  ]       n2 = n-n1
            
                                       [ A11 ]
 The subroutine calls itself to factor [ --- ],
                                       [ A12 ]
                 [ A12 ]
 do the swaps on [ --- ], solve A12, update A22,
                 [ A22 ]

 then calls itself to factor A22 and do the swaps on A21.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix to be factored.
          On exit, the factors L and U from the factorization
          A = P*L*U; the unit diagonal elements of L are not stored.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]IPIV
          IPIV is INTEGER array, dimension (min(M,N))
          The pivot indices; for 1 <= i <= min(M,N), row i of the
          matrix was interchanged with row IPIV(i).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
                has been completed, but the factor U is exactly
                singular, and division by zero will occur if it is used
                to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 115 of file zgetrf2.f.

115 *
116 * -- LAPACK computational routine (version 3.6.0) --
117 * -- LAPACK is a software package provided by Univ. of Tennessee, --
118 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
119 * November 2015
120 *
121 * .. Scalar Arguments ..
122  INTEGER info, lda, m, n
123 * ..
124 * .. Array Arguments ..
125  INTEGER ipiv( * )
126  COMPLEX*16 a( lda, * )
127 * ..
128 *
129 * =====================================================================
130 *
131 * .. Parameters ..
132  COMPLEX*16 one, zero
133  parameter( one = ( 1.0d+0, 0.0d+0 ),
134  $ zero = ( 0.0d+0, 0.0d+0 ) )
135 * ..
136 * .. Local Scalars ..
137  DOUBLE PRECISION sfmin
138  COMPLEX*16 temp
139  INTEGER i, iinfo, n1, n2
140 * ..
141 * .. External Functions ..
142  DOUBLE PRECISION dlamch
143  INTEGER izamax
144  EXTERNAL dlamch, izamax
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL zgemm, zscal, zlaswp, ztrsm, zerbla
148 * ..
149 * .. Intrinsic Functions ..
150  INTRINSIC max, min
151 * ..
152 * .. Executable Statements ..
153 *
154 * Test the input parameters
155 *
156  info = 0
157  IF( m.LT.0 ) THEN
158  info = -1
159  ELSE IF( n.LT.0 ) THEN
160  info = -2
161  ELSE IF( lda.LT.max( 1, m ) ) THEN
162  info = -4
163  END IF
164  IF( info.NE.0 ) THEN
165  CALL xerbla( 'ZGETRF2', -info )
166  RETURN
167  END IF
168 *
169 * Quick return if possible
170 *
171  IF( m.EQ.0 .OR. n.EQ.0 )
172  $ RETURN
173 
174  IF ( m.EQ.1 ) THEN
175 *
176 * Use unblocked code for one row case
177 * Just need to handle IPIV and INFO
178 *
179  ipiv( 1 ) = 1
180  IF ( a(1,1).EQ.zero )
181  $ info = 1
182 *
183  ELSE IF( n.EQ.1 ) THEN
184 *
185 * Use unblocked code for one column case
186 *
187 *
188 * Compute machine safe minimum
189 *
190  sfmin = dlamch('S')
191 *
192 * Find pivot and test for singularity
193 *
194  i = izamax( m, a( 1, 1 ), 1 )
195  ipiv( 1 ) = i
196  IF( a( i, 1 ).NE.zero ) THEN
197 *
198 * Apply the interchange
199 *
200  IF( i.NE.1 ) THEN
201  temp = a( 1, 1 )
202  a( 1, 1 ) = a( i, 1 )
203  a( i, 1 ) = temp
204  END IF
205 *
206 * Compute elements 2:M of the column
207 *
208  IF( abs(a( 1, 1 )) .GE. sfmin ) THEN
209  CALL zscal( m-1, one / a( 1, 1 ), a( 2, 1 ), 1 )
210  ELSE
211  DO 10 i = 1, m-1
212  a( 1+i, 1 ) = a( 1+i, 1 ) / a( 1, 1 )
213  10 CONTINUE
214  END IF
215 *
216  ELSE
217  info = 1
218  END IF
219 
220  ELSE
221 *
222 * Use recursive code
223 *
224  n1 = min( m, n ) / 2
225  n2 = n-n1
226 *
227 * [ A11 ]
228 * Factor [ --- ]
229 * [ A21 ]
230 *
231  CALL zgetrf2( m, n1, a, lda, ipiv, iinfo )
232 
233  IF ( info.EQ.0 .AND. iinfo.GT.0 )
234  $ info = iinfo
235 *
236 * [ A12 ]
237 * Apply interchanges to [ --- ]
238 * [ A22 ]
239 *
240  CALL zlaswp( n2, a( 1, n1+1 ), lda, 1, n1, ipiv, 1 )
241 *
242 * Solve A12
243 *
244  CALL ztrsm( 'L', 'L', 'N', 'U', n1, n2, one, a, lda,
245  $ a( 1, n1+1 ), lda )
246 *
247 * Update A22
248 *
249  CALL zgemm( 'N', 'N', m-n1, n2, n1, -one, a( n1+1, 1 ), lda,
250  $ a( 1, n1+1 ), lda, one, a( n1+1, n1+1 ), lda )
251 *
252 * Factor A22
253 *
254  CALL zgetrf2( m-n1, n2, a( n1+1, n1+1 ), lda, ipiv( n1+1 ),
255  $ iinfo )
256 *
257 * Adjust INFO and the pivot indices
258 *
259  IF ( info.EQ.0 .AND. iinfo.GT.0 )
260  $ info = iinfo + n1
261  DO 20 i = n1+1, min( m, n )
262  ipiv( i ) = ipiv( i ) + n1
263  20 CONTINUE
264 *
265 * Apply interchanges to A21
266 *
267  CALL zlaswp( n1, a( 1, 1 ), lda, n1+1, min( m, n), ipiv, 1 )
268 *
269  END IF
270  RETURN
271 *
272 * End of ZGETRF2
273 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
subroutine zlaswp(N, A, LDA, K1, K2, IPIV, INCX)
ZLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: zlaswp.f:116
recursive subroutine zgetrf2(M, N, A, LDA, IPIV, INFO)
ZGETRF2
Definition: zgetrf2.f:115

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgetri ( integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGETRI

Download ZGETRI + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGETRI computes the inverse of a matrix using the LU factorization
 computed by ZGETRF.

 This method inverts U and then computes inv(A) by solving the system
 inv(A)*L = inv(U) for inv(A).
Parameters
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the factors L and U from the factorization
          A = P*L*U as computed by ZGETRF.
          On exit, if INFO = 0, the inverse of the original matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices from ZGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,N).
          For optimal performance LWORK >= N*NB, where NB is
          the optimal blocksize returned by ILAENV.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is
                singular and its inverse could not be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 116 of file zgetri.f.

116 *
117 * -- LAPACK computational routine (version 3.4.0) --
118 * -- LAPACK is a software package provided by Univ. of Tennessee, --
119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120 * November 2011
121 *
122 * .. Scalar Arguments ..
123  INTEGER info, lda, lwork, n
124 * ..
125 * .. Array Arguments ..
126  INTEGER ipiv( * )
127  COMPLEX*16 a( lda, * ), work( * )
128 * ..
129 *
130 * =====================================================================
131 *
132 * .. Parameters ..
133  COMPLEX*16 zero, one
134  parameter( zero = ( 0.0d+0, 0.0d+0 ),
135  $ one = ( 1.0d+0, 0.0d+0 ) )
136 * ..
137 * .. Local Scalars ..
138  LOGICAL lquery
139  INTEGER i, iws, j, jb, jj, jp, ldwork, lwkopt, nb,
140  $ nbmin, nn
141 * ..
142 * .. External Functions ..
143  INTEGER ilaenv
144  EXTERNAL ilaenv
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL xerbla, zgemm, zgemv, zswap, ztrsm, ztrtri
148 * ..
149 * .. Intrinsic Functions ..
150  INTRINSIC max, min
151 * ..
152 * .. Executable Statements ..
153 *
154 * Test the input parameters.
155 *
156  info = 0
157  nb = ilaenv( 1, 'ZGETRI', ' ', n, -1, -1, -1 )
158  lwkopt = n*nb
159  work( 1 ) = lwkopt
160  lquery = ( lwork.EQ.-1 )
161  IF( n.LT.0 ) THEN
162  info = -1
163  ELSE IF( lda.LT.max( 1, n ) ) THEN
164  info = -3
165  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
166  info = -6
167  END IF
168  IF( info.NE.0 ) THEN
169  CALL xerbla( 'ZGETRI', -info )
170  RETURN
171  ELSE IF( lquery ) THEN
172  RETURN
173  END IF
174 *
175 * Quick return if possible
176 *
177  IF( n.EQ.0 )
178  $ RETURN
179 *
180 * Form inv(U). If INFO > 0 from ZTRTRI, then U is singular,
181 * and the inverse is not computed.
182 *
183  CALL ztrtri( 'Upper', 'Non-unit', n, a, lda, info )
184  IF( info.GT.0 )
185  $ RETURN
186 *
187  nbmin = 2
188  ldwork = n
189  IF( nb.GT.1 .AND. nb.LT.n ) THEN
190  iws = max( ldwork*nb, 1 )
191  IF( lwork.LT.iws ) THEN
192  nb = lwork / ldwork
193  nbmin = max( 2, ilaenv( 2, 'ZGETRI', ' ', n, -1, -1, -1 ) )
194  END IF
195  ELSE
196  iws = n
197  END IF
198 *
199 * Solve the equation inv(A)*L = inv(U) for inv(A).
200 *
201  IF( nb.LT.nbmin .OR. nb.GE.n ) THEN
202 *
203 * Use unblocked code.
204 *
205  DO 20 j = n, 1, -1
206 *
207 * Copy current column of L to WORK and replace with zeros.
208 *
209  DO 10 i = j + 1, n
210  work( i ) = a( i, j )
211  a( i, j ) = zero
212  10 CONTINUE
213 *
214 * Compute current column of inv(A).
215 *
216  IF( j.LT.n )
217  $ CALL zgemv( 'No transpose', n, n-j, -one, a( 1, j+1 ),
218  $ lda, work( j+1 ), 1, one, a( 1, j ), 1 )
219  20 CONTINUE
220  ELSE
221 *
222 * Use blocked code.
223 *
224  nn = ( ( n-1 ) / nb )*nb + 1
225  DO 50 j = nn, 1, -nb
226  jb = min( nb, n-j+1 )
227 *
228 * Copy current block column of L to WORK and replace with
229 * zeros.
230 *
231  DO 40 jj = j, j + jb - 1
232  DO 30 i = jj + 1, n
233  work( i+( jj-j )*ldwork ) = a( i, jj )
234  a( i, jj ) = zero
235  30 CONTINUE
236  40 CONTINUE
237 *
238 * Compute current block column of inv(A).
239 *
240  IF( j+jb.LE.n )
241  $ CALL zgemm( 'No transpose', 'No transpose', n, jb,
242  $ n-j-jb+1, -one, a( 1, j+jb ), lda,
243  $ work( j+jb ), ldwork, one, a( 1, j ), lda )
244  CALL ztrsm( 'Right', 'Lower', 'No transpose', 'Unit', n, jb,
245  $ one, work( j ), ldwork, a( 1, j ), lda )
246  50 CONTINUE
247  END IF
248 *
249 * Apply column interchanges.
250 *
251  DO 60 j = n - 1, 1, -1
252  jp = ipiv( j )
253  IF( jp.NE.j )
254  $ CALL zswap( n, a( 1, j ), 1, a( 1, jp ), 1 )
255  60 CONTINUE
256 *
257  work( 1 ) = iws
258  RETURN
259 *
260 * End of ZGETRI
261 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine ztrtri(UPLO, DIAG, N, A, LDA, INFO)
ZTRTRI
Definition: ztrtri.f:111
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgetrs ( character  TRANS,
integer  N,
integer  NRHS,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

ZGETRS

Download ZGETRS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGETRS solves a system of linear equations
    A * X = B,  A**T * X = B,  or  A**H * X = B
 with a general N-by-N matrix A using the LU factorization computed
 by ZGETRF.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          Specifies the form of the system of equations:
          = 'N':  A * X = B     (No transpose)
          = 'T':  A**T * X = B  (Transpose)
          = 'C':  A**H * X = B  (Conjugate transpose)
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.
[in]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The factors L and U from the factorization A = P*L*U
          as computed by ZGETRF.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices from ZGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB,NRHS)
          On entry, the right hand side matrix B.
          On exit, the solution matrix X.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 123 of file zgetrs.f.

123 *
124 * -- LAPACK computational routine (version 3.4.0) --
125 * -- LAPACK is a software package provided by Univ. of Tennessee, --
126 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
127 * November 2011
128 *
129 * .. Scalar Arguments ..
130  CHARACTER trans
131  INTEGER info, lda, ldb, n, nrhs
132 * ..
133 * .. Array Arguments ..
134  INTEGER ipiv( * )
135  COMPLEX*16 a( lda, * ), b( ldb, * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  COMPLEX*16 one
142  parameter( one = ( 1.0d+0, 0.0d+0 ) )
143 * ..
144 * .. Local Scalars ..
145  LOGICAL notran
146 * ..
147 * .. External Functions ..
148  LOGICAL lsame
149  EXTERNAL lsame
150 * ..
151 * .. External Subroutines ..
152  EXTERNAL xerbla, zlaswp, ztrsm
153 * ..
154 * .. Intrinsic Functions ..
155  INTRINSIC max
156 * ..
157 * .. Executable Statements ..
158 *
159 * Test the input parameters.
160 *
161  info = 0
162  notran = lsame( trans, 'N' )
163  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
164  $ lsame( trans, 'C' ) ) THEN
165  info = -1
166  ELSE IF( n.LT.0 ) THEN
167  info = -2
168  ELSE IF( nrhs.LT.0 ) THEN
169  info = -3
170  ELSE IF( lda.LT.max( 1, n ) ) THEN
171  info = -5
172  ELSE IF( ldb.LT.max( 1, n ) ) THEN
173  info = -8
174  END IF
175  IF( info.NE.0 ) THEN
176  CALL xerbla( 'ZGETRS', -info )
177  RETURN
178  END IF
179 *
180 * Quick return if possible
181 *
182  IF( n.EQ.0 .OR. nrhs.EQ.0 )
183  $ RETURN
184 *
185  IF( notran ) THEN
186 *
187 * Solve A * X = B.
188 *
189 * Apply row interchanges to the right hand sides.
190 *
191  CALL zlaswp( nrhs, b, ldb, 1, n, ipiv, 1 )
192 *
193 * Solve L*X = B, overwriting B with X.
194 *
195  CALL ztrsm( 'Left', 'Lower', 'No transpose', 'Unit', n, nrhs,
196  $ one, a, lda, b, ldb )
197 *
198 * Solve U*X = B, overwriting B with X.
199 *
200  CALL ztrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', n,
201  $ nrhs, one, a, lda, b, ldb )
202  ELSE
203 *
204 * Solve A**T * X = B or A**H * X = B.
205 *
206 * Solve U**T *X = B or U**H *X = B, overwriting B with X.
207 *
208  CALL ztrsm( 'Left', 'Upper', trans, 'Non-unit', n, nrhs, one,
209  $ a, lda, b, ldb )
210 *
211 * Solve L**T *X = B, or L**H *X = B overwriting B with X.
212 *
213  CALL ztrsm( 'Left', 'Lower', trans, 'Unit', n, nrhs, one, a,
214  $ lda, b, ldb )
215 *
216 * Apply row interchanges to the solution vectors.
217 *
218  CALL zlaswp( nrhs, b, ldb, 1, n, ipiv, -1 )
219  END IF
220 *
221  RETURN
222 *
223 * End of ZGETRS
224 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
subroutine zlaswp(N, A, LDA, K1, K2, IPIV, INCX)
ZLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: zlaswp.f:116

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zhgeqz ( character  JOB,
character  COMPQ,
character  COMPZ,
integer  N,
integer  ILO,
integer  IHI,
complex*16, dimension( ldh, * )  H,
integer  LDH,
complex*16, dimension( ldt, * )  T,
integer  LDT,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldq, * )  Q,
integer  LDQ,
complex*16, dimension( ldz, * )  Z,
integer  LDZ,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZHGEQZ

Download ZHGEQZ + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZHGEQZ computes the eigenvalues of a complex matrix pair (H,T),
 where H is an upper Hessenberg matrix and T is upper triangular,
 using the single-shift QZ method.
 Matrix pairs of this type are produced by the reduction to
 generalized upper Hessenberg form of a complex matrix pair (A,B):
 
    A = Q1*H*Z1**H,  B = Q1*T*Z1**H,
 
 as computed by ZGGHRD.
 
 If JOB='S', then the Hessenberg-triangular pair (H,T) is
 also reduced to generalized Schur form,
 
    H = Q*S*Z**H,  T = Q*P*Z**H,
 
 where Q and Z are unitary matrices and S and P are upper triangular.
 
 Optionally, the unitary matrix Q from the generalized Schur
 factorization may be postmultiplied into an input matrix Q1, and the
 unitary matrix Z may be postmultiplied into an input matrix Z1.
 If Q1 and Z1 are the unitary matrices from ZGGHRD that reduced
 the matrix pair (A,B) to generalized Hessenberg form, then the output
 matrices Q1*Q and Z1*Z are the unitary factors from the generalized
 Schur factorization of (A,B):
 
    A = (Q1*Q)*S*(Z1*Z)**H,  B = (Q1*Q)*P*(Z1*Z)**H.
 
 To avoid overflow, eigenvalues of the matrix pair (H,T)
 (equivalently, of (A,B)) are computed as a pair of complex values
 (alpha,beta).  If beta is nonzero, lambda = alpha / beta is an
 eigenvalue of the generali