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

Functions

subroutine cgebak (JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
 CGEBAK More...
 
subroutine cgebal (JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
 CGEBAL More...
 
subroutine cgebd2 (M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO)
 CGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm. More...
 
subroutine cgebrd (M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
 CGEBRD More...
 
subroutine cgecon (NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
 CGECON More...
 
subroutine cgeequ (M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
 CGEEQU More...
 
subroutine cgeequb (M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
 CGEEQUB More...
 
subroutine cgehd2 (N, ILO, IHI, A, LDA, TAU, WORK, INFO)
 CGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm. More...
 
subroutine cgehrd (N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
 CGEHRD More...
 
subroutine cgelq2 (M, N, A, LDA, TAU, WORK, INFO)
 CGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine cgelqf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 CGELQF More...
 
subroutine cgemqrt (SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
 CGEMQRT More...
 
subroutine cgeql2 (M, N, A, LDA, TAU, WORK, INFO)
 CGEQL2 computes the QL factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine cgeqlf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 CGEQLF More...
 
subroutine cgeqp3 (M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
 CGEQP3 More...
 
subroutine cgeqr2 (M, N, A, LDA, TAU, WORK, INFO)
 CGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine cgeqr2p (M, N, A, LDA, TAU, WORK, INFO)
 CGEQR2P computes the QR factorization of a general rectangular matrix with non-negative diagonal elements using an unblocked algorithm. More...
 
subroutine cgeqrf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 CGEQRF More...
 
subroutine cgeqrfp (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 CGEQRFP More...
 
subroutine cgeqrt (M, N, NB, A, LDA, T, LDT, WORK, INFO)
 CGEQRT More...
 
subroutine cgeqrt2 (M, N, A, LDA, T, LDT, INFO)
 CGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY representation of Q. More...
 
recursive subroutine cgeqrt3 (M, N, A, LDA, T, LDT, INFO)
 CGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact WY representation of Q. More...
 
subroutine cgerfs (TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
 CGERFS More...
 
subroutine cgerfsx (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)
 CGERFSX More...
 
subroutine cgerq2 (M, N, A, LDA, TAU, WORK, INFO)
 CGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine cgerqf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 CGERQF More...
 
subroutine cgesvj (JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, CWORK, LWORK, RWORK, LRWORK, INFO)
 CGESVJ More...
 
subroutine cgetf2 (M, N, A, LDA, IPIV, INFO)
 CGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row interchanges (unblocked algorithm). More...
 
subroutine cgetrf (M, N, A, LDA, IPIV, INFO)
 CGETRF More...
 
recursive subroutine cgetrf2 (M, N, A, LDA, IPIV, INFO)
 CGETRF2 More...
 
subroutine cgetri (N, A, LDA, IPIV, WORK, LWORK, INFO)
 CGETRI More...
 
subroutine cgetrs (TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 CGETRS More...
 
subroutine chgeqz (JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
 CHGEQZ More...
 
subroutine cla_geamv (TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
 CLA_GEAMV computes a matrix-vector product using a general matrix to calculate error bounds. More...
 
real function cla_gercond_c (TRANS, N, A, LDA, AF, LDAF, IPIV, C, CAPPLY, INFO, WORK, RWORK)
 CLA_GERCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general matrices. More...
 
real function cla_gercond_x (TRANS, N, A, LDA, AF, LDAF, IPIV, X, INFO, WORK, RWORK)
 CLA_GERCOND_X computes the infinity norm condition number of op(A)*diag(x) for general matrices. More...
 
subroutine cla_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)
 CLA_GERFSX_EXTENDED More...
 
real function cla_gerpvgrw (N, NCOLS, A, LDA, AF, LDAF)
 CLA_GERPVGRW multiplies a square real matrix by a complex matrix. More...
 
subroutine ctgevc (SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
 CTGEVC More...
 
subroutine ctgexc (WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
 CTGEXC More...
 
subroutine cgeqpf (M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO)
 CGEQPF More...
 

Detailed Description

This is the group of complex computational functions for GE matrices

Function Documentation

subroutine cgebak ( character  JOB,
character  SIDE,
integer  N,
integer  ILO,
integer  IHI,
real, dimension( * )  SCALE,
integer  M,
complex, dimension( ldv, * )  V,
integer  LDV,
integer  INFO 
)

CGEBAK

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

Purpose:
 CGEBAK forms the right or left eigenvectors of a complex general
 matrix by backward transformation on the computed eigenvectors of the
 balanced matrix output by CGEBAL.
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 CGEBAL.
[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 CGEBAL.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
[in]SCALE
          SCALE is REAL array, dimension (N)
          Details of the permutation and scaling factors, as returned
          by CGEBAL.
[in]M
          M is INTEGER
          The number of columns of the matrix V.  M >= 0.
[in,out]V
          V is COMPLEX array, dimension (LDV,M)
          On entry, the matrix of right or left eigenvectors to be
          transformed, as returned by CHSEIN or CTREVC.
          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 cgebak.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  REAL scale( * )
145  COMPLEX v( ldv, * )
146 * ..
147 *
148 * =====================================================================
149 *
150 * .. Parameters ..
151  REAL one
152  parameter( one = 1.0e+0 )
153 * ..
154 * .. Local Scalars ..
155  LOGICAL leftv, rightv
156  INTEGER i, ii, k
157  REAL s
158 * ..
159 * .. External Functions ..
160  LOGICAL lsame
161  EXTERNAL lsame
162 * ..
163 * .. External Subroutines ..
164  EXTERNAL csscal, cswap, xerbla
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( 'CGEBAK', -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 csscal( 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 csscal( 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 cswap( 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 cswap( 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 CGEBAK
269 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine cgebal ( character  JOB,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
integer  ILO,
integer  IHI,
real, dimension( * )  SCALE,
integer  INFO 
)

CGEBAL

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

Purpose:
 CGEBAL 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 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
          ILO is INTEGER
[out]IHI
          IHI is INTEGER
          ILO and IHI are set to integers 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 REAL 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 163 of file cgebal.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine cgebd2 ( integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  D,
real, dimension( * )  E,
complex, dimension( * )  TAUQ,
complex, dimension( * )  TAUP,
complex, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 CGEBD2 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 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 REAL array, dimension (min(M,N))
          The diagonal elements of the bidiagonal matrix B:
          D(i) = A(i,i).
[out]E
          E is REAL 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 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 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 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 192 of file cgebd2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine cgebrd ( integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  D,
real, dimension( * )  E,
complex, dimension( * )  TAUQ,
complex, dimension( * )  TAUP,
complex, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

CGEBRD

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

Purpose:
 CGEBRD 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 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 REAL array, dimension (min(M,N))
          The diagonal elements of the bidiagonal matrix B:
          D(i) = A(i,i).
[out]E
          E is REAL 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 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 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 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 208 of file cgebrd.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine cgecon ( character  NORM,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real  ANORM,
real  RCOND,
complex, dimension( * )  WORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CGECON

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

Purpose:
 CGECON 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 CGETRF.

 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 array, dimension (LDA,N)
          The factors L and U from the factorization A = P*L*U
          as computed by CGETRF.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]ANORM
          ANORM is REAL
          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 REAL
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(norm(A) * norm(inv(A))).
[out]WORK
          WORK is COMPLEX array, dimension (2*N)
[out]RWORK
          RWORK is REAL 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 cgecon.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  REAL anorm, rcond
136 * ..
137 * .. Array Arguments ..
138  REAL rwork( * )
139  COMPLEX a( lda, * ), work( * )
140 * ..
141 *
142 * =====================================================================
143 *
144 * .. Parameters ..
145  REAL one, zero
146  parameter( one = 1.0e+0, zero = 0.0e+0 )
147 * ..
148 * .. Local Scalars ..
149  LOGICAL onenrm
150  CHARACTER normin
151  INTEGER ix, kase, kase1
152  REAL ainvnm, scale, sl, smlnum, su
153  COMPLEX zdum
154 * ..
155 * .. Local Arrays ..
156  INTEGER isave( 3 )
157 * ..
158 * .. External Functions ..
159  LOGICAL lsame
160  INTEGER icamax
161  REAL slamch
162  EXTERNAL lsame, icamax, slamch
163 * ..
164 * .. External Subroutines ..
165  EXTERNAL clacn2, clatrs, csrscl, xerbla
166 * ..
167 * .. Intrinsic Functions ..
168  INTRINSIC abs, aimag, max, real
169 * ..
170 * .. Statement Functions ..
171  REAL cabs1
172 * ..
173 * .. Statement Function definitions ..
174  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( 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( 'CGECON', -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 = slamch( '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 clacn2( 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 clatrs( 'Lower', 'No transpose', 'Unit', normin, n, a,
226  $ lda, work, sl, rwork, info )
227 *
228 * Multiply by inv(U).
229 *
230  CALL clatrs( '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 clatrs( '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 clatrs( '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 = icamax( n, work, 1 )
252  IF( scale.LT.cabs1( work( ix ) )*smlnum .OR. scale.EQ.zero )
253  $ GO TO 20
254  CALL csrscl( 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 CGECON
268 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:53
subroutine csrscl(N, SA, SX, INCX)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: csrscl.f:86
subroutine clatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
CLATRS solves a triangular system of equations with the scale factor set to prevent overflow...
Definition: clatrs.f:241
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine clacn2(N, V, X, EST, KASE, ISAVE)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: clacn2.f:135

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine cgeequ ( integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  R,
real, dimension( * )  C,
real  ROWCND,
real  COLCND,
real  AMAX,
integer  INFO 
)

CGEEQU

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

Purpose:
 CGEEQU 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 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 REAL array, dimension (M)
          If INFO = 0 or INFO > M, R contains the row scale factors
          for A.
[out]C
          C is REAL array, dimension (N)
          If INFO = 0,  C contains the column scale factors for A.
[out]ROWCND
          ROWCND is REAL
          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 REAL
          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 REAL
          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 cgeequ.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  REAL amax, colcnd, rowcnd
151 * ..
152 * .. Array Arguments ..
153  REAL c( * ), r( * )
154  COMPLEX a( lda, * )
155 * ..
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160  REAL one, zero
161  parameter( one = 1.0e+0, zero = 0.0e+0 )
162 * ..
163 * .. Local Scalars ..
164  INTEGER i, j
165  REAL bignum, rcmax, rcmin, smlnum
166  COMPLEX zdum
167 * ..
168 * .. External Functions ..
169  REAL slamch
170  EXTERNAL slamch
171 * ..
172 * .. External Subroutines ..
173  EXTERNAL xerbla
174 * ..
175 * .. Intrinsic Functions ..
176  INTRINSIC abs, aimag, max, min, real
177 * ..
178 * .. Statement Functions ..
179  REAL cabs1
180 * ..
181 * .. Statement Function definitions ..
182  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( 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( 'CGEEQU', -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 = slamch( '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 CGEEQU
312 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine cgeequb ( integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  R,
real, dimension( * )  C,
real  ROWCND,
real  COLCND,
real  AMAX,
integer  INFO 
)

CGEEQUB

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

Purpose:
 CGEEQUB 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 CGEEQU 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 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 REAL array, dimension (M)
          If INFO = 0 or INFO > M, R contains the row scale factors
          for A.
[out]C
          C is REAL array, dimension (N)
          If INFO = 0,  C contains the column scale factors for A.
[out]ROWCND
          ROWCND is REAL
          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 REAL
          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 REAL
          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 cgeequb.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  REAL amax, colcnd, rowcnd
158 * ..
159 * .. Array Arguments ..
160  REAL c( * ), r( * )
161  COMPLEX a( lda, * )
162 * ..
163 *
164 * =====================================================================
165 *
166 * .. Parameters ..
167  REAL one, zero
168  parameter( one = 1.0e+0, zero = 0.0e+0 )
169 * ..
170 * .. Local Scalars ..
171  INTEGER i, j
172  REAL bignum, rcmax, rcmin, smlnum, radix, logrdx
173  COMPLEX zdum
174 * ..
175 * .. External Functions ..
176  REAL slamch
177  EXTERNAL slamch
178 * ..
179 * .. External Subroutines ..
180  EXTERNAL xerbla
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC abs, max, min, log, REAL, aimag
184 * ..
185 * .. Statement Functions ..
186  REAL cabs1
187 * ..
188 * .. Statement Function definitions ..
189  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( 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( 'CGEEQUB', -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 = slamch( 'S' )
220  bignum = one / smlnum
221  radix = slamch( '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 CGEEQUB
329 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 CGEHD2 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 CGEBAL; 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 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 array, dimension (N-1)
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgehd2.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 a( lda, * ), tau( * ), work( * )
162 * ..
163 *
164 * =====================================================================
165 *
166 * .. Parameters ..
167  COMPLEX one
168  parameter( one = ( 1.0e+0, 0.0e+0 ) )
169 * ..
170 * .. Local Scalars ..
171  INTEGER i
172  COMPLEX alpha
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL clarf, clarfg, xerbla
176 * ..
177 * .. Intrinsic Functions ..
178  INTRINSIC conjg, 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( 'CGEHD2', -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 clarfg( 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 clarf( '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 clarf( 'Left', ihi-i, n-i, a( i+1, i ), 1,
215  $ conjg( 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 CGEHD2
223 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:130

Here is the call graph for this function:

Here is the caller graph for this function:

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

CGEHRD

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

Purpose:
 CGEHRD 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 CGEBAL; 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 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 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 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 cgehrd.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 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 zero, one
189  parameter( zero = ( 0.0e+0, 0.0e+0 ),
190  $ one = ( 1.0e+0, 0.0e+0 ) )
191 * ..
192 * .. Local Scalars ..
193  LOGICAL lquery
194  INTEGER i, ib, iinfo, iwt, j, ldwork, lwkopt, nb,
195  $ nbmin, nh, nx
196  COMPLEX ei
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL caxpy, cgehd2, cgemm, clahr2, clarfb, ctrmm,
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, 'CGEHRD', ' ', n, ilo, ihi, -1 ) )
232  lwkopt = n*nb + tsize
233  work( 1 ) = lwkopt
234  END IF
235 *
236  IF( info.NE.0 ) THEN
237  CALL xerbla( 'CGEHRD', -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, 'CGEHRD', ' ', 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, 'CGEHRD', ' ', 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, 'CGEHRD', ' ', 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 clahr2( 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 cgemm( '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 ctrmm( '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 caxpy( 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 clarfb( '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 cgehd2( n, i, ihi, a, lda, tau, work, iinfo )
350  work( 1 ) = lwkopt
351 *
352  RETURN
353 *
354 * End of CGEHRD
355 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: clarfb.f:197
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
subroutine cgehd2(N, ILO, IHI, A, LDA, TAU, WORK, INFO)
CGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm...
Definition: cgehd2.f:151
subroutine clahr2(N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)
CLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elemen...
Definition: clahr2.f:183
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:179

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 CGELQ2 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgelq2.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 a( lda, * ), tau( * ), work( * )
134 * ..
135 *
136 * =====================================================================
137 *
138 * .. Parameters ..
139  COMPLEX one
140  parameter( one = ( 1.0e+0, 0.0e+0 ) )
141 * ..
142 * .. Local Scalars ..
143  INTEGER i, k
144  COMPLEX alpha
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL clacgv, clarf, clarfg, xerbla
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( 'CGELQ2', -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 clacgv( n-i+1, a( i, i ), lda )
176  alpha = a( i, i )
177  CALL clarfg( 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 clarf( '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 clacgv( n-i+1, a( i, i ), lda )
189  10 CONTINUE
190  RETURN
191 *
192 * End of CGELQ2
193 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:76
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:130

Here is the call graph for this function:

Here is the caller graph for this function:

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

CGELQF

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

Purpose:
 CGELQF 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgelqf.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 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 cgelq2, clarfb, clarft, xerbla
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, 'CGELQF', ' ', 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( 'CGELQF', -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, 'CGELQF', ' ', 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, 'CGELQF', ' ', 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 cgelq2( 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 clarft( '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 clarfb( '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 cgelq2( 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 CGELQF
268 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: clarfb.f:197
subroutine clarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
CLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: clarft.f:165
subroutine cgelq2(M, N, A, LDA, TAU, WORK, INFO)
CGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: cgelq2.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 cgemqrt ( character  SIDE,
character  TRANS,
integer  M,
integer  N,
integer  K,
integer  NB,
complex, dimension( ldv, * )  V,
integer  LDV,
complex, dimension( ldt, * )  T,
integer  LDT,
complex, dimension( ldc, * )  C,
integer  LDC,
complex, dimension( * )  WORK,
integer  INFO 
)

CGEMQRT

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

Purpose:
 CGEMQRT 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 CGEQRT. 

 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 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 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 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 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 cgemqrt.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 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, clarfb
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( 'CGEMQRT', -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 clarfb( '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 clarfb( '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 clarfb( '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 clarfb( '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 CGEMQRT
290 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: clarfb.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 cgeql2 ( integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( * )  TAU,
complex, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 CGEQL2 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgeql2.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 a( lda, * ), tau( * ), work( * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  COMPLEX one
142  parameter( one = ( 1.0e+0, 0.0e+0 ) )
143 * ..
144 * .. Local Scalars ..
145  INTEGER i, k
146  COMPLEX alpha
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL clarf, clarfg, xerbla
150 * ..
151 * .. Intrinsic Functions ..
152  INTRINSIC conjg, 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( 'CGEQL2', -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 clarfg( 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 clarf( 'Left', m-k+i, n-k+i-1, a( 1, n-k+i ), 1,
185  $ conjg( tau( i ) ), a, lda, work )
186  a( m-k+i, n-k+i ) = alpha
187  10 CONTINUE
188  RETURN
189 *
190 * End of CGEQL2
191 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:130

Here is the call graph for this function:

Here is the caller graph for this function:

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

CGEQLF

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

Purpose:
 CGEQLF 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgeqlf.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 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 cgeql2, clarfb, clarft, xerbla
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, 'CGEQLF', ' ', 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( 'CGEQLF', -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, 'CGEQLF', ' ', 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, 'CGEQLF', ' ', 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 cgeql2( 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 clarft( '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 clarfb( '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 cgeql2( mu, nu, a, lda, tau, work, iinfo )
281 *
282  work( 1 ) = iws
283  RETURN
284 *
285 * End of CGEQLF
286 *
subroutine cgeql2(M, N, A, LDA, TAU, WORK, INFO)
CGEQL2 computes the QL factorization of a general rectangular matrix using an unblocked algorithm...
Definition: cgeql2.f:125
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: clarfb.f:197
subroutine clarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
CLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: clarft.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 cgeqp3 ( integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
complex, dimension( * )  TAU,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CGEQP3

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

Purpose:
 CGEQP3 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[out]WORK
          WORK is COMPLEX 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 REAL 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 cgeqp3.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  REAL rwork( * )
173  COMPLEX 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 cgeqrf, claqp2, claqps, cswap, cunmqr, xerbla
189 * ..
190 * .. External Functions ..
191  INTEGER ilaenv
192  REAL scnrm2
193  EXTERNAL ilaenv, scnrm2
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, 'CGEQRF', ' ', m, n, -1, -1 )
221  lwkopt = ( n + 1 )*nb
222  END IF
223  work( 1 ) = cmplx( 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( 'CGEQP3', -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 cswap( 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 CGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
265  CALL cgeqrf( m, na, a, lda, tau, work, lwork, info )
266  iws = max( iws, int( work( 1 ) ) )
267  IF( na.LT.n ) THEN
268 *CC CALL CUNM2R( 'Left', 'Conjugate Transpose', M, N-NA,
269 *CC $ NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK,
270 *CC $ INFO )
271  CALL cunmqr( '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, 'CGEQRF', ' ', 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, 'CGEQRF', ' ', 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, 'CGEQRF', ' ', 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 ) = scnrm2( 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 claqps( 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 claqp2( 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 ) = cmplx( lwkopt )
368  RETURN
369 *
370 * End of CGEQP3
371 *
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:138
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine claqp2(M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)
CLAQP2 computes a QR factorization with column pivoting of the matrix block.
Definition: claqp2.f:151
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:170
subroutine claqps(M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)
CLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BL...
Definition: claqps.f:180
real function scnrm2(N, X, INCX)
SCNRM2
Definition: scnrm2.f:56
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine cgeqpf ( integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
complex, dimension( * )  TAU,
complex, dimension( * )  WORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CGEQPF

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

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

 CGEQPF 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[out]WORK
          WORK is COMPLEX array, dimension (N)
[out]RWORK
          RWORK is REAL 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 cgeqpf.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  REAL rwork( * )
162  COMPLEX a( lda, * ), tau( * ), work( * )
163 * ..
164 *
165 * =====================================================================
166 *
167 * .. Parameters ..
168  REAL zero, one
169  parameter( zero = 0.0e+0, one = 1.0e+0 )
170 * ..
171 * .. Local Scalars ..
172  INTEGER i, itemp, j, ma, mn, pvt
173  REAL temp, temp2, tol3z
174  COMPLEX aii
175 * ..
176 * .. External Subroutines ..
177  EXTERNAL cgeqr2, clarf, clarfg, cswap, cunm2r, xerbla
178 * ..
179 * .. Intrinsic Functions ..
180  INTRINSIC abs, cmplx, conjg, max, min, sqrt
181 * ..
182 * .. External Functions ..
183  INTEGER isamax
184  REAL scnrm2, slamch
185  EXTERNAL isamax, scnrm2, slamch
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( 'CGEQPF', -info )
201  RETURN
202  END IF
203 *
204  mn = min( m, n )
205  tol3z = sqrt(slamch('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 cswap( 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 cgeqr2( m, ma, a, lda, tau, work, info )
231  IF( ma.LT.n ) THEN
232  CALL cunm2r( '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 ) = scnrm2( 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 ) + isamax( n-i+1, rwork( i ), 1 )
254 *
255  IF( pvt.NE.i ) THEN
256  CALL cswap( 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 clarfg( 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 ) = cmplx( one )
277  CALL clarf( 'Left', m-i+1, n-i, a( i, i ), 1,
278  $ conjg( 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 ) = scnrm2( 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 CGEQPF
312 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
real function scnrm2(N, X, INCX)
SCNRM2
Definition: scnrm2.f:56
subroutine cgeqr2(M, N, A, LDA, TAU, WORK, INFO)
CGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: cgeqr2.f:123
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:130
subroutine cunm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
CUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition: cunm2r.f:161
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 CGEQR2 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgeqr2.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 a( lda, * ), tau( * ), work( * )
134 * ..
135 *
136 * =====================================================================
137 *
138 * .. Parameters ..
139  COMPLEX one
140  parameter( one = ( 1.0e+0, 0.0e+0 ) )
141 * ..
142 * .. Local Scalars ..
143  INTEGER i, k
144  COMPLEX alpha
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL clarf, clarfg, xerbla
148 * ..
149 * .. Intrinsic Functions ..
150  INTRINSIC conjg, 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( 'CGEQR2', -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 clarfg( 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 clarf( 'Left', m-i+1, n-i, a( i, i ), 1,
184  $ conjg( 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 CGEQR2
191 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:130

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 CGEQR2P 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgeqr2p.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 a( lda, * ), tau( * ), work( * )
137 * ..
138 *
139 * =====================================================================
140 *
141 * .. Parameters ..
142  COMPLEX one
143  parameter( one = ( 1.0e+0, 0.0e+0 ) )
144 * ..
145 * .. Local Scalars ..
146  INTEGER i, k
147  COMPLEX alpha
148 * ..
149 * .. External Subroutines ..
150  EXTERNAL clarf, clarfgp, xerbla
151 * ..
152 * .. Intrinsic Functions ..
153  INTRINSIC conjg, 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( 'CGEQR2P', -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 clarfgp( 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 clarf( 'Left', m-i+1, n-i, a( i, i ), 1,
187  $ conjg( 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 CGEQR2P
194 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfgp(N, ALPHA, X, INCX, TAU)
CLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: clarfgp.f:106
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:130

Here is the call graph for this function:

Here is the caller graph for this function:

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

CGEQRF

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

Purpose:
 CGEQRF 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgeqrf.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 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 cgeqr2, clarfb, clarft, xerbla
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, 'CGEQRF', ' ', 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( 'CGEQRF', -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, 'CGEQRF', ' ', 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, 'CGEQRF', ' ', 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 cgeqr2( 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 clarft( '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 clarfb( '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 cgeqr2( 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 CGEQRF
269 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: clarfb.f:197
subroutine clarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
CLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: clarft.f:165
subroutine cgeqr2(M, N, A, LDA, TAU, WORK, INFO)
CGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: cgeqr2.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 cgeqrfp ( integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( * )  TAU,
complex, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

CGEQRFP

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

Purpose:
 CGEQRFP 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgeqrfp.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 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 cgeqr2p, clarfb, clarft, xerbla
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, 'CGEQRF', ' ', 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( 'CGEQRFP', -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, 'CGEQRF', ' ', 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, 'CGEQRF', ' ', 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 cgeqr2p( 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 clarft( '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 clarfb( '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 cgeqr2p( 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 CGEQRFP
272 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: clarfb.f:197
subroutine cgeqr2p(M, N, A, LDA, TAU, WORK, INFO)
CGEQR2P computes the QR factorization of a general rectangular matrix with non-negative diagonal elem...
Definition: cgeqr2p.f:126
subroutine clarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
CLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: clarft.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 cgeqrt ( integer  M,
integer  N,
integer  NB,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldt, * )  T,
integer  LDT,
complex, dimension( * )  WORK,
integer  INFO 
)

CGEQRT

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

Purpose:
 CGEQRT 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 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 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 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 cgeqrt.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 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 cgeqrt2, cgeqrt3, clarfb, 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( 'CGEQRT', -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 cgeqrt3( m-i+1, ib, a(i,i), lda, t(1,i), ldt, iinfo )
202  ELSE
203  CALL cgeqrt2( 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 clarfb( '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 CGEQRT
217 *
recursive subroutine cgeqrt3(M, N, A, LDA, T, LDT, INFO)
CGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact...
Definition: cgeqrt3.f:134
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: clarfb.f:197
subroutine cgeqrt2(M, N, A, LDA, T, LDT, INFO)
CGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY represen...
Definition: cgeqrt2.f:129

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 CGEQRT2 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 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 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 cgeqrt2.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 a( lda, * ), t( ldt, * )
140 * ..
141 *
142 * =====================================================================
143 *
144 * .. Parameters ..
145  COMPLEX one, zero
146  parameter( one = (1.0,0.0), zero = (0.0,0.0) )
147 * ..
148 * .. Local Scalars ..
149  INTEGER i, k
150  COMPLEX aii, alpha
151 * ..
152 * .. External Subroutines ..
153  EXTERNAL clarfg, cgemv, cgerc, ctrmv, 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( 'CGEQRT2', -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 clarfg( 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 cgemv( '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 cgerc( 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 cgemv( '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 ctrmv( '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 CGEQRT2
226 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine ctrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRMV
Definition: ctrmv.f:149
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
subroutine cgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERC
Definition: cgerc.f:132

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 CGEQRT3 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 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 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 cgeqrt3.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 a( lda, * ), t( ldt, * )
145 * ..
146 *
147 * =====================================================================
148 *
149 * .. Parameters ..
150  COMPLEX one
151  parameter( one = (1.0,0.0) )
152 * ..
153 * .. Local Scalars ..
154  INTEGER i, i1, j, j1, n1, n2, iinfo
155 * ..
156 * .. External Subroutines ..
157  EXTERNAL clarfg, ctrmm, cgemm, 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( 'CGEQRT3', -info )
173  RETURN
174  END IF
175 *
176  IF( n.EQ.1 ) THEN
177 *
178 * Compute Householder transform when N=1
179 *
180  CALL clarfg( 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 cgeqrt3( 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 ctrmm( 'L', 'L', 'C', 'U', n1, n2, one,
203  & a, lda, t( 1, j1 ), ldt )
204 *
205  CALL cgemm( 'C', 'N', n1, n2, m-n1, one, a( j1, 1 ), lda,
206  & a( j1, j1 ), lda, one, t( 1, j1 ), ldt)
207 *
208  CALL ctrmm( 'L', 'U', 'C', 'N', n1, n2, one,
209  & t, ldt, t( 1, j1 ), ldt )
210 *
211  CALL cgemm( 'N', 'N', m-n1, n2, n1, -one, a( j1, 1 ), lda,
212  & t( 1, j1 ), ldt, one, a( j1, j1 ), lda )
213 *
214  CALL ctrmm( '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 cgeqrt3( 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 ctrmm( 'R', 'L', 'N', 'U', n1, n2, one,
237  & a( j1, j1 ), lda, t( 1, j1 ), ldt )
238 *
239  CALL cgemm( 'C', 'N', n1, n2, m-n, one, a( i1, 1 ), lda,
240  & a( i1, j1 ), lda, one, t( 1, j1 ), ldt )
241 *
242  CALL ctrmm( 'L', 'U', 'N', 'N', n1, n2, -one, t, ldt,
243  & t( 1, j1 ), ldt )
244 *
245  CALL ctrmm( '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 CGEQRT3
256 *
recursive subroutine cgeqrt3(M, N, A, LDA, T, LDT, INFO)
CGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact...
Definition: cgeqrt3.f:134
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:179

Here is the call graph for this function:

Here is the caller graph for this function:

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

CGERFS

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

Purpose:
 CGERFS 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 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 array, dimension (LDAF,N)
          The factors L and U from the factorization A = P*L*U
          as computed by CGETRF.
[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 CGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[in]B
          B is COMPLEX 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 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by CGETRS.
          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 REAL 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 REAL 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 array, dimension (2*N)
[out]RWORK
          RWORK is REAL 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 cgerfs.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  REAL berr( * ), ferr( * ), rwork( * )
201  COMPLEX a( lda, * ), af( ldaf, * ), b( ldb, * ),
202  $ work( * ), x( ldx, * )
203 * ..
204 *
205 * =====================================================================
206 *
207 * .. Parameters ..
208  INTEGER itmax
209  parameter( itmax = 5 )
210  REAL zero
211  parameter( zero = 0.0e+0 )
212  COMPLEX one
213  parameter( one = ( 1.0e+0, 0.0e+0 ) )
214  REAL two
215  parameter( two = 2.0e+0 )
216  REAL three
217  parameter( three = 3.0e+0 )
218 * ..
219 * .. Local Scalars ..
220  LOGICAL notran
221  CHARACTER transn, transt
222  INTEGER count, i, j, k, kase, nz
223  REAL eps, lstres, s, safe1, safe2, safmin, xk
224  COMPLEX zdum
225 * ..
226 * .. Local Arrays ..
227  INTEGER isave( 3 )
228 * ..
229 * .. External Functions ..
230  LOGICAL lsame
231  REAL slamch
232  EXTERNAL lsame, slamch
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL caxpy, ccopy, cgemv, cgetrs, clacn2, xerbla
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, aimag, max, real
239 * ..
240 * .. Statement Functions ..
241  REAL cabs1
242 * ..
243 * .. Statement Function definitions ..
244  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( 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( 'CGERFS', -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 = slamch( 'Epsilon' )
295  safmin = slamch( '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 ccopy( n, b( 1, j ), 1, work, 1 )
313  CALL cgemv( 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 cgetrs( trans, n, 1, af, ldaf, ipiv, work, n, info )
370  CALL caxpy( 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 CLACN2 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 clacn2( 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 cgetrs( 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 cgetrs( 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 CGERFS
447 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine cgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CGETRS
Definition: cgetrs.f:123
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
subroutine clacn2(N, V, X, EST, KASE, ISAVE)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: clacn2.f:135

Here is the call graph for this function:

Here is the caller graph for this function:

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

CGERFSX

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

Purpose:
    CGERFSX 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 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 array, dimension (LDAF,N)
     The factors L and U from the factorization A = P*L*U
     as computed by CGETRF.
[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 CGETRF; for 1<=i<=N, row i of the
     matrix was interchanged with row IPIV(i).
[in]R
          R is REAL 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 REAL 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 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 array, dimension (LDX,NRHS)
     On entry, the solution matrix X, as computed by CGETRS.
     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 REAL
     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 REAL 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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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.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 array, dimension (2*N)
[out]RWORK
          RWORK is REAL 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 cgerfsx.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  REAL rcond
427 * ..
428 * .. Array Arguments ..
429  INTEGER ipiv( * )
430  COMPLEX a( lda, * ), af( ldaf, * ), b( ldb, * ),
431  $ x( ldx , * ), work( * )
432  REAL r( * ), c( * ), params( * ), berr( * ),
433  $ err_bnds_norm( nrhs, * ),
434  $ err_bnds_comp( nrhs, * ), rwork( * )
435 * ..
436 *
437 * ==================================================================
438 *
439 * .. Parameters ..
440  REAL zero, one
441  parameter( zero = 0.0e+0, one = 1.0e+0 )
442  REAL itref_default, ithresh_default,
443  $ componentwise_default
444  REAL rthresh_default, dzthresh_default
445  parameter( itref_default = 1.0 )
446  parameter( ithresh_default = 10.0 )
447  parameter( componentwise_default = 1.0 )
448  parameter( rthresh_default = 0.5 )
449  parameter( dzthresh_default = 0.25 )
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  REAL anorm, rcond_tmp
466  REAL illrcond_thresh, err_lbnd, cwise_wrong
467  LOGICAL ignore_cwise
468  INTEGER ithresh
469  REAL 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
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.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 = REAL( N ) * slamch( 'Epsilon' )
503  ithresh = int( ithresh_default )
504  rthresh = rthresh_default
505  unstable_thresh = dzthresh_default
506  ignore_cwise = componentwise_default .EQ. 0.0
507 *
508  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
509  IF ( params( la_linrx_ithresh_i ).LT.0.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.0 ) THEN
517  IF ( ignore_cwise ) THEN
518  params( la_linrx_cwise_i ) = 0.0
519  ELSE
520  params( la_linrx_cwise_i ) = 1.0
521  END IF
522  ELSE
523  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.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( 'CGERFSX', -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.0
567  DO j = 1, nrhs
568  berr( j ) = 0.0
569  IF ( n_err_bnds .GE. 1 ) THEN
570  err_bnds_norm( j, la_linrx_trust_i ) = 1.0
571  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
572  END IF
573  IF ( n_err_bnds .GE. 2 ) THEN
574  err_bnds_norm( j, la_linrx_err_i ) = 0.0
575  err_bnds_comp( j, la_linrx_err_i ) = 0.0
576  END IF
577  IF ( n_err_bnds .GE. 3 ) THEN
578  err_bnds_norm( j, la_linrx_rcond_i ) = 1.0
579  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0
580  END IF
581  END DO
582  RETURN
583  END IF
584 *
585 * Default to failure.
586 *
587  rcond = 0.0
588  DO j = 1, nrhs
589  berr( j ) = 1.0
590  IF ( n_err_bnds .GE. 1 ) THEN
591  err_bnds_norm( j, la_linrx_trust_i ) = 1.0
592  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
593  END IF
594  IF ( n_err_bnds .GE. 2 ) THEN
595  err_bnds_norm( j, la_linrx_err_i ) = 1.0
596  err_bnds_comp( j, la_linrx_err_i ) = 1.0
597  END IF
598  IF ( n_err_bnds .GE. 3 ) THEN
599  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0
600  err_bnds_comp( j, la_linrx_rcond_i ) = 0.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 = clange( norm, n, n, a, lda, rwork )
613  CALL cgecon( 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( 'D' )
620 
621  IF ( notran ) THEN
622  CALL cla_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 cla_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.0, sqrt( REAL( N ) ) ) * slamch( '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 = cla_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 = cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv,
650  $ r, .true., info, work, rwork )
651  ELSE
652  rcond_tmp = cla_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.0 )
661  $ err_bnds_norm( j, la_linrx_err_i ) = 1.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.0
667  err_bnds_norm( j, la_linrx_trust_i ) = 0.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.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( slamch( 'Epsilon' ) )
694  DO j = 1, nrhs
695  IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
696  $ THEN
697  rcond_tmp = cla_gercond_x( trans, n, a, lda, af, ldaf,
698  $ ipiv, x(1,j), info, work, rwork )
699  ELSE
700  rcond_tmp = 0.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.0 )
707  $ err_bnds_comp( j, la_linrx_err_i ) = 1.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.0
713  err_bnds_comp( j, la_linrx_trust_i ) = 0.0
714  IF ( params( la_linrx_cwise_i ) .EQ. 1.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.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 CGERFSX
734 *
real function cla_gercond_x(TRANS, N, A, LDA, AF, LDAF, IPIV, X, INFO, WORK, RWORK)
CLA_GERCOND_X computes the infinity norm condition number of op(A)*diag(x) for general matrices...
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:60
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine cgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
CGECON
Definition: cgecon.f:126
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
real function cla_gercond_c(TRANS, N, A, LDA, AF, LDAF, IPIV, C, CAPPLY, INFO, WORK, RWORK)
CLA_GERCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general matrices...
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:60
subroutine cla_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)
CLA_GERFSX_EXTENDED
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:117

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 CGERQ2 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgerq2.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 a( lda, * ), tau( * ), work( * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  COMPLEX one
142  parameter( one = ( 1.0e+0, 0.0e+0 ) )
143 * ..
144 * .. Local Scalars ..
145  INTEGER i, k
146  COMPLEX alpha
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL clacgv, clarf, clarfg, xerbla
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( 'CGERQ2', -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 clacgv( n-k+i, a( m-k+i, 1 ), lda )
179  alpha = a( m-k+i, n-k+i )
180  CALL clarfg( n-k+i, alpha, a( m-k+i, 1 ), lda,
181  $ tau( i ) )
182 *
183 * Apply H(i) to A(1:m-k+i-1,1:n-k+i) from the right
184 *
185  a( m-k+i, n-k+i ) = one
186  CALL clarf( 'Right', m-k+i-1, n-k+i, a( m-k+i, 1 ), lda,
187  $ tau( i ), a, lda, work )
188  a( m-k+i, n-k+i ) = alpha
189  CALL clacgv( n-k+i-1, a( m-k+i, 1 ), lda )
190  10 CONTINUE
191  RETURN
192 *
193 * End of CGERQ2
194 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:76
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:130

Here is the call graph for this function:

Here is the caller graph for this function:

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

CGERQF

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

Purpose:
 CGERQF 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 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 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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 cgerqf.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 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 cgerq2, clarfb, clarft, xerbla
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, 'CGERQF', ' ', 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( 'CGERQF', -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, 'CGERQF', ' ', 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, 'CGERQF', ' ', 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 cgerq2( 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 clarft( '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 clarfb( '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 cgerq2( mu, nu, a, lda, tau, work, iinfo )
281 *
282  work( 1 ) = iws
283  RETURN
284 *
285 * End of CGERQF
286 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: clarfb.f:197
subroutine clarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
CLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: clarft.f:165
subroutine cgerq2(M, N, A, LDA, TAU, WORK, INFO)
CGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: cgerq2.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 cgesvj ( character*1  JOBA,
character*1  JOBU,
character*1  JOBV,
integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real, dimension( n )  SVA,
integer  MV,
complex, dimension( ldv, * )  V,
integer  LDV,
complex, dimension( lwork )  CWORK,
integer  LWORK,
real, dimension( lrwork )  RWORK,
integer  LRWORK,
integer  INFO 
)

CGESVJ

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

Purpose:
 
Parameters
[in]JOBA
          JOBA is CHARACTER* 1
          Specifies the structure of A.
          = 'L': The input matrix A is lower triangular;
          = 'U': The input matrix A is upper triangular;
          = 'G': The input matrix A is general M-by-N matrix, M >= N.
[in]JOBU
          JOBU is CHARACTER*1
          Specifies whether to compute the left singular vectors
          (columns of U):
          = 'U': The left singular vectors corresponding to the nonzero
                 singular values are computed and returned in the leading
                 columns of A. See more details in the description of A.
                 The default numerical orthogonality threshold is set to
                 approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=SLAMCH('E').
          = 'C': Analogous to JOBU='U', except that user can control the
                 level of numerical orthogonality of the computed left
                 singular vectors. TOL can be set to TOL = CTOL*EPS, where
                 CTOL is given on input in the array WORK.
                 No CTOL smaller than ONE is allowed. CTOL greater
                 than 1 / EPS is meaningless. The option 'C'
                 can be used if M*EPS is satisfactory orthogonality
                 of the computed left singular vectors, so CTOL=M could
                 save few sweeps of Jacobi rotations.
                 See the descriptions of A and WORK(1).
          = 'N': The matrix U is not computed. However, see the
                 description of A.
[in]JOBV
          JOBV is CHARACTER*1
          Specifies whether to compute the right singular vectors, that
          is, the matrix V:
          = 'V' : the matrix V is computed and returned in the array V
          = 'A' : the Jacobi rotations are applied to the MV-by-N
                  array V. In other words, the right singular vector
                  matrix V is not computed explicitly; instead it is
                  applied to an MV-by-N matrix initially stored in the
                  first MV rows of V.
          = 'N' : the matrix V is not computed and the array V is not
                  referenced
[in]M
          M is INTEGER
          The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0.
[in]N
          N is INTEGER
          The number of columns of the input matrix A.
          M >= N >= 0.
[in,out]A
          A is COMPLEX array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit,
          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C':
                 If INFO .EQ. 0 :
                 RANKA orthonormal columns of U are returned in the
                 leading RANKA columns of the array A. Here RANKA <= N
                 is the number of computed singular values of A that are
                 above the underflow threshold SLAMCH('S'). The singular
                 vectors corresponding to underflowed or zero singular
                 values are not computed. The value of RANKA is returned
                 in the array RWORK as RANKA=NINT(RWORK(2)). Also see the
                 descriptions of SVA and RWORK. The computed columns of U
                 are mutually numerically orthogonal up to approximately
                 TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
                 see the description of JOBU.
                 If INFO .GT. 0,
                 the procedure CGESVJ did not converge in the given number
                 of iterations (sweeps). In that case, the computed
                 columns of U may not be orthogonal up to TOL. The output
                 U (stored in A), SIGMA (given by the computed singular
                 values in SVA(1:N)) and V is still a decomposition of the
                 input matrix A in the sense that the residual
                 || A - SCALE * U * SIGMA * V^* ||_2 / ||A||_2 is small.
          If JOBU .EQ. 'N':
                 If INFO .EQ. 0 :
                 Note that the left singular vectors are 'for free' in the
                 one-sided Jacobi SVD algorithm. However, if only the
                 singular values are needed, the level of numerical
                 orthogonality of U is not an issue and iterations are
                 stopped when the columns of the iterated matrix are
                 numerically orthogonal up to approximately M*EPS. Thus,
                 on exit, A contains the columns of U scaled with the
                 corresponding singular values.
                 If INFO .GT. 0 :
                 the procedure CGESVJ did not converge in the given number
                 of iterations (sweeps).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]SVA
          SVA is REAL array, dimension (N)
          On exit,
          If INFO .EQ. 0 :
          depending on the value SCALE = RWORK(1), we have:
                 If SCALE .EQ. ONE:
                 SVA(1:N) contains the computed singular values of A.
                 During the computation SVA contains the Euclidean column
                 norms of the iterated matrices in the array A.
                 If SCALE .NE. ONE:
                 The singular values of A are SCALE*SVA(1:N), and this
                 factored representation is due to the fact that some of the
                 singular values of A might underflow or overflow.

          If INFO .GT. 0 :
          the procedure CGESVJ did not converge in the given number of
          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
[in]MV
          MV is INTEGER
          If JOBV .EQ. 'A', then the product of Jacobi rotations in CGESVJ
          is applied to the first MV rows of V. See the description of JOBV.
[in,out]V
          V is COMPLEX array, dimension (LDV,N)
          If JOBV = 'V', then V contains on exit the N-by-N matrix of
                         the right singular vectors;
          If JOBV = 'A', then V contains the product of the computed right
                         singular vector matrix and the initial matrix in
                         the array V.
          If JOBV = 'N', then V is not referenced.
[in]LDV
          LDV is INTEGER
          The leading dimension of the array V, LDV .GE. 1.
          If JOBV .EQ. 'V', then LDV .GE. max(1,N).
          If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .
[in,out]CWORKCWORK is COMPLEX array, dimension M+N. Used as work space.
[in]LWORK
          LWORK is INTEGER
          Length of CWORK, LWORK >= M+N.

 \param[in,out] RWORK
          RWORK is REAL array, dimension max(6,M+N).
          On entry,
          If JOBU .EQ. 'C' :
          RWORK(1) = CTOL, where CTOL defines the threshold for convergence.
                    The process stops if all columns of A are mutually
                    orthogonal up to CTOL*EPS, EPS=SLAMCH('E').
                    It is required that CTOL >= ONE, i.e. it is not
                    allowed to force the routine to obtain orthogonality
                    below EPSILON.
          On exit,
          RWORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
                    are the computed singular values of A.
                    (See description of SVA().)
          RWORK(2) = NINT(RWORK(2)) is the number of the computed nonzero
                    singular values.
          RWORK(3) = NINT(RWORK(3)) is the number of the computed singular
                    values that are larger than the underflow threshold.
          RWORK(4) = NINT(RWORK(4)) is the number of sweeps of Jacobi
                    rotations needed for numerical convergence.
          RWORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
                    This is useful information in cases when CGESVJ did
                    not converge, as it can be used to estimate whether
                    the output is stil useful and for post festum analysis.
          RWORK(6) = the largest absolute value over all sines of the
                    Jacobi rotation angles in the last sweep. It can be
                    useful for a post festum analysis.
[in]LRWORKLRWORK is INTEGER Length of RWORK, LRWORK >= MAX(6,N).
[out]INFO
          INFO is INTEGER
          = 0 : successful exit.
          < 0 : if INFO = -i, then the i-th argument had an illegal value
          > 0 : CGESVJ did not converge in the maximal allowed number 
                (NSWEEP=30) of sweeps. The output may still be useful. 
                See the description of RWORK.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Further Details:
The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane rotations. In the case of underflow of the tangent of the Jacobi angle, a modified Jacobi transformation of Drmac [3] is used. Pivot strategy uses column interchanges of de Rijk [1]. The relative accuracy of the computed singular values and the accuracy of the computed singular vectors (in angle metric) is as guaranteed by the theory of Demmel and Veselic [2]. The condition number that determines the accuracy in the full rank case is essentially min_{D=diag} kappa(A*D), where kappa(.) is the spectral condition number. The best performance of this Jacobi SVD procedure is achieved if used in an accelerated version of Drmac and Veselic [4,5], and it is the kernel routine in the SIGMA library [6]. Some tunning parameters (marked with [TP]) are available for the implementer. The computational range for the nonzero singular values is the machine number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even denormalized singular values can be computed with the corresponding gradual loss of accurate digits.
Contributors:
Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
References:
[1] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the singular value decomposition on a vector computer. SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. [2] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. [3] Z. Drmac: Implementation of Jacobi rotations for accurate singular value computation in floating point arithmetic. SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222. [4] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. LAPACK Working note 169. [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. LAPACK Working note 170. [6] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, QSVD, (H,K)-SVD computations. Department of Mathematics, University of Zagreb, 2008, 2015.
Bugs, Examples and Comments:
Please report all bugs and send interesting test examples and comments to drmac.nosp@m.@mat.nosp@m.h.hr. Thank you.

Definition at line 328 of file cgesvj.f.

328 *
329 * -- LAPACK computational routine (version 3.6.0) --
330 * -- LAPACK is a software package provided by Univ. of Tennessee, --
331 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
332 * November 2015
333 *
334  IMPLICIT NONE
335 * .. Scalar Arguments ..
336  INTEGER info, lda, ldv, lwork, lrwork, m, mv, n
337  CHARACTER*1 joba, jobu, jobv
338 * ..
339 * .. Array Arguments ..
340  COMPLEX a( lda, * ), v( ldv, * ), cwork( lwork )
341  REAL rwork( lrwork ), sva( n )
342 * ..
343 *
344 * =====================================================================
345 *
346 * .. Local Parameters ..
347  REAL zero, half, one
348  parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
349  COMPLEX czero, cone
350  parameter( czero = (0.0e0, 0.0e0), cone = (1.0e0, 0.0e0) )
351  INTEGER nsweep
352  parameter( nsweep = 30 )
353 * ..
354 * .. Local Scalars ..
355  COMPLEX aapq, ompq
356  REAL aapp, aapp0, aapq1, aaqq, apoaq, aqoap, big,
357  $ bigtheta, cs, ctol, epsln, large, mxaapq,
358  $ mxsinj, rootbig, rooteps, rootsfmin, roottol,
359  $ skl, sfmin, small, sn, t, temp1, theta, thsign, tol
360  INTEGER blskip, emptsw, i, ibr, ierr, igl, ijblsk, ir1,
361  $ iswrot, jbc, jgl, kbl, lkahead, mvl, n2, n34,
362  $ n4, nbl, notrot, p, pskipped, q, rowskip, swband
363  LOGICAL applv, goscale, lower, lsvec, noscale, rotok,
364  $ rsvec, uctol, upper
365 * ..
366 * ..
367 * .. Intrinsic Functions ..
368  INTRINSIC abs, amax1, amin1, conjg, float, min0, max0,
369  $ sign, sqrt
370 * ..
371 * .. External Functions ..
372 * ..
373 * from BLAS
374  REAL scnrm2
375  COMPLEX cdotc
376  EXTERNAL cdotc, scnrm2
377  INTEGER isamax
378  EXTERNAL isamax
379 * from LAPACK
380  REAL slamch
381  EXTERNAL slamch
382  LOGICAL lsame
383  EXTERNAL lsame
384 * ..
385 * .. External Subroutines ..
386 * ..
387 * from BLAS
388  EXTERNAL ccopy, crot, csscal, cswap
389 * from LAPACK
390  EXTERNAL clascl, claset, classq, xerbla
391  EXTERNAL cgsvj0, cgsvj1
392 * ..
393 * .. Executable Statements ..
394 *
395 * Test the input arguments
396 *
397  lsvec = lsame( jobu, 'U' )
398  uctol = lsame( jobu, 'C' )
399  rsvec = lsame( jobv, 'V' )
400  applv = lsame( jobv, 'A' )
401  upper = lsame( joba, 'U' )
402  lower = lsame( joba, 'L' )
403 *
404  IF( .NOT.( upper .OR. lower .OR. lsame( joba, 'G' ) ) ) THEN
405  info = -1
406  ELSE IF( .NOT.( lsvec .OR. uctol .OR. lsame( jobu, 'N' ) ) ) THEN
407  info = -2
408  ELSE IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
409  info = -3
410  ELSE IF( m.LT.0 ) THEN
411  info = -4
412  ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
413  info = -5
414  ELSE IF( lda.LT.m ) THEN
415  info = -7
416  ELSE IF( mv.LT.0 ) THEN
417  info = -9
418  ELSE IF( ( rsvec .AND. ( ldv.LT.n ) ) .OR.
419  $ ( applv .AND. ( ldv.LT.mv ) ) ) THEN
420  info = -11
421  ELSE IF( uctol .AND. ( rwork( 1 ).LE.one ) ) THEN
422  info = -12
423  ELSE IF( lwork.LT.( m+n ) ) THEN
424  info = -13
425  ELSE IF( lrwork.LT.max0( n, 6 ) ) THEN
426  info = -15
427  ELSE
428  info = 0
429  END IF
430 *
431 * #:(
432  IF( info.NE.0 ) THEN
433  CALL xerbla( 'CGESVJ', -info )
434  RETURN
435  END IF
436 *
437 * #:) Quick return for void matrix
438 *
439  IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) )RETURN
440 *
441 * Set numerical parameters
442 * The stopping criterion for Jacobi rotations is
443 *
444 * max_{i<>j}|A(:,i)^* * A(:,j)| / (||A(:,i)||*||A(:,j)||) < CTOL*EPS
445 *
446 * where EPS is the round-off and CTOL is defined as follows:
447 *
448  IF( uctol ) THEN
449 * ... user controlled
450  ctol = rwork( 1 )
451  ELSE
452 * ... default
453  IF( lsvec .OR. rsvec .OR. applv ) THEN
454  ctol = sqrt( float( m ) )
455  ELSE
456  ctol = float( m )
457  END IF
458  END IF
459 * ... and the machine dependent parameters are
460 *[!] (Make sure that SLAMCH() works properly on the target machine.)
461 *
462  epsln = slamch( 'Epsilon' )
463  rooteps = sqrt( epsln )
464  sfmin = slamch( 'SafeMinimum' )
465  rootsfmin = sqrt( sfmin )
466  small = sfmin / epsln
467  big = slamch( 'Overflow' )
468 * BIG = ONE / SFMIN
469  rootbig = one / rootsfmin
470  large = big / sqrt( float( m*n ) )
471  bigtheta = one / rooteps
472 *
473  tol = ctol*epsln
474  roottol = sqrt( tol )
475 *
476  IF( float( m )*epsln.GE.one ) THEN
477  info = -4
478  CALL xerbla( 'CGESVJ', -info )
479  RETURN
480  END IF
481 *
482 * Initialize the right singular vector matrix.
483 *
484  IF( rsvec ) THEN
485  mvl = n
486  CALL claset( 'A', mvl, n, czero, cone, v, ldv )
487  ELSE IF( applv ) THEN
488  mvl = mv
489  END IF
490  rsvec = rsvec .OR. applv
491 *
492 * Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N )
493 *(!) If necessary, scale A to protect the largest singular value
494 * from overflow. It is possible that saving the largest singular
495 * value destroys the information about the small ones.
496 * This initial scaling is almost minimal in the sense that the
497 * goal is to make sure that no column norm overflows, and that
498 * SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
499 * in A are detected, the procedure returns with INFO=-6.
500 *
501  skl = one / sqrt( float( m )*float( n ) )
502  noscale = .true.
503  goscale = .true.
504 *
505  IF( lower ) THEN
506 * the input matrix is M-by-N lower triangular (trapezoidal)
507  DO 1874 p = 1, n
508  aapp = zero
509  aaqq = one
510  CALL classq( m-p+1, a( p, p ), 1, aapp, aaqq )
511  IF( aapp.GT.big ) THEN
512  info = -6
513  CALL xerbla( 'CGESVJ', -info )
514  RETURN
515  END IF
516  aaqq = sqrt( aaqq )
517  IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale ) THEN
518  sva( p ) = aapp*aaqq
519  ELSE
520  noscale = .false.
521  sva( p ) = aapp*( aaqq*skl )
522  IF( goscale ) THEN
523  goscale = .false.
524  DO 1873 q = 1, p - 1
525  sva( q ) = sva( q )*skl
526  1873 CONTINUE
527  END IF
528  END IF
529  1874 CONTINUE
530  ELSE IF( upper ) THEN
531 * the input matrix is M-by-N upper triangular (trapezoidal)
532  DO 2874 p = 1, n
533  aapp = zero
534  aaqq = one
535  CALL classq( p, a( 1, p ), 1, aapp, aaqq )
536  IF( aapp.GT.big ) THEN
537  info = -6
538  CALL xerbla( 'CGESVJ', -info )
539  RETURN
540  END IF
541  aaqq = sqrt( aaqq )
542  IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale ) THEN
543  sva( p ) = aapp*aaqq
544  ELSE
545  noscale = .false.
546  sva( p ) = aapp*( aaqq*skl )
547  IF( goscale ) THEN
548  goscale = .false.
549  DO 2873 q = 1, p - 1
550  sva( q ) = sva( q )*skl
551  2873 CONTINUE
552  END IF
553  END IF
554  2874 CONTINUE
555  ELSE
556 * the input matrix is M-by-N general dense
557  DO 3874 p = 1, n
558  aapp = zero
559  aaqq = one
560  CALL classq( m, a( 1, p ), 1, aapp, aaqq )
561  IF( aapp.GT.big ) THEN
562  info = -6
563  CALL xerbla( 'CGESVJ', -info )
564  RETURN
565  END IF
566  aaqq = sqrt( aaqq )
567  IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale ) THEN
568  sva( p ) = aapp*aaqq
569  ELSE
570  noscale = .false.
571  sva( p ) = aapp*( aaqq*skl )
572  IF( goscale ) THEN
573  goscale = .false.
574  DO 3873 q = 1, p - 1
575  sva( q ) = sva( q )*skl
576  3873 CONTINUE
577  END IF
578  END IF
579  3874 CONTINUE
580  END IF
581 *
582  IF( noscale )skl = one
583 *
584 * Move the smaller part of the spectrum from the underflow threshold
585 *(!) Start by determining the position of the nonzero entries of the
586 * array SVA() relative to ( SFMIN, BIG ).
587 *
588  aapp = zero
589  aaqq = big
590  DO 4781 p = 1, n
591  IF( sva( p ).NE.zero )aaqq = amin1( aaqq, sva( p ) )
592  aapp = amax1( aapp, sva( p ) )
593  4781 CONTINUE
594 *
595 * #:) Quick return for zero matrix
596 *
597  IF( aapp.EQ.zero ) THEN
598  IF( lsvec )CALL claset( 'G', m, n, czero, cone, a, lda )
599  rwork( 1 ) = one
600  rwork( 2 ) = zero
601  rwork( 3 ) = zero
602  rwork( 4 ) = zero
603  rwork( 5 ) = zero
604  rwork( 6 ) = zero
605  RETURN
606  END IF
607 *
608 * #:) Quick return for one-column matrix
609 *
610  IF( n.EQ.1 ) THEN
611  IF( lsvec )CALL clascl( 'G', 0, 0, sva( 1 ), skl, m, 1,
612  $ a( 1, 1 ), lda, ierr )
613  rwork( 1 ) = one / skl
614  IF( sva( 1 ).GE.sfmin ) THEN
615  rwork( 2 ) = one
616  ELSE
617  rwork( 2 ) = zero
618  END IF
619  rwork( 3 ) = zero
620  rwork( 4 ) = zero
621  rwork( 5 ) = zero
622  rwork( 6 ) = zero
623  RETURN
624  END IF
625 *
626 * Protect small singular values from underflow, and try to
627 * avoid underflows/overflows in computing Jacobi rotations.
628 *
629  sn = sqrt( sfmin / epsln )
630  temp1 = sqrt( big / float( n ) )
631  IF( ( aapp.LE.sn ) .OR. ( aaqq.GE.temp1 ) .OR.
632  $ ( ( sn.LE.aaqq ) .AND. ( aapp.LE.temp1 ) ) ) THEN
633  temp1 = amin1( big, temp1 / aapp )
634 * AAQQ = AAQQ*TEMP1
635 * AAPP = AAPP*TEMP1
636  ELSE IF( ( aaqq.LE.sn ) .AND. ( aapp.LE.temp1 ) ) THEN
637  temp1 = amin1( sn / aaqq, big / ( aapp*sqrt( float( n ) ) ) )
638 * AAQQ = AAQQ*TEMP1
639 * AAPP = AAPP*TEMP1
640  ELSE IF( ( aaqq.GE.sn ) .AND. ( aapp.GE.temp1 ) ) THEN
641  temp1 = amax1( sn / aaqq, temp1 / aapp )
642 * AAQQ = AAQQ*TEMP1
643 * AAPP = AAPP*TEMP1
644  ELSE IF( ( aaqq.LE.sn ) .AND. ( aapp.GE.temp1 ) ) THEN
645  temp1 = amin1( sn / aaqq, big / ( sqrt( float( n ) )*aapp ) )
646 * AAQQ = AAQQ*TEMP1
647 * AAPP = AAPP*TEMP1
648  ELSE
649  temp1 = one
650  END IF
651 *
652 * Scale, if necessary
653 *
654  IF( temp1.NE.one ) THEN
655  CALL slascl( 'G', 0, 0, one, temp1, n, 1, sva, n, ierr )
656  END IF
657  skl = temp1*skl
658  IF( skl.NE.one ) THEN
659  CALL clascl( joba, 0, 0, one, skl, m, n, a, lda, ierr )
660  skl = one / skl
661  END IF
662 *
663 * Row-cyclic Jacobi SVD algorithm with column pivoting
664 *
665  emptsw = ( n*( n-1 ) ) / 2
666  notrot = 0
667 
668  DO 1868 q = 1, n
669  cwork( q ) = cone
670  1868 CONTINUE
671 *
672 *
673 *
674  swband = 3
675 *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
676 * if CGESVJ is used as a computational routine in the preconditioned
677 * Jacobi SVD algorithm CGEJSV. For sweeps i=1:SWBAND the procedure
678 * works on pivots inside a band-like region around the diagonal.
679 * The boundaries are determined dynamically, based on the number of
680 * pivots above a threshold.
681 *
682  kbl = min0( 8, n )
683 *[TP] KBL is a tuning parameter that defines the tile size in the
684 * tiling of the p-q loops of pivot pairs. In general, an optimal
685 * value of KBL depends on the matrix dimensions and on the
686 * parameters of the computer's memory.
687 *
688  nbl = n / kbl
689  IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
690 *
691  blskip = kbl**2
692 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
693 *
694  rowskip = min0( 5, kbl )
695 *[TP] ROWSKIP is a tuning parameter.
696 *
697  lkahead = 1
698 *[TP] LKAHEAD is a tuning parameter.
699 *
700 * Quasi block transformations, using the lower (upper) triangular
701 * structure of the input matrix. The quasi-block-cycling usually
702 * invokes cubic convergence. Big part of this cycle is done inside
703 * canonical subspaces of dimensions less than M.
704 *
705  IF( ( lower .OR. upper ) .AND. ( n.GT.max0( 64, 4*kbl ) ) ) THEN
706 *[TP] The number of partition levels and the actual partition are
707 * tuning parameters.
708  n4 = n / 4
709  n2 = n / 2
710  n34 = 3*n4
711  IF( applv ) THEN
712  q = 0
713  ELSE
714  q = 1
715  END IF
716 *
717  IF( lower ) THEN
718 *
719 * This works very well on lower triangular matrices, in particular
720 * in the framework of the preconditioned Jacobi SVD (xGEJSV).
721 * The idea is simple:
722 * [+ 0 0 0] Note that Jacobi transformations of [0 0]
723 * [+ + 0 0] [0 0]
724 * [+ + x 0] actually work on [x 0] [x 0]
725 * [+ + x x] [x x]. [x x]
726 *
727  CALL cgsvj0( jobv, m-n34, n-n34, a( n34+1, n34+1 ), lda,
728  $ cwork( n34+1 ), sva( n34+1 ), mvl,
729  $ v( n34*q+1, n34+1 ), ldv, epsln, sfmin, tol,
730  $ 2, cwork( n+1 ), lwork-n, ierr )
731 
732  CALL cgsvj0( jobv, m-n2, n34-n2, a( n2+1, n2+1 ), lda,
733  $ cwork( n2+1 ), sva( n2+1 ), mvl,
734  $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 2,
735  $ cwork( n+1 ), lwork-n, ierr )
736 
737  CALL cgsvj1( jobv, m-n2, n-n2, n4, a( n2+1, n2+1 ), lda,
738  $ cwork( n2+1 ), sva( n2+1 ), mvl,
739  $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 1,
740  $ cwork( n+1 ), lwork-n, ierr )
741 *
742  CALL cgsvj0( jobv, m-n4, n2-n4, a( n4+1, n4+1 ), lda,
743  $ cwork( n4+1 ), sva( n4+1 ), mvl,
744  $ v( n4*q+1, n4+1 ), ldv, epsln, sfmin, tol, 1,
745  $ cwork( n+1 ), lwork-n, ierr )
746 *
747  CALL cgsvj0( jobv, m, n4, a, lda, cwork, sva, mvl, v, ldv,
748  $ epsln, sfmin, tol, 1, cwork( n+1 ), lwork-n,
749  $ ierr )
750 *
751  CALL cgsvj1( jobv, m, n2, n4, a, lda, cwork, sva, mvl, v,
752  $ ldv, epsln, sfmin, tol, 1, cwork( n+1 ),
753  $ lwork-n, ierr )
754 *
755 *
756  ELSE IF( upper ) THEN
757 *
758 *
759  CALL cgsvj0( jobv, n4, n4, a, lda, cwork, sva, mvl, v, ldv,
760  $ epsln, sfmin, tol, 2, cwork( n+1 ), lwork-n,
761  $ ierr )
762 *
763  CALL cgsvj0( jobv, n2, n4, a( 1, n4+1 ), lda, cwork( n4+1 ),
764  $ sva( n4+1 ), mvl, v( n4*q+1, n4+1 ), ldv,
765  $ epsln, sfmin, tol, 1, cwork( n+1 ), lwork-n,
766  $ ierr )
767 *
768  CALL cgsvj1( jobv, n2, n2, n4, a, lda, cwork, sva, mvl, v,
769  $ ldv, epsln, sfmin, tol, 1, cwork( n+1 ),
770  $ lwork-n, ierr )
771 *
772  CALL cgsvj0( jobv, n2+n4, n4, a( 1, n2+1 ), lda,
773  $ cwork( n2+1 ), sva( n2+1 ), mvl,
774  $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 1,
775  $ cwork( n+1 ), lwork-n, ierr )
776 
777  END IF
778 *
779  END IF
780 *
781 * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
782 *
783  DO 1993 i = 1, nsweep
784 *
785 * .. go go go ...
786 *
787  mxaapq = zero
788  mxsinj = zero
789  iswrot = 0
790 *
791  notrot = 0
792  pskipped = 0
793 *
794 * Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
795 * 1 <= p < q <= N. This is the first step toward a blocked implementation
796 * of the rotations. New implementation, based on block transformations,
797 * is under development.
798 *
799  DO 2000 ibr = 1, nbl
800 *
801  igl = ( ibr-1 )*kbl + 1
802 *
803  DO 1002 ir1 = 0, min0( lkahead, nbl-ibr )
804 *
805  igl = igl + ir1*kbl
806 *
807  DO 2001 p = igl, min0( igl+kbl-1, n-1 )
808 *
809 * .. de Rijk's pivoting
810 *
811  q = isamax( n-p+1, sva( p ), 1 ) + p - 1
812  IF( p.NE.q ) THEN
813  CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
814  IF( rsvec )CALL cswap( mvl, v( 1, p ), 1,
815  $ v( 1, q ), 1 )
816  temp1 = sva( p )
817  sva( p ) = sva( q )
818  sva( q ) = temp1
819  aapq = cwork(p)
820  cwork(p) = cwork(q)
821  cwork(q) = aapq
822  END IF
823 *
824  IF( ir1.EQ.0 ) THEN
825 *
826 * Column norms are periodically updated by explicit
827 * norm computation.
828 *[!] Caveat:
829 * Unfortunately, some BLAS implementations compute SCNRM2(M,A(1,p),1)
830 * as SQRT(S=CDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to
831 * overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
832 * underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
833 * Hence, SCNRM2 cannot be trusted, not even in the case when
834 * the true norm is far from the under(over)flow boundaries.
835 * If properly implemented SCNRM2 is available, the IF-THEN-ELSE-END IF
836 * below should be replaced with "AAPP = SCNRM2( M, A(1,p), 1 )".
837 *
838  IF( ( sva( p ).LT.rootbig ) .AND.
839  $ ( sva( p ).GT.rootsfmin ) ) THEN
840  sva( p ) = scnrm2( m, a( 1, p ), 1 )
841  ELSE
842  temp1 = zero
843  aapp = one
844  CALL classq( m, a( 1, p ), 1, temp1, aapp )
845  sva( p ) = temp1*sqrt( aapp )
846  END IF
847  aapp = sva( p )
848  ELSE
849  aapp = sva( p )
850  END IF
851 *
852  IF( aapp.GT.zero ) THEN
853 *
854  pskipped = 0
855 *
856  DO 2002 q = p + 1, min0( igl+kbl-1, n )
857 *
858  aaqq = sva( q )
859 *
860  IF( aaqq.GT.zero ) THEN
861 *
862  aapp0 = aapp
863  IF( aaqq.GE.one ) THEN
864  rotok = ( small*aapp ).LE.aaqq
865  IF( aapp.LT.( big / aaqq ) ) THEN
866  aapq = ( cdotc( m, a( 1, p ), 1,
867  $ a( 1, q ), 1 ) / aaqq ) / aapp
868  ELSE
869  CALL ccopy( m, a( 1, p ), 1,
870  $ cwork(n+1), 1 )
871  CALL clascl( 'G', 0, 0, aapp, one,
872  $ m, 1, cwork(n+1), lda, ierr )
873  aapq = cdotc( m, cwork(n+1), 1,
874  $ a( 1, q ), 1 ) / aaqq
875  END IF
876  ELSE
877  rotok = aapp.LE.( aaqq / small )
878  IF( aapp.GT.( small / aaqq ) ) THEN
879  aapq = ( cdotc( m, a( 1, p ), 1,
880  $ a( 1, q ), 1 ) / aaqq ) / aapp
881  ELSE
882  CALL ccopy( m, a( 1, q ), 1,
883  $ cwork(n+1), 1 )
884  CALL clascl( 'G', 0, 0, aaqq,
885  $ one, m, 1,
886  $ cwork(n+1), lda, ierr )
887  aapq = cdotc( m, a(1, p ), 1,
888  $ cwork(n+1), 1 ) / aapp
889  END IF
890  END IF
891 *
892  ompq = aapq / abs(aapq)
893 * AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
894  aapq1 = -abs(aapq)
895  mxaapq = amax1( mxaapq, -aapq1 )
896 *
897 * TO rotate or NOT to rotate, THAT is the question ...
898 *
899  IF( abs( aapq1 ).GT.tol ) THEN
900 *
901 * .. rotate
902 *[RTD] ROTATED = ROTATED + ONE
903 *
904  IF( ir1.EQ.0 ) THEN
905  notrot = 0
906  pskipped = 0
907  iswrot = iswrot + 1
908  END IF
909 *
910  IF( rotok ) THEN
911 *
912  aqoap = aaqq / aapp
913  apoaq = aapp / aaqq
914  theta = -half*abs( aqoap-apoaq )/aapq1
915 *
916  IF( abs( theta ).GT.bigtheta ) THEN
917 *
918  t = half / theta
919  cs = one
920 
921  CALL crot( m, a(1,p), 1, a(1,q), 1,
922  $ cs, conjg(ompq)*t )
923  IF ( rsvec ) THEN
924  CALL crot( mvl, v(1,p), 1,
925  $ v(1,q), 1, cs, conjg(ompq)*t )
926  END IF
927 
928  sva( q ) = aaqq*sqrt( amax1( zero,
929  $ one+t*apoaq*aapq1 ) )
930  aapp = aapp*sqrt( amax1( zero,
931  $ one-t*aqoap*aapq1 ) )
932  mxsinj = amax1( mxsinj, abs( t ) )
933 *
934  ELSE
935 *
936 * .. choose correct signum for THETA and rotate
937 *
938  thsign = -sign( one, aapq1 )
939  t = one / ( theta+thsign*
940  $ sqrt( one+theta*theta ) )
941  cs = sqrt( one / ( one+t*t ) )
942  sn = t*cs
943 *
944  mxsinj = amax1( mxsinj, abs( sn ) )
945  sva( q ) = aaqq*sqrt( amax1( zero,
946  $ one+t*apoaq*aapq1 ) )
947  aapp = aapp*sqrt( amax1( zero,
948  $ one-t*aqoap*aapq1 ) )
949 *
950  CALL crot( m, a(1,p), 1, a(1,q), 1,
951  $ cs, conjg(ompq)*sn )
952  IF ( rsvec ) THEN
953  CALL crot( mvl, v(1,p), 1,
954  $ v(1,q), 1, cs, conjg(ompq)*sn )
955  END IF
956  END IF
957  cwork(p) = -cwork(q) * ompq
958 *
959  ELSE
960 * .. have to use modified Gram-Schmidt like transformation
961  CALL ccopy( m, a( 1, p ), 1,
962  $ cwork(n+1), 1 )
963  CALL clascl( 'G', 0, 0, aapp, one, m,
964  $ 1, cwork(n+1), lda,
965  $ ierr )
966  CALL clascl( 'G', 0, 0, aaqq, one, m,
967  $ 1, a( 1, q ), lda, ierr )
968  CALL caxpy( m, -aapq, cwork(n+1), 1,
969  $ a( 1, q ), 1 )
970  CALL clascl( 'G', 0, 0, one, aaqq, m,
971  $ 1, a( 1, q ), lda, ierr )
972  sva( q ) = aaqq*sqrt( amax1( zero,
973  $ one-aapq1*aapq1 ) )
974  mxsinj = amax1( mxsinj, sfmin )
975  END IF
976 * END IF ROTOK THEN ... ELSE
977 *
978 * In the case of cancellation in updating SVA(q), SVA(p)
979 * recompute SVA(q), SVA(p).
980 *
981  IF( ( sva( q ) / aaqq )**2.LE.rooteps )
982  $ THEN
983  IF( ( aaqq.LT.rootbig ) .AND.
984  $ ( aaqq.GT.rootsfmin ) ) THEN
985  sva( q ) = scnrm2( m, a( 1, q ), 1 )
986  ELSE
987  t = zero
988  aaqq = one
989  CALL classq( m, a( 1, q ), 1, t,
990  $ aaqq )
991  sva( q ) = t*sqrt( aaqq )
992  END IF
993  END IF
994  IF( ( aapp / aapp0 ).LE.rooteps ) THEN
995  IF( ( aapp.LT.rootbig ) .AND.
996  $ ( aapp.GT.rootsfmin ) ) THEN
997  aapp = scnrm2( m, a( 1, p ), 1 )
998  ELSE
999  t = zero
1000  aapp = one
1001  CALL classq( m, a( 1, p ), 1, t,
1002  $ aapp )
1003  aapp = t*sqrt( aapp )
1004  END IF
1005  sva( p ) = aapp
1006  END IF
1007 *
1008  ELSE
1009 * A(:,p) and A(:,q) already numerically orthogonal
1010  IF( ir1.EQ.0 )notrot = notrot + 1
1011 *[RTD] SKIPPED = SKIPPED + 1
1012  pskipped = pskipped + 1
1013  END IF
1014  ELSE
1015 * A(:,q) is zero column
1016  IF( ir1.EQ.0 )notrot = notrot + 1
1017  pskipped = pskipped + 1
1018  END IF
1019 *
1020  IF( ( i.LE.swband ) .AND.
1021  $ ( pskipped.GT.rowskip ) ) THEN
1022  IF( ir1.EQ.0 )aapp = -aapp
1023  notrot = 0
1024  GO TO 2103
1025  END IF
1026 *
1027  2002 CONTINUE
1028 * END q-LOOP
1029 *
1030  2103 CONTINUE
1031 * bailed out of q-loop
1032 *
1033  sva( p ) = aapp
1034 *
1035  ELSE
1036  sva( p ) = aapp
1037  IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
1038  $ notrot = notrot + min0( igl+kbl-1, n ) - p
1039  END IF
1040 *
1041  2001 CONTINUE
1042 * end of the p-loop
1043 * end of doing the block ( ibr, ibr )
1044  1002 CONTINUE
1045 * end of ir1-loop
1046 *
1047 * ... go to the off diagonal blocks
1048 *
1049  igl = ( ibr-1 )*kbl + 1
1050 *
1051  DO 2010 jbc = ibr + 1, nbl
1052 *
1053  jgl = ( jbc-1 )*kbl + 1
1054 *
1055 * doing the block at ( ibr, jbc )
1056 *
1057  ijblsk = 0
1058  DO 2100 p = igl, min0( igl+kbl-1, n )
1059 *
1060  aapp = sva( p )
1061  IF( aapp.GT.zero ) THEN
1062 *
1063  pskipped = 0
1064 *
1065  DO 2200 q = jgl, min0( jgl+kbl-1, n )
1066 *
1067  aaqq = sva( q )
1068  IF( aaqq.GT.zero ) THEN
1069  aapp0 = aapp
1070 *
1071 * .. M x 2 Jacobi SVD ..
1072 *
1073 * Safe Gram matrix computation
1074 *
1075  IF( aaqq.GE.one ) THEN
1076  IF( aapp.GE.aaqq ) THEN
1077  rotok = ( small*aapp ).LE.aaqq
1078  ELSE
1079  rotok = ( small*aaqq ).LE.aapp
1080  END IF
1081  IF( aapp.LT.( big / aaqq ) ) THEN
1082  aapq = ( cdotc( m, a( 1, p ), 1,
1083  $ a( 1, q ), 1 ) / aaqq ) / aapp
1084  ELSE
1085  CALL ccopy( m, a( 1, p ), 1,
1086  $ cwork(n+1), 1 )
1087  CALL clascl( 'G', 0, 0, aapp,
1088  $ one, m, 1,
1089  $ cwork(n+1), lda, ierr )
1090  aapq = cdotc( m, cwork(n+1), 1,
1091  $ a( 1, q ), 1 ) / aaqq
1092  END IF
1093  ELSE
1094  IF( aapp.GE.aaqq ) THEN
1095  rotok = aapp.LE.( aaqq / small )
1096  ELSE
1097  rotok = aaqq.LE.( aapp / small )
1098  END IF
1099  IF( aapp.GT.( small / aaqq ) ) THEN
1100  aapq = ( cdotc( m, a( 1, p ), 1,
1101  $ a( 1, q ), 1 ) / aaqq ) / aapp
1102  ELSE
1103  CALL ccopy( m, a( 1, q ), 1,
1104  $ cwork(n+1), 1 )
1105  CALL clascl( 'G', 0, 0, aaqq,
1106  $ one, m, 1,
1107  $ cwork(n+1), lda, ierr )
1108  aapq = cdotc( m, a( 1, p ), 1,
1109  $ cwork(n+1), 1 ) / aapp
1110  END IF
1111  END IF
1112 *
1113  ompq = aapq / abs(aapq)
1114 * AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
1115  aapq1 = -abs(aapq)
1116  mxaapq = a