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

Functions

subroutine dgeqpf (M, N, A, LDA, JPVT, TAU, WORK, INFO)
 DGEQPF More...
 
subroutine dgebak (JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
 DGEBAK More...
 
subroutine dgebal (JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
 DGEBAL More...
 
subroutine dgebd2 (M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO)
 DGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm. More...
 
subroutine dgebrd (M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
 DGEBRD More...
 
subroutine dgecon (NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
 DGECON More...
 
subroutine dgeequ (M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
 DGEEQU More...
 
subroutine dgeequb (M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
 DGEEQUB More...
 
subroutine dgehd2 (N, ILO, IHI, A, LDA, TAU, WORK, INFO)
 DGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm. More...
 
subroutine dgehrd (N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
 DGEHRD More...
 
subroutine dgelq2 (M, N, A, LDA, TAU, WORK, INFO)
 DGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine dgelqf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 DGELQF More...
 
subroutine dgemqrt (SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
 DGEMQRT More...
 
subroutine dgeql2 (M, N, A, LDA, TAU, WORK, INFO)
 DGEQL2 computes the QL factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine dgeqlf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 DGEQLF More...
 
subroutine dgeqp3 (M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
 DGEQP3 More...
 
subroutine dgeqr2 (M, N, A, LDA, TAU, WORK, INFO)
 DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine dgeqr2p (M, N, A, LDA, TAU, WORK, INFO)
 DGEQR2P computes the QR factorization of a general rectangular matrix with non-negative diagonal elements using an unblocked algorithm. More...
 
subroutine dgeqrf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 DGEQRF More...
 
subroutine dgeqrfp (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 DGEQRFP More...
 
subroutine dgeqrt (M, N, NB, A, LDA, T, LDT, WORK, INFO)
 DGEQRT More...
 
subroutine dgeqrt2 (M, N, A, LDA, T, LDT, INFO)
 DGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY representation of Q. More...
 
recursive subroutine dgeqrt3 (M, N, A, LDA, T, LDT, INFO)
 DGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact WY representation of Q. More...
 
subroutine dgerfs (TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
 DGERFS More...
 
subroutine dgerfsx (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, IWORK, INFO)
 DGERFSX More...
 
subroutine dgerq2 (M, N, A, LDA, TAU, WORK, INFO)
 DGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine dgerqf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 DGERQF More...
 
subroutine dgesvj (JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
 DGESVJ More...
 
subroutine dgetf2 (M, N, A, LDA, IPIV, INFO)
 DGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row interchanges (unblocked algorithm). More...
 
subroutine dgetrf (M, N, A, LDA, IPIV, INFO)
 DGETRF More...
 
recursive subroutine dgetrf2 (M, N, A, LDA, IPIV, INFO)
 DGETRF2 More...
 
subroutine dgetri (N, A, LDA, IPIV, WORK, LWORK, INFO)
 DGETRI More...
 
subroutine dgetrs (TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 DGETRS More...
 
subroutine dhgeqz (JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
 DHGEQZ More...
 
subroutine dla_geamv (TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
 DLA_GEAMV computes a matrix-vector product using a general matrix to calculate error bounds. More...
 
double precision function dla_gercond (TRANS, N, A, LDA, AF, LDAF, IPIV, CMODE, C, INFO, WORK, IWORK)
 DLA_GERCOND estimates the Skeel condition number for a general matrix. More...
 
subroutine dla_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)
 DLA_GERFSX_EXTENDED improves the computed solution to a system of linear equations for general matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution. More...
 
double precision function dla_gerpvgrw (N, NCOLS, A, LDA, AF, LDAF)
 DLA_GERPVGRW More...
 
subroutine dtgevc (SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
 DTGEVC More...
 
subroutine dtgexc (WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
 DTGEXC More...
 
subroutine zgesvj (JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, CWORK, LWORK, RWORK, LRWORK, INFO)
 ZGESVJ More...
 

Detailed Description

This is the group of double computational functions for GE matrices

Function Documentation

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

DGEBAK

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

Purpose:
 DGEBAK forms the right or left eigenvectors of a real general matrix
 by backward transformation on the computed eigenvectors of the
 balanced matrix output by DGEBAL.
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 DGEBAL.
[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 DGEBAL.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
[in]SCALE
          SCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutation and scaling factors, as returned
          by DGEBAL.
[in]M
          M is INTEGER
          The number of columns of the matrix V.  M >= 0.
[in,out]V
          V is DOUBLE PRECISION array, dimension (LDV,M)
          On entry, the matrix of right or left eigenvectors to be
          transformed, as returned by DHSEIN or DTREVC.
          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 132 of file dgebak.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

DGEBAL

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

Purpose:
 DGEBAL balances a general real 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 DOUBLE 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 DOUBLE 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 BALANC.

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

Definition at line 162 of file dgebal.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 DGEBD2 reduces a real general m by n matrix A to upper or lower
 bidiagonal form B by an orthogonal transformation: Q**T * 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 DOUBLE PRECISION 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 orthogonal matrix Q as a product of elementary
            reflectors, and the elements above the first superdiagonal,
            with the array TAUP, represent the orthogonal 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 orthogonal matrix Q as a product of
            elementary reflectors, and the elements above the diagonal,
            with the array TAUP, represent the orthogonal matrix P as
            a product of elementary reflectors.
          See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]D
          D is DOUBLE PRECISION array, dimension (min(M,N))
          The diagonal elements of the bidiagonal matrix B:
          D(i) = A(i,i).
[out]E
          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
          The off-diagonal elements of the bidiagonal matrix B:
          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
[out]TAUQ
          TAUQ is DOUBLE PRECISION array dimension (min(M,N))
          The scalar factors of the elementary reflectors which
          represent the orthogonal matrix Q. See Further Details.
[out]TAUP
          TAUP is DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors which
          represent the orthogonal matrix P. See Further Details.
[out]WORK
          WORK is DOUBLE PRECISION 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**T  and G(i) = I - taup * u * u**T

  where tauq and taup are real scalars, and v and u are real 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**T  and G(i) = I - taup * u * u**T

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

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

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

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

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

Definition at line 191 of file dgebd2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

DGEBRD

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

Purpose:
 DGEBRD reduces a general real M-by-N matrix A to upper or lower
 bidiagonal form B by an orthogonal transformation: Q**T * 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 DOUBLE PRECISION 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 orthogonal matrix Q as a product of elementary
            reflectors, and the elements above the first superdiagonal,
            with the array TAUP, represent the orthogonal 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 orthogonal matrix Q as a product of
            elementary reflectors, and the elements above the diagonal,
            with the array TAUP, represent the orthogonal matrix P as
            a product of elementary reflectors.
          See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]D
          D is DOUBLE PRECISION array, dimension (min(M,N))
          The diagonal elements of the bidiagonal matrix B:
          D(i) = A(i,i).
[out]E
          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
          The off-diagonal elements of the bidiagonal matrix B:
          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
[out]TAUQ
          TAUQ is DOUBLE PRECISION array dimension (min(M,N))
          The scalar factors of the elementary reflectors which
          represent the orthogonal matrix Q. See Further Details.
[out]TAUP
          TAUP is DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors which
          represent the orthogonal matrix P. See Further Details.
[out]WORK
          WORK is DOUBLE PRECISION 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**T  and G(i) = I - taup * u * u**T

  where tauq and taup are real scalars, and v and u are real 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**T  and G(i) = I - taup * u * u**T

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

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

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

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

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

Definition at line 207 of file dgebrd.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgecon ( character  NORM,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision  ANORM,
double precision  RCOND,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DGECON

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

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

 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 DOUBLE PRECISION array, dimension (LDA,N)
          The factors L and U from the factorization A = P*L*U
          as computed by DGETRF.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]ANORM
          ANORM is DOUBLE PRECISION
          If NORM = '1' or 'O', the 1-norm of the original matrix A.
          If NORM = 'I', the infinity-norm of the original matrix A.
[out]RCOND
          RCOND is DOUBLE PRECISION
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(norm(A) * norm(inv(A))).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (4*N)
[out]IWORK
          IWORK is INTEGER 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 2011

Definition at line 126 of file dgecon.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

DGEEQU

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

Purpose:
 DGEEQU 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 DOUBLE PRECISION array, dimension (LDA,N)
          The M-by-N matrix whose equilibration factors are
          to be computed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]R
          R is DOUBLE PRECISION array, dimension (M)
          If INFO = 0 or INFO > M, R contains the row scale factors
          for A.
[out]C
          C is DOUBLE PRECISION array, dimension (N)
          If INFO = 0,  C contains the column scale factors for A.
[out]ROWCND
          ROWCND is DOUBLE PRECISION
          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
          AMAX is neither too large nor too small, it is not worth
          scaling by R.
[out]COLCND
          COLCND is DOUBLE PRECISION
          If INFO = 0, COLCND contains the ratio of the smallest
          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
          worth scaling by C.
[out]AMAX
          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i,  and i is
                <= M:  the i-th row of A is exactly zero
                >  M:  the (i-M)-th column of A is exactly zero
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 141 of file dgeequ.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

DGEEQUB

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

Purpose:
 DGEEQUB 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 DGEEQU 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 DOUBLE PRECISION array, dimension (LDA,N)
          The M-by-N matrix whose equilibration factors are
          to be computed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]R
          R is DOUBLE PRECISION array, dimension (M)
          If INFO = 0 or INFO > M, R contains the row scale factors
          for A.
[out]C
          C is DOUBLE PRECISION array, dimension (N)
          If INFO = 0,  C contains the column scale factors for A.
[out]ROWCND
          ROWCND is DOUBLE PRECISION
          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
          AMAX is neither too large nor too small, it is not worth
          scaling by R.
[out]COLCND
          COLCND is DOUBLE PRECISION
          If INFO = 0, COLCND contains the ratio of the smallest
          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
          worth scaling by C.
[out]AMAX
          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i,  and i is
                <= M:  the i-th row of A is exactly zero
                >  M:  the (i-M)-th column of A is exactly zero
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 148 of file dgeequb.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgehd2 ( integer  N,
integer  ILO,
integer  IHI,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 DGEHD2 reduces a real general matrix A to upper Hessenberg form H by
 an orthogonal similarity transformation:  Q**T * 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 DGEBAL; otherwise they should be
          set to 1 and N respectively. See Further Details.
          1 <= ILO <= IHI <= max(1,N).
[in,out]A
          A is DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION array, dimension (N-1)
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
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**T

  where tau is a real scalar, and v is a real 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 dgehd2.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  DOUBLE PRECISION a( lda, * ), tau( * ), work( * )
162 * ..
163 *
164 * =====================================================================
165 *
166 * .. Parameters ..
167  DOUBLE PRECISION one
168  parameter( one = 1.0d+0 )
169 * ..
170 * .. Local Scalars ..
171  INTEGER i
172  DOUBLE PRECISION aii
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL dlarf, dlarfg, xerbla
176 * ..
177 * .. Intrinsic Functions ..
178  INTRINSIC 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( 'DGEHD2', -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  CALL dlarfg( ihi-i, a( i+1, i ), a( min( i+2, n ), i ), 1,
204  $ tau( i ) )
205  aii = a( i+1, i )
206  a( i+1, i ) = one
207 *
208 * Apply H(i) to A(1:ihi,i+1:ihi) from the right
209 *
210  CALL dlarf( 'Right', ihi, ihi-i, a( i+1, i ), 1, tau( i ),
211  $ a( 1, i+1 ), lda, work )
212 *
213 * Apply H(i) to A(i+1:ihi,i+1:n) from the left
214 *
215  CALL dlarf( 'Left', ihi-i, n-i, a( i+1, i ), 1, tau( i ),
216  $ a( i+1, i+1 ), lda, work )
217 *
218  a( i+1, i ) = aii
219  10 CONTINUE
220 *
221  RETURN
222 *
223 * End of DGEHD2
224 *
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:126

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgehrd ( integer  N,
integer  ILO,
integer  IHI,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGEHRD

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

Purpose:
 DGEHRD reduces a real general matrix A to upper Hessenberg form H by
 an orthogonal similarity transformation:  Q**T * 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 DGEBAL; 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 DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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**T

  where tau is a real scalar, and v is a real 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 dgehrd.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  DOUBLE PRECISION 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  DOUBLE PRECISION zero, one
189  parameter( zero = 0.0d+0,
190  $ one = 1.0d+0 )
191 * ..
192 * .. Local Scalars ..
193  LOGICAL lquery
194  INTEGER i, ib, iinfo, iwt, j, ldwork, lwkopt, nb,
195  $ nbmin, nh, nx
196  DOUBLE PRECISION ei
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL daxpy, dgehd2, dgemm, dlahr2, dlarfb, dtrmm,
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, 'DGEHRD', ' ', 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( 'DGEHRD', -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, 'DGEHRD', ' ', 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, 'DGEHRD', ' ', 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, 'DGEHRD', ' ', 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**T
308 * which performs the reduction, and also the matrix Y = A*V*T
309 *
310  CALL dlahr2( 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**T. 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 dgemm( 'No transpose', '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 dtrmm( 'Right', 'Lower', 'Transpose',
329  $ 'Unit', i, ib-1,
330  $ one, a( i+1, i ), lda, work, ldwork )
331  DO 30 j = 0, ib-2
332  CALL daxpy( 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 dlarfb( 'Left', '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 dgehd2( n, i, ihi, a, lda, tau, work, iinfo )
350  work( 1 ) = lwkopt
351 *
352  RETURN
353 *
354 * End of DGEHRD
355 *
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM
Definition: dtrmm.f:179
subroutine dgehd2(N, ILO, IHI, A, LDA, TAU, WORK, INFO)
DGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm...
Definition: dgehd2.f:151
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: dlarfb.f:197
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dlahr2(N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)
DLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elemen...
Definition: dlahr2.f:183
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgelq2 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 DGELQ2 computes an LQ factorization of a real 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 DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION 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(2) H(1), where k = min(m,n).

  Each H(i) has the form

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

  where tau is a real scalar, and v is a real vector with
  v(1:i-1) = 0 and v(i) = 1; 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 dgelq2.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  DOUBLE PRECISION a( lda, * ), tau( * ), work( * )
134 * ..
135 *
136 * =====================================================================
137 *
138 * .. Parameters ..
139  DOUBLE PRECISION one
140  parameter( one = 1.0d+0 )
141 * ..
142 * .. Local Scalars ..
143  INTEGER i, k
144  DOUBLE PRECISION aii
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL dlarf, dlarfg, 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( 'DGELQ2', -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 dlarfg( n-i+1, a( i, i ), a( i, min( i+1, n ) ), lda,
176  $ tau( i ) )
177  IF( i.LT.m ) THEN
178 *
179 * Apply H(i) to A(i+1:m,i:n) from the right
180 *
181  aii = a( i, i )
182  a( i, i ) = one
183  CALL dlarf( 'Right', m-i, n-i+1, a( i, i ), lda, tau( i ),
184  $ a( i+1, i ), lda, work )
185  a( i, i ) = aii
186  END IF
187  10 CONTINUE
188  RETURN
189 *
190 * End of DGELQ2
191 *
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:126

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgelqf ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGELQF

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

Purpose:
 DGELQF computes an LQ factorization of a real 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 DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION 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(2) H(1), where k = min(m,n).

  Each H(i) has the form

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

  where tau is a real scalar, and v is a real vector with
  v(1:i-1) = 0 and v(i) = 1; 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 dgelqf.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  DOUBLE PRECISION 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 dgelq2, dlarfb, dlarft, 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, 'DGELQF', ' ', 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( 'DGELQF', -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, 'DGELQF', ' ', 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, 'DGELQF', ' ', 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 dgelq2( 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 dlarft( '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 dlarfb( '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 dgelq2( 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 DGELQF
268 *
subroutine dgelq2(M, N, A, LDA, TAU, WORK, INFO)
DGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: dgelq2.f:123
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
DLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: dlarft.f:165
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: dlarfb.f:197
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 dgemqrt ( character  SIDE,
character  TRANS,
integer  M,
integer  N,
integer  K,
integer  NB,
double precision, dimension( ldv, * )  V,
integer  LDV,
double precision, dimension( ldt, * )  T,
integer  LDT,
double precision, dimension( ldc, * )  C,
integer  LDC,
double precision, dimension( * )  WORK,
integer  INFO 
)

DGEMQRT

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

Purpose:
 DGEMQRT overwrites the general real M-by-N matrix C with

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

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

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

 generated using the compact WY representation as returned by DGEQRT. 

 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**T from the Left;
          = 'R': apply Q or Q**T from the Right.
[in]TRANS
          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'C':  Transpose, apply Q**T.
[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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q C, Q**T C, C Q**T or C Q.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).
[out]WORK
          WORK is DOUBLE PRECISION 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 dgemqrt.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  DOUBLE PRECISION 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, dlarfb
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, 'T' )
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( 'DGEMQRT', -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 dlarfb( 'L', 'T', '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 dlarfb( '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 dlarfb( '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 dlarfb( 'R', 'T', '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 DGEMQRT
290 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: dlarfb.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 dgeql2 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 DGEQL2 computes a QL factorization of a real 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 DOUBLE PRECISION 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
          orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
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**T

  where tau is a real scalar, and v is a real 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 dgeql2.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  DOUBLE PRECISION a( lda, * ), tau( * ), work( * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  DOUBLE PRECISION one
142  parameter( one = 1.0d+0 )
143 * ..
144 * .. Local Scalars ..
145  INTEGER i, k
146  DOUBLE PRECISION aii
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL dlarf, dlarfg, 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( 'DGEQL2', -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  CALL dlarfg( m-k+i, a( m-k+i, n-k+i ), a( 1, n-k+i ), 1,
179  $ tau( i ) )
180 *
181 * Apply H(i) to A(1:m-k+i,1:n-k+i-1) from the left
182 *
183  aii = a( m-k+i, n-k+i )
184  a( m-k+i, n-k+i ) = one
185  CALL dlarf( 'Left', m-k+i, n-k+i-1, a( 1, n-k+i ), 1, tau( i ),
186  $ a, lda, work )
187  a( m-k+i, n-k+i ) = aii
188  10 CONTINUE
189  RETURN
190 *
191 * End of DGEQL2
192 *
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:126

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgeqlf ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGEQLF

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

Purpose:
 DGEQLF computes a QL factorization of a real 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 DOUBLE PRECISION 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
          orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION 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**T

  where tau is a real scalar, and v is a real 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 dgeqlf.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  DOUBLE PRECISION 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 dgeql2, dlarfb, dlarft, 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, 'DGEQLF', ' ', 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( 'DGEQLF', -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, 'DGEQLF', ' ', 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, 'DGEQLF', ' ', 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 dgeql2( 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 dlarft( 'Backward', 'Columnwise', m-k+i+ib-1, ib,
260  $ a( 1, n-k+i ), lda, tau( i ), work, ldwork )
261 *
262 * Apply H**T to A(1:m-k+i+ib-1,1:n-k+i-1) from the left
263 *
264  CALL dlarfb( 'Left', '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 dgeql2( mu, nu, a, lda, tau, work, iinfo )
281 *
282  work( 1 ) = iws
283  RETURN
284 *
285 * End of DGEQLF
286 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
DLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: dlarft.f:165
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: dlarfb.f:197
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine dgeql2(M, N, A, LDA, TAU, WORK, INFO)
DGEQL2 computes the QL factorization of a general rectangular matrix using an unblocked algorithm...
Definition: dgeql2.f:125

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgeqp3 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGEQP3

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

Purpose:
 DGEQP3 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 DOUBLE PRECISION 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
          orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[out]WORK
          WORK is DOUBLE PRECISION 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 >= 3*N+1.
          For optimal performance LWORK >= 2*N+( 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]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**T

  where tau is a real 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 153 of file dgeqp3.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgeqpf ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  INFO 
)

DGEQPF

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

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

 DGEQPF computes a QR factorization with column pivoting of a
 real 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 DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (3*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**T

  where tau is a real scalar, and v is a real 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 144 of file dgeqpf.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgeqr2 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 DGEQR2 computes a QR factorization of a real 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 DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
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**T

  where tau is a real scalar, and v is a real 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 dgeqr2.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  DOUBLE PRECISION a( lda, * ), tau( * ), work( * )
134 * ..
135 *
136 * =====================================================================
137 *
138 * .. Parameters ..
139  DOUBLE PRECISION one
140  parameter( one = 1.0d+0 )
141 * ..
142 * .. Local Scalars ..
143  INTEGER i, k
144  DOUBLE PRECISION aii
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL dlarf, dlarfg, 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( 'DGEQR2', -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 dlarfg( 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) to A(i:m,i+1:n) from the left
180 *
181  aii = a( i, i )
182  a( i, i ) = one
183  CALL dlarf( 'Left', m-i+1, n-i, a( i, i ), 1, tau( i ),
184  $ a( i, i+1 ), lda, work )
185  a( i, i ) = aii
186  END IF
187  10 CONTINUE
188  RETURN
189 *
190 * End of DGEQR2
191 *
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:126

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgeqr2p ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 DGEQR2 computes a QR factorization of a real m by n matrix A:
 A = Q * R. The diagonal entries of R are 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 DOUBLE PRECISION 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
          nonnegative; the elements below the diagonal,
          with the array TAU, represent the orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
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**T

  where tau is a real scalar, and v is a real 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 dgeqr2p.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  DOUBLE PRECISION a( lda, * ), tau( * ), work( * )
137 * ..
138 *
139 * =====================================================================
140 *
141 * .. Parameters ..
142  DOUBLE PRECISION one
143  parameter( one = 1.0d+0 )
144 * ..
145 * .. Local Scalars ..
146  INTEGER i, k
147  DOUBLE PRECISION aii
148 * ..
149 * .. External Subroutines ..
150  EXTERNAL dlarf, dlarfgp, xerbla
151 * ..
152 * .. Intrinsic Functions ..
153  INTRINSIC 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( 'DGEQR2P', -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 dlarfgp( 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) to A(i:m,i+1:n) from the left
183 *
184  aii = a( i, i )
185  a( i, i ) = one
186  CALL dlarf( 'Left', m-i+1, n-i, a( i, i ), 1, tau( i ),
187  $ a( i, i+1 ), lda, work )
188  a( i, i ) = aii
189  END IF
190  10 CONTINUE
191  RETURN
192 *
193 * End of DGEQR2P
194 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarfgp(N, ALPHA, X, INCX, TAU)
DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: dlarfgp.f:106
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:126

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgeqrf ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGEQRF

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

Purpose:
 DGEQRF computes a QR factorization of a real 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 DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION 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**T

  where tau is a real scalar, and v is a real 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 dgeqrf.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  DOUBLE PRECISION 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 dgeqr2, dlarfb, dlarft, 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, 'DGEQRF', ' ', 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( 'DGEQRF', -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, 'DGEQRF', ' ', 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, 'DGEQRF', ' ', 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 dgeqr2( 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 dlarft( 'Forward', 'Columnwise', m-i+1, ib,
245  $ a( i, i ), lda, tau( i ), work, ldwork )
246 *
247 * Apply H**T to A(i:m,i+ib:n) from the left
248 *
249  CALL dlarfb( 'Left', '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 dgeqr2( 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 DGEQRF
269 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
DLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: dlarft.f:165
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: dlarfb.f:197
subroutine dgeqr2(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: dgeqr2.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 dgeqrfp ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGEQRFP

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

Purpose:
 DGEQRFP computes a QR factorization of a real M-by-N matrix A:
 A = Q * R. The diagonal entries of R are 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 DOUBLE PRECISION 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 nonnegative; the elements below the diagonal,
          with the array TAU, represent the orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION 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**T

  where tau is a real scalar, and v is a real 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 dgeqrfp.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  DOUBLE PRECISION 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 dgeqr2p, dlarfb, dlarft, 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, 'DGEQRF', ' ', 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( 'DGEQRFP', -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, 'DGEQRF', ' ', 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, 'DGEQRF', ' ', 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 dgeqr2p( 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 dlarft( 'Forward', 'Columnwise', m-i+1, ib,
248  $ a( i, i ), lda, tau( i ), work, ldwork )
249 *
250 * Apply H**T to A(i:m,i+ib:n) from the left
251 *
252  CALL dlarfb( 'Left', '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 dgeqr2p( 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 DGEQRFP
272 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
DLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: dlarft.f:165
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: dlarfb.f:197
subroutine dgeqr2p(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2P computes the QR factorization of a general rectangular matrix with non-negative diagonal elem...
Definition: dgeqr2p.f:126
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 dgeqrt ( integer  M,
integer  N,
integer  NB,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldt, * )  T,
integer  LDT,
double precision, dimension( * )  WORK,
integer  INFO 
)

DGEQRT

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

Purpose:
 DGEQRT computes a blocked QR factorization of a real 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dgeqrt.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  DOUBLE PRECISION 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 dgeqrt2, dgeqrt3, dlarfb, 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( 'DGEQRT', -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 dgeqrt3( m-i+1, ib, a(i,i), lda, t(1,i), ldt, iinfo )
202  ELSE
203  CALL dgeqrt2( 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**T to A(I:M,I+IB:N) from the left
208 *
209  CALL dlarfb( 'L', 'T', '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 DGEQRT
217 *
recursive subroutine dgeqrt3(M, N, A, LDA, T, LDT, INFO)
DGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact...
Definition: dgeqrt3.f:134
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: dlarfb.f:197
subroutine dgeqrt2(M, N, A, LDA, T, LDT, INFO)
DGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY represen...
Definition: dgeqrt2.f:129

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgeqrt2 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldt, * )  T,
integer  LDT,
integer  INFO 
)

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

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

Purpose:
 DGEQRT2 computes a QR factorization of a real 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 DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the real 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 DOUBLE PRECISION 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**T

  where V**T is the transpose of V.

Definition at line 129 of file dgeqrt2.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  DOUBLE PRECISION a( lda, * ), t( ldt, * )
140 * ..
141 *
142 * =====================================================================
143 *
144 * .. Parameters ..
145  DOUBLE PRECISION one, zero
146  parameter( one = 1.0d+00, zero = 0.0d+00 )
147 * ..
148 * .. Local Scalars ..
149  INTEGER i, k
150  DOUBLE PRECISION aii, alpha
151 * ..
152 * .. External Subroutines ..
153  EXTERNAL dlarfg, dgemv, dger, dtrmv, 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( 'DGEQRT2', -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 dlarfg( 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 dgemv( 'T',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 = -(t( i, 1 ))
197  CALL dger( 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)**T * A(I:M,I)
208 *
209  alpha = -t( i, 1 )
210  CALL dgemv( 'T', 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 dtrmv( '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 DGEQRT2
226 *
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
DTRMV
Definition: dtrmv.f:149

Here is the call graph for this function:

Here is the caller graph for this function:

recursive subroutine dgeqrt3 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldt, * )  T,
integer  LDT,
integer  INFO 
)

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

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

Purpose:
 DGEQRT3 recursively computes a QR factorization of a real 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 DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the real 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 DOUBLE PRECISION 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**T

  where V**T is the transpose of V.

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

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

DGERFS

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

Purpose:
 DGERFS 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 = 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDAF,N)
          The factors L and U from the factorization A = P*L*U
          as computed by DGETRF.
[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 DGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[in]B
          B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by DGETRS.
          On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]FERR
          FERR is DOUBLE PRECISION array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (3*N)
[out]IWORK
          IWORK is INTEGER 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 187 of file dgerfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

DGERFSX

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

Purpose:
    DGERFSX 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDAF,N)
     The factors L and U from the factorization A = P*L*U
     as computed by DGETRF.
[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 DGETRF; for 1<=i<=N, row i of the
     matrix was interchanged with row IPIV(i).
[in]R
          R is DOUBLE PRECISION array, dimension (N)
     The row scale factors for A.  If EQUED = 'R' or 'B', A is
     multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
     is not accessed.  
     If R is accessed, each element of R should be a power of the radix
     to ensure a reliable solution and error estimates. Scaling by
     powers of the radix does not cause rounding errors unless the
     result underflows or overflows. Rounding errors during scaling
     lead to refining with a matrix that is not equivalent to the
     input matrix, producing error estimates that may not be
     reliable.
[in]C
          C is DOUBLE PRECISION array, dimension (N)
     The column scale factors for A.  If EQUED = 'C' or 'B', A is
     multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
     is not accessed. 
     If C is accessed, each element of C should be a power of the radix
     to ensure a reliable solution and error estimates. Scaling by
     powers of the radix does not cause rounding errors unless the
     result underflows or overflows. Rounding errors during scaling
     lead to refining with a matrix that is not equivalent to the
     input matrix, producing error estimates that may not be
     reliable.
[in]B
          B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
     On entry, the solution matrix X, as computed by DGETRS.
     On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
     The leading dimension of the array X.  LDX >= max(1,N).
[out]RCOND
          RCOND is DOUBLE PRECISION
     Reciprocal scaled condition number.  This is an estimate of the
     reciprocal Skeel condition number of the matrix A after
     equilibration (if done).  If this is less than the machine
     precision (in particular, if it is zero), the matrix is singular
     to working precision.  Note that the error may still be small even
     if this number is very small and the matrix appears ill-
     conditioned.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
     Componentwise relative backward error.  This is the
     componentwise relative backward error of each solution vector X(j)
     (i.e., the smallest relative change in any element of A or B that
     makes X(j) an exact solution).
[in]N_ERR_BNDS
          N_ERR_BNDS is INTEGER
     Number of error bounds to return for each right hand side
     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
     ERR_BNDS_COMP below.
[out]ERR_BNDS_NORM
          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     normwise relative error, which is defined as follows:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
            will attempt to find a solution with small componentwise
            relative error in the double-precision algorithm.  Positive
            is true, 0.0 is false.
         Default: 1.0 (attempt componentwise convergence)
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (4*N)
[out]IWORK
          IWORK is INTEGER array, dimension (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 dgerfsx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgerq2 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 DGERQ2 computes an RQ factorization of a real 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 DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION 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(2) . . . H(k), where k = min(m,n).

  Each H(i) has the form

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

  where tau is a real scalar, and v is a real vector with
  v(n-k+i+1:n) = 0 and v(n-k+i) = 1; 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 dgerq2.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  DOUBLE PRECISION a( lda, * ), tau( * ), work( * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  DOUBLE PRECISION one
142  parameter( one = 1.0d+0 )
143 * ..
144 * .. Local Scalars ..
145  INTEGER i, k
146  DOUBLE PRECISION aii
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL dlarf, dlarfg, 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( 'DGERQ2', -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 dlarfg( n-k+i, a( m-k+i, n-k+i ), a( m-k+i, 1 ), lda,
179  $ tau( i ) )
180 *
181 * Apply H(i) to A(1:m-k+i-1,1:n-k+i) from the right
182 *
183  aii = a( m-k+i, n-k+i )
184  a( m-k+i, n-k+i ) = one
185  CALL dlarf( 'Right', m-k+i-1, n-k+i, a( m-k+i, 1 ), lda,
186  $ tau( i ), a, lda, work )
187  a( m-k+i, n-k+i ) = aii
188  10 CONTINUE
189  RETURN
190 *
191 * End of DGERQ2
192 *
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:126

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgerqf ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGERQF

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

Purpose:
 DGERQF computes an RQ factorization of a real 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 DOUBLE PRECISION 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
          orthogonal 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is DOUBLE PRECISION 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(2) . . . H(k), where k = min(m,n).

  Each H(i) has the form

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

  where tau is a real scalar, and v is a real vector with
  v(n-k+i+1:n) = 0 and v(n-k+i) = 1; 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 dgerqf.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  DOUBLE PRECISION 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 dgerq2, dlarfb, dlarft, 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, 'DGERQF', ' ', 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( 'DGERQF', -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, 'DGERQF', ' ', 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, 'DGERQF', ' ', 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 dgerq2( 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 dlarft( '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 dlarfb( '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 dgerq2( mu, nu, a, lda, tau, work, iinfo )
281 *
282  work( 1 ) = iws
283  RETURN
284 *
285 * End of DGERQF
286 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
DLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: dlarft.f:165
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: dlarfb.f:197
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine dgerq2(M, N, A, LDA, TAU, WORK, INFO)
DGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: dgerq2.f:125

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgesvj ( character*1  JOBA,
character*1  JOBU,
character*1  JOBV,
integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( n )  SVA,
integer  MV,
double precision, dimension( ldv, * )  V,
integer  LDV,
double precision, dimension( lwork )  WORK,
integer  LWORK,
integer  INFO 
)

DGESVJ

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

Purpose:
 DGESVJ computes the singular value decomposition (SVD) of a real
 M-by-N matrix A, where M >= N. The SVD of A is written as
                                    [++]   [xx]   [x0]   [xx]
              A = U * SIGMA * V^t,  [++] = [xx] * [ox] * [xx]
                                    [++]   [xx]
 where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
 matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
 of SIGMA are the singular values of A. The columns of U and V are the
 left and the right singular vectors of A, respectively.
 DGESVJ can sometimes compute tiny singular values and their singular vectors much
 more accurately than other SVD routines, see below under Further Details.
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=DSQRT(M), EPS=DLAMCH('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/DLAMCH('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 DOUBLE PRECISION 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 DLAMCH('S'). The singular
                 vectors corresponding to underflowed or zero singular
                 values are not computed. The value of RANKA is returned
                 in the array WORK as RANKA=NINT(WORK(2)). Also see the
                 descriptions of SVA and WORK. The computed columns of U
                 are mutually numerically orthogonal up to approximately
                 TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
                 see the description of JOBU.
                 If INFO .GT. 0 :
                 the procedure DGESVJ 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^T||_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 DGESVJ 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 DOUBLE PRECISION array, dimension (N)
          On exit :
          If INFO .EQ. 0 :
          depending on the value SCALE = WORK(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 DGESVJ 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 DGESVJ
          is applied to the first MV rows of V. See the description of JOBV.
[in,out]V
          V is DOUBLE PRECISION 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]WORK
          WORK is DOUBLE PRECISION array, dimension max(4,M+N).
          On entry :
          If JOBU .EQ. 'C' :
          WORK(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=DLAMCH('E').
                    It is required that CTOL >= ONE, i.e. it is not
                    allowed to force the routine to obtain orthogonality
                    below EPS.
          On exit :
          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
                    are the computed singular values of A.
                    (See description of SVA().)
          WORK(2) = NINT(WORK(2)) is the number of the computed nonzero
                    singular values.
          WORK(3) = NINT(WORK(3)) is the number of the computed singular
                    values that are larger than the underflow threshold.
          WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi
                    rotations needed for numerical convergence.
          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
                    This is useful information in cases when DGESVJ did
                    not converge, as it can be used to estimate whether
                    the output is stil useful and for post festum analysis.
          WORK(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]LWORK
          LWORK is INTEGER
          length of WORK, WORK >= MAX(6,M+N)
[out]INFO
          INFO is INTEGER
          = 0 : successful exit.
          < 0 : if INFO = -i, then the i-th argument had an illegal value
          > 0 : DGESVJ did not converge in the maximal allowed number (30)
                of sweeps. The output may still be useful. See the
                description of WORK.
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. The rotations are implemented as fast scaled rotations of
  Anda and Park [1]. In the case of underflow of the Jacobi angle, a
  modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
  column interchanges of de Rijk [2]. 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 [3].
  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 [5,6], and it is the kernel routine in the SIGMA library [7].
  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] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
     SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
 [2] 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.
 [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
 [4] 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.
 [5] 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.
 [6] 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.
 [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
     QSVD, (H,K)-SVD computations.
     Department of Mathematics, University of Zagreb, 2008.
Bugs, examples and comments:
  ===========================
  Please report all bugs and send interesting test examples and comments to
  drmac@math.hr. Thank you.

Definition at line 339 of file dgesvj.f.

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