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

Functions

subroutine sgeqpf (M, N, A, LDA, JPVT, TAU, WORK, INFO)
 SGEQPF More...
 
subroutine sgebak (JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
 SGEBAK More...
 
subroutine sgebal (JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
 SGEBAL More...
 
subroutine sgebd2 (M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO)
 SGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm. More...
 
subroutine sgebrd (M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
 SGEBRD More...
 
subroutine sgecon (NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
 SGECON More...
 
subroutine sgeequ (M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
 SGEEQU More...
 
subroutine sgeequb (M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
 SGEEQUB More...
 
subroutine sgehd2 (N, ILO, IHI, A, LDA, TAU, WORK, INFO)
 SGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm. More...
 
subroutine sgehrd (N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
 SGEHRD More...
 
subroutine sgelq2 (M, N, A, LDA, TAU, WORK, INFO)
 SGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine sgelqf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 SGELQF More...
 
subroutine sgemqrt (SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
 SGEMQRT More...
 
subroutine sgeql2 (M, N, A, LDA, TAU, WORK, INFO)
 SGEQL2 computes the QL factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine sgeqlf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 SGEQLF More...
 
subroutine sgeqp3 (M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
 SGEQP3 More...
 
subroutine sgeqr2 (M, N, A, LDA, TAU, WORK, INFO)
 SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine sgeqr2p (M, N, A, LDA, TAU, WORK, INFO)
 SGEQR2P computes the QR factorization of a general rectangular matrix with non-negative diagonal elements using an unblocked algorithm. More...
 
subroutine sgeqrf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 SGEQRF More...
 
subroutine sgeqrfp (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 SGEQRFP More...
 
subroutine sgeqrt (M, N, NB, A, LDA, T, LDT, WORK, INFO)
 SGEQRT More...
 
subroutine sgeqrt2 (M, N, A, LDA, T, LDT, INFO)
 SGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY representation of Q. More...
 
recursive subroutine sgeqrt3 (M, N, A, LDA, T, LDT, INFO)
 SGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact WY representation of Q. More...
 
subroutine sgerfs (TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
 SGERFS More...
 
subroutine sgerfsx (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)
 SGERFSX More...
 
subroutine sgerq2 (M, N, A, LDA, TAU, WORK, INFO)
 SGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm. More...
 
subroutine sgerqf (M, N, A, LDA, TAU, WORK, LWORK, INFO)
 SGERQF More...
 
subroutine sgesvj (JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
 SGESVJ More...
 
subroutine sgetf2 (M, N, A, LDA, IPIV, INFO)
 SGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row interchanges (unblocked algorithm). More...
 
subroutine sgetrf (M, N, A, LDA, IPIV, INFO)
 SGETRF More...
 
recursive subroutine sgetrf2 (M, N, A, LDA, IPIV, INFO)
 SGETRF2 More...
 
subroutine sgetri (N, A, LDA, IPIV, WORK, LWORK, INFO)
 SGETRI More...
 
subroutine sgetrs (TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 SGETRS More...
 
subroutine shgeqz (JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
 SHGEQZ More...
 
subroutine sla_geamv (TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
 SLA_GEAMV computes a matrix-vector product using a general matrix to calculate error bounds. More...
 
real function sla_gercond (TRANS, N, A, LDA, AF, LDAF, IPIV, CMODE, C, INFO, WORK, IWORK)
 SLA_GERCOND estimates the Skeel condition number for a general matrix. More...
 
subroutine sla_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)
 SLA_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...
 
real function sla_gerpvgrw (N, NCOLS, A, LDA, AF, LDAF)
 SLA_GERPVGRW More...
 
subroutine stgevc (SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
 STGEVC More...
 
subroutine stgexc (WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
 STGEXC More...
 

Detailed Description

This is the group of real computational functions for GE matrices

Function Documentation

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

SGEBAK

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

Purpose:
 SGEBAK forms the right or left eigenvectors of a real general matrix
 by backward transformation on the computed eigenvectors of the
 balanced matrix output by SGEBAL.
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 SGEBAL.
[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 SGEBAL.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
[in]SCALE
          SCALE is REAL array, dimension (N)
          Details of the permutation and scaling factors, as returned
          by SGEBAL.
[in]M
          M is INTEGER
          The number of columns of the matrix V.  M >= 0.
[in,out]V
          V is REAL array, dimension (LDV,M)
          On entry, the matrix of right or left eigenvectors to be
          transformed, as returned by SHSEIN or STREVC.
          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 sgebak.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  REAL v( ldv, * ), scale( * )
144 * ..
145 *
146 * =====================================================================
147 *
148 * .. Parameters ..
149  REAL one
150  parameter( one = 1.0e+0 )
151 * ..
152 * .. Local Scalars ..
153  LOGICAL leftv, rightv
154  INTEGER i, ii, k
155  REAL s
156 * ..
157 * .. External Functions ..
158  LOGICAL lsame
159  EXTERNAL lsame
160 * ..
161 * .. External Subroutines ..
162  EXTERNAL sscal, sswap, 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( 'SGEBAK', -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 sscal( 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 sscal( 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 sswap( 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 sswap( 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 SGEBAK
267 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGEBAL

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

Purpose:
 SGEBAL 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 REAL array, dimension (LDA,N)
          On entry, the input matrix A.
          On exit,  A is overwritten by the balanced matrix.
          If JOB = 'N', A is not referenced.
          See Further Details.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]ILO
          ILO is INTEGER
[out]IHI
          IHI is INTEGER
          ILO and IHI are set to integers such that on exit
          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
[out]SCALE
          SCALE is REAL array, dimension (N)
          Details of the permutations and scaling factors applied to
          A.  If P(j) is the index of the row and column interchanged
          with row and column j and D(j) is the scaling factor
          applied to row and column j, then
          SCALE(j) = P(j)    for j = 1,...,ILO-1
                   = D(j)    for j = ILO,...,IHI
                   = P(j)    for j = IHI+1,...,N.
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Further Details:
  The permutations consist of row and column interchanges which put
  the matrix in the form

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

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

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

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

  This subroutine is based on the EISPACK routine BALANC.

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

Definition at line 162 of file sgebal.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  REAL a( lda, * ), scale( * )
174 * ..
175 *
176 * =====================================================================
177 *
178 * .. Parameters ..
179  REAL zero, one
180  parameter( zero = 0.0e+0, one = 1.0e+0 )
181  REAL sclfac
182  parameter( sclfac = 2.0e+0 )
183  REAL factor
184  parameter( factor = 0.95e+0 )
185 * ..
186 * .. Local Scalars ..
187  LOGICAL noconv
188  INTEGER i, ica, iexc, ira, j, k, l, m
189  REAL c, ca, f, g, r, ra, s, sfmax1, sfmax2, sfmin1,
190  $ sfmin2
191 * ..
192 * .. External Functions ..
193  LOGICAL sisnan, lsame
194  INTEGER isamax
195  REAL slamch, snrm2
196  EXTERNAL sisnan, lsame, isamax, slamch, snrm2
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL sscal, sswap, 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( 'SGEBAL', -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 sswap( l, a( 1, j ), 1, a( 1, m ), 1 )
248  CALL sswap( 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 = slamch( 'S' ) / slamch( 'P' )
310  sfmax1 = one / sfmin1
311  sfmin2 = sfmin1*sclfac
312  sfmax2 = one / sfmin2
313  140 CONTINUE
314  noconv = .false.
315 *
316  DO 200 i = k, l
317 *
318  c = snrm2( l-k+1, a( k, i ), 1 )
319  r = snrm2( l-k+1, a( i, k ), lda )
320  ica = isamax( l, a( 1, i ), 1 )
321  ca = abs( a( ica, i ) )
322  ira = isamax( n-k+1, a( i, k ), lda )
323  ra = abs( a( i, ira+k-1 ) )
324 *
325 * Guard against zero C or R due to underflow.
326 *
327  IF( c.EQ.zero .OR. r.EQ.zero )
328  $ GO TO 200
329  g = r / sclfac
330  f = one
331  s = c + r
332  160 CONTINUE
333  IF( c.GE.g .OR. max( f, c, ca ).GE.sfmax2 .OR.
334  $ min( r, g, ra ).LE.sfmin2 )GO TO 170
335  f = f*sclfac
336  c = c*sclfac
337  ca = ca*sclfac
338  r = r / sclfac
339  g = g / sclfac
340  ra = ra / sclfac
341  GO TO 160
342 *
343  170 CONTINUE
344  g = c / sclfac
345  180 CONTINUE
346  IF( g.LT.r .OR. max( r, ra ).GE.sfmax2 .OR.
347  $ min( f, c, g, ca ).LE.sfmin2 )GO TO 190
348  IF( sisnan( c+f+ca+r+g+ra ) ) THEN
349 *
350 * Exit if NaN to avoid infinite loop
351 *
352  info = -3
353  CALL xerbla( 'SGEBAL', -info )
354  RETURN
355  END IF
356  f = f / sclfac
357  c = c / sclfac
358  g = g / sclfac
359  ca = ca / sclfac
360  r = r*sclfac
361  ra = ra*sclfac
362  GO TO 180
363 *
364 * Now balance.
365 *
366  190 CONTINUE
367  IF( ( c+r ).GE.factor*s )
368  $ GO TO 200
369  IF( f.LT.one .AND. scale( i ).LT.one ) THEN
370  IF( f*scale( i ).LE.sfmin1 )
371  $ GO TO 200
372  END IF
373  IF( f.GT.one .AND. scale( i ).GT.one ) THEN
374  IF( scale( i ).GE.sfmax1 / f )
375  $ GO TO 200
376  END IF
377  g = one / f
378  scale( i ) = scale( i )*f
379  noconv = .true.
380 *
381  CALL sscal( n-k+1, g, a( i, k ), lda )
382  CALL sscal( l, f, a( 1, i ), 1 )
383 *
384  200 CONTINUE
385 *
386  IF( noconv )
387  $ GO TO 140
388 *
389  210 CONTINUE
390  ilo = k
391  ihi = l
392 *
393  RETURN
394 *
395 * End of SGEBAL
396 *
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:56
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
logical function sisnan(SIN)
SISNAN tests input for NaN.
Definition: sisnan.f:61
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 SGEBD2 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 REAL 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 REAL array, dimension (min(M,N))
          The diagonal elements of the bidiagonal matrix B:
          D(i) = A(i,i).
[out]E
          E is REAL array, dimension (min(M,N)-1)
          The off-diagonal elements of the bidiagonal matrix B:
          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
[out]TAUQ
          TAUQ is REAL 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 REAL 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 REAL 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 sgebd2.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  REAL a( lda, * ), d( * ), e( * ), taup( * ),
202  $ tauq( * ), work( * )
203 * ..
204 *
205 * =====================================================================
206 *
207 * .. Parameters ..
208  REAL zero, one
209  parameter( zero = 0.0e+0, one = 1.0e+0 )
210 * ..
211 * .. Local Scalars ..
212  INTEGER i
213 * ..
214 * .. External Subroutines ..
215  EXTERNAL slarf, slarfg, 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( 'SGEBD2', -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 slarfg( 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 slarf( '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 slarfg( 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 slarf( '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 slarfg( 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 slarf( '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 slarfg( 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 slarf( '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 SGEBD2
319 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:126
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGEBRD

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

Purpose:
 SGEBRD 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 REAL 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 REAL array, dimension (min(M,N))
          The diagonal elements of the bidiagonal matrix B:
          D(i) = A(i,i).
[out]E
          E is REAL array, dimension (min(M,N)-1)
          The off-diagonal elements of the bidiagonal matrix B:
          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
[out]TAUQ
          TAUQ is REAL 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 REAL 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 REAL 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 sgebrd.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  REAL a( lda, * ), d( * ), e( * ), taup( * ),
218  $ tauq( * ), work( * )
219 * ..
220 *
221 * =====================================================================
222 *
223 * .. Parameters ..
224  REAL one
225  parameter( one = 1.0e+0 )
226 * ..
227 * .. Local Scalars ..
228  LOGICAL lquery
229  INTEGER i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb,
230  $ nbmin, nx
231  REAL ws
232 * ..
233 * .. External Subroutines ..
234  EXTERNAL sgebd2, sgemm, slabrd, xerbla
235 * ..
236 * .. Intrinsic Functions ..
237  INTRINSIC max, min, real
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, 'SGEBRD', ' ', m, n, -1, -1 ) )
249  lwkopt = ( m+n )*nb
250  work( 1 ) = REAL( 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( 'SGEBRD', -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, 'SGEBRD', ' ', 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, 'SGEBRD', ' ', 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 slabrd( 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 sgemm( '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 sgemm( '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 sgebd2( 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 SGEBRD
352 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgebd2(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO)
SGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm.
Definition: sgebd2.f:191
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine slabrd(M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY)
SLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form.
Definition: slabrd.f:212

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sgecon ( character  NORM,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real  ANORM,
real  RCOND,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

SGECON

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

Purpose:
 SGECON 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 SGETRF.

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

126 *
127 * -- LAPACK computational routine (version 3.4.0) --
128 * -- LAPACK is a software package provided by Univ. of Tennessee, --
129 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
130 * November 2011
131 *
132 * .. Scalar Arguments ..
133  CHARACTER norm
134  INTEGER info, lda, n
135  REAL anorm, rcond
136 * ..
137 * .. Array Arguments ..
138  INTEGER iwork( * )
139  REAL a( lda, * ), work( * )
140 * ..
141 *
142 * =====================================================================
143 *
144 * .. Parameters ..
145  REAL one, zero
146  parameter( one = 1.0e+0, zero = 0.0e+0 )
147 * ..
148 * .. Local Scalars ..
149  LOGICAL onenrm
150  CHARACTER normin
151  INTEGER ix, kase, kase1
152  REAL ainvnm, scale, sl, smlnum, su
153 * ..
154 * .. Local Arrays ..
155  INTEGER isave( 3 )
156 * ..
157 * .. External Functions ..
158  LOGICAL lsame
159  INTEGER isamax
160  REAL slamch
161  EXTERNAL lsame, isamax, slamch
162 * ..
163 * .. External Subroutines ..
164  EXTERNAL slacn2, slatrs, srscl, 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( 'SGECON', -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 = slamch( '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 slacn2( 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 slatrs( '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 slatrs( '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 slatrs( '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 slatrs( '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 = isamax( n, work, 1 )
244  IF( scale.LT.abs( work( ix ) )*smlnum .OR. scale.EQ.zero )
245  $ GO TO 20
246  CALL srscl( 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 SGECON
260 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
SLATRS solves a triangular system of equations with the scale factor set to prevent overflow...
Definition: slatrs.f:240
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine srscl(N, SA, SX, INCX)
SRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: srscl.f:86
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: slacn2.f:138

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGEEQU

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

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

Definition at line 141 of file sgeequ.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  REAL amax, colcnd, rowcnd
150 * ..
151 * .. Array Arguments ..
152  REAL a( lda, * ), c( * ), r( * )
153 * ..
154 *
155 * =====================================================================
156 *
157 * .. Parameters ..
158  REAL one, zero
159  parameter( one = 1.0e+0, zero = 0.0e+0 )
160 * ..
161 * .. Local Scalars ..
162  INTEGER i, j
163  REAL bignum, rcmax, rcmin, smlnum
164 * ..
165 * .. External Functions ..
166  REAL slamch
167  EXTERNAL slamch
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( 'SGEEQU', -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 = slamch( '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 SGEEQU
303 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGEEQUB

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

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

Definition at line 148 of file sgeequb.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  REAL amax, colcnd, rowcnd
157 * ..
158 * .. Array Arguments ..
159  REAL a( lda, * ), c( * ), r( * )
160 * ..
161 *
162 * =====================================================================
163 *
164 * .. Parameters ..
165  REAL one, zero
166  parameter( one = 1.0e+0, zero = 0.0e+0 )
167 * ..
168 * .. Local Scalars ..
169  INTEGER i, j
170  REAL bignum, rcmax, rcmin, smlnum, radix, logrdx
171 * ..
172 * .. External Functions ..
173  REAL slamch
174  EXTERNAL slamch
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( 'SGEEQUB', -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 = slamch( 'S' )
211  bignum = one / smlnum
212  radix = slamch( '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 SGEEQUB
320 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 SGEHD2 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 SGEBAL; otherwise they should be
          set to 1 and N respectively. See Further Details.
          1 <= ILO <= IHI <= max(1,N).
[in,out]A
          A is REAL 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 REAL array, dimension (N-1)
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
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 sgehd2.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  REAL a( lda, * ), tau( * ), work( * )
162 * ..
163 *
164 * =====================================================================
165 *
166 * .. Parameters ..
167  REAL one
168  parameter( one = 1.0e+0 )
169 * ..
170 * .. Local Scalars ..
171  INTEGER i
172  REAL aii
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL slarf, slarfg, 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( 'SGEHD2', -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 slarfg( 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 slarf( '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 slarf( '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 SGEHD2
224 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:126
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGEHRD

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

Purpose:
 SGEHRD 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 SGEBAL; 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 REAL 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 REAL 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 REAL 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 sgehrd.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  REAL 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  REAL zero, one
189  parameter( zero = 0.0e+0,
190  $ one = 1.0e+0 )
191 * ..
192 * .. Local Scalars ..
193  LOGICAL lquery
194  INTEGER i, ib, iinfo, iwt, j, ldwork, lwkopt, nb,
195  $ nbmin, nh, nx
196  REAL ei
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL saxpy, sgehd2, sgemm, slahr2, slarfb, strmm,
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, 'SGEHRD', ' ', 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( 'SGEHRD', -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, 'SGEHRD', ' ', 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, 'SGEHRD', ' ', 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, 'SGEHRD', ' ', 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 slahr2( 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 sgemm( '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 strmm( '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 saxpy( 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 slarfb( '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 sgehd2( n, i, ihi, a, lda, tau, work, iinfo )
350  work( 1 ) = lwkopt
351 *
352  RETURN
353 *
354 * End of SGEHRD
355 *
subroutine slahr2(N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)
SLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elemen...
Definition: slahr2.f:183
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:179
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
subroutine sgehd2(N, ILO, IHI, A, LDA, TAU, WORK, INFO)
SGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm...
Definition: sgehd2.f:151
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine slarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: slarfb.f:197

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sgelq2 ( integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 SGELQ2 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL 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 sgelq2.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  REAL a( lda, * ), tau( * ), work( * )
134 * ..
135 *
136 * =====================================================================
137 *
138 * .. Parameters ..
139  REAL one
140  parameter( one = 1.0e+0 )
141 * ..
142 * .. Local Scalars ..
143  INTEGER i, k
144  REAL aii
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL slarf, slarfg, 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( 'SGELQ2', -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 slarfg( 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 slarf( '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 SGELQ2
191 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:126
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGELQF

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

Purpose:
 SGELQF 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL 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 sgelqf.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  REAL 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 sgelq2, slarfb, slarft, 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, 'SGELQF', ' ', 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( 'SGELQF', -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, 'SGELQF', ' ', 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, 'SGELQF', ' ', 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 sgelq2( 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 slarft( '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 slarfb( '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 sgelq2( 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 SGELQF
268 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgelq2(M, N, A, LDA, TAU, WORK, INFO)
SGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: sgelq2.f:123
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine slarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: slarfb.f:197
subroutine slarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: slarft.f:165

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGEMQRT

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

Purpose:
 SGEMQRT 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 SGEQRT. 

 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;
          = 'T':  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 REAL 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 REAL 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 REAL 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 REAL 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 sgemqrt.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  REAL 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, slarfb
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( 'SGEMQRT', -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 slarfb( '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 slarfb( '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 slarfb( '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 slarfb( '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 SGEMQRT
290 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine slarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: slarfb.f:197

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sgeql2 ( integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 SGEQL2 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
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 sgeql2.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  REAL a( lda, * ), tau( * ), work( * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  REAL one
142  parameter( one = 1.0e+0 )
143 * ..
144 * .. Local Scalars ..
145  INTEGER i, k
146  REAL aii
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL slarf, slarfg, 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( 'SGEQL2', -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 slarfg( 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 slarf( '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 SGEQL2
192 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:126
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGEQLF

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

Purpose:
 SGEQLF 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL 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 sgeqlf.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  REAL 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 sgeql2, slarfb, slarft, 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, 'SGEQLF', ' ', 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( 'SGEQLF', -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, 'SGEQLF', ' ', 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, 'SGEQLF', ' ', 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 sgeql2( 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 slarft( '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 slarfb( '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 sgeql2( mu, nu, a, lda, tau, work, iinfo )
281 *
282  work( 1 ) = iws
283  RETURN
284 *
285 * End of SGEQLF
286 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgeql2(M, N, A, LDA, TAU, WORK, INFO)
SGEQL2 computes the QL factorization of a general rectangular matrix using an unblocked algorithm...
Definition: sgeql2.f:125
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine slarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: slarfb.f:197
subroutine slarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: slarft.f:165

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sgeqp3 ( integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

SGEQP3

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

Purpose:
 SGEQP3 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[out]WORK
          WORK is REAL 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 sgeqp3.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  REAL 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 sgeqrf, slaqp2, slaqps, sormqr, sswap, xerbla
180 * ..
181 * .. External Functions ..
182  INTEGER ilaenv
183  REAL snrm2
184  EXTERNAL ilaenv, snrm2
185 * ..
186 * .. Intrinsic Functions ..
187  INTRINSIC int, max, min
188 * Test input arguments
189 * ====================
190 *
191  info = 0
192  lquery = ( lwork.EQ.-1 )
193  IF( m.LT.0 ) THEN
194  info = -1
195  ELSE IF( n.LT.0 ) THEN
196  info = -2
197  ELSE IF( lda.LT.max( 1, m ) ) THEN
198  info = -4
199  END IF
200 *
201  IF( info.EQ.0 ) THEN
202  minmn = min( m, n )
203  IF( minmn.EQ.0 ) THEN
204  iws = 1
205  lwkopt = 1
206  ELSE
207  iws = 3*n + 1
208  nb = ilaenv( inb, 'SGEQRF', ' ', m, n, -1, -1 )
209  lwkopt = 2*n + ( n + 1 )*nb
210  END IF
211  work( 1 ) = lwkopt
212 *
213  IF( ( lwork.LT.iws ) .AND. .NOT.lquery ) THEN
214  info = -8
215  END IF
216  END IF
217 *
218  IF( info.NE.0 ) THEN
219  CALL xerbla( 'SGEQP3', -info )
220  RETURN
221  ELSE IF( lquery ) THEN
222  RETURN
223  END IF
224 *
225 * Move initial columns up front.
226 *
227  nfxd = 1
228  DO 10 j = 1, n
229  IF( jpvt( j ).NE.0 ) THEN
230  IF( j.NE.nfxd ) THEN
231  CALL sswap( m, a( 1, j ), 1, a( 1, nfxd ), 1 )
232  jpvt( j ) = jpvt( nfxd )
233  jpvt( nfxd ) = j
234  ELSE
235  jpvt( j ) = j
236  END IF
237  nfxd = nfxd + 1
238  ELSE
239  jpvt( j ) = j
240  END IF
241  10 CONTINUE
242  nfxd = nfxd - 1
243 *
244 * Factorize fixed columns
245 * =======================
246 *
247 * Compute the QR factorization of fixed columns and update
248 * remaining columns.
249 *
250  IF( nfxd.GT.0 ) THEN
251  na = min( m, nfxd )
252 *CC CALL SGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
253  CALL sgeqrf( m, na, a, lda, tau, work, lwork, info )
254  iws = max( iws, int( work( 1 ) ) )
255  IF( na.LT.n ) THEN
256 *CC CALL SORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
257 *CC $ TAU, A( 1, NA+1 ), LDA, WORK, INFO )
258  CALL sormqr( 'Left', 'Transpose', m, n-na, na, a, lda, tau,
259  $ a( 1, na+1 ), lda, work, lwork, info )
260  iws = max( iws, int( work( 1 ) ) )
261  END IF
262  END IF
263 *
264 * Factorize free columns
265 * ======================
266 *
267  IF( nfxd.LT.minmn ) THEN
268 *
269  sm = m - nfxd
270  sn = n - nfxd
271  sminmn = minmn - nfxd
272 *
273 * Determine the block size.
274 *
275  nb = ilaenv( inb, 'SGEQRF', ' ', sm, sn, -1, -1 )
276  nbmin = 2
277  nx = 0
278 *
279  IF( ( nb.GT.1 ) .AND. ( nb.LT.sminmn ) ) THEN
280 *
281 * Determine when to cross over from blocked to unblocked code.
282 *
283  nx = max( 0, ilaenv( ixover, 'SGEQRF', ' ', sm, sn, -1,
284  $ -1 ) )
285 *
286 *
287  IF( nx.LT.sminmn ) THEN
288 *
289 * Determine if workspace is large enough for blocked code.
290 *
291  minws = 2*sn + ( sn+1 )*nb
292  iws = max( iws, minws )
293  IF( lwork.LT.minws ) THEN
294 *
295 * Not enough workspace to use optimal NB: Reduce NB and
296 * determine the minimum value of NB.
297 *
298  nb = ( lwork-2*sn ) / ( sn+1 )
299  nbmin = max( 2, ilaenv( inbmin, 'SGEQRF', ' ', sm, sn,
300  $ -1, -1 ) )
301 *
302 *
303  END IF
304  END IF
305  END IF
306 *
307 * Initialize partial column norms. The first N elements of work
308 * store the exact column norms.
309 *
310  DO 20 j = nfxd + 1, n
311  work( j ) = snrm2( sm, a( nfxd+1, j ), 1 )
312  work( n+j ) = work( j )
313  20 CONTINUE
314 *
315  IF( ( nb.GE.nbmin ) .AND. ( nb.LT.sminmn ) .AND.
316  $ ( nx.LT.sminmn ) ) THEN
317 *
318 * Use blocked code initially.
319 *
320  j = nfxd + 1
321 *
322 * Compute factorization: while loop.
323 *
324 *
325  topbmn = minmn - nx
326  30 CONTINUE
327  IF( j.LE.topbmn ) THEN
328  jb = min( nb, topbmn-j+1 )
329 *
330 * Factorize JB columns among columns J:N.
331 *
332  CALL slaqps( m, n-j+1, j-1, jb, fjb, a( 1, j ), lda,
333  $ jpvt( j ), tau( j ), work( j ), work( n+j ),
334  $ work( 2*n+1 ), work( 2*n+jb+1 ), n-j+1 )
335 *
336  j = j + fjb
337  GO TO 30
338  END IF
339  ELSE
340  j = nfxd + 1
341  END IF
342 *
343 * Use unblocked code to factor the last or only block.
344 *
345 *
346  IF( j.LE.minmn )
347  $ CALL slaqp2( m, n-j+1, j-1, a( 1, j ), lda, jpvt( j ),
348  $ tau( j ), work( j ), work( n+j ),
349  $ work( 2*n+1 ) )
350 *
351  END IF
352 *
353  work( 1 ) = iws
354  RETURN
355 *
356 * End of SGEQP3
357 *
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:56
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:170
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:138
subroutine slaqp2(M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)
SLAQP2 computes a QR factorization with column pivoting of the matrix block.
Definition: slaqp2.f:151
subroutine slaqps(M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)
SLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BL...
Definition: slaqps.f:180
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
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 sgeqpf ( integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  INFO 
)

SGEQPF

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

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

 SGEQPF 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[out]WORK
          WORK is REAL 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 sgeqpf.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  REAL a( lda, * ), tau( * ), work( * )
156 * ..
157 *
158 * =====================================================================
159 *
160 * .. Parameters ..
161  REAL zero, one
162  parameter( zero = 0.0e+0, one = 1.0e+0 )
163 * ..
164 * .. Local Scalars ..
165  INTEGER i, itemp, j, ma, mn, pvt
166  REAL aii, temp, temp2, tol3z
167 * ..
168 * .. External Subroutines ..
169  EXTERNAL sgeqr2, slarf, slarfg, sorm2r, sswap, xerbla
170 * ..
171 * .. Intrinsic Functions ..
172  INTRINSIC abs, max, min, sqrt
173 * ..
174 * .. External Functions ..
175  INTEGER isamax
176  REAL slamch, snrm2
177  EXTERNAL isamax, slamch, snrm2
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( 'SGEQPF', -info )
193  RETURN
194  END IF
195 *
196  mn = min( m, n )
197  tol3z = sqrt(slamch('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 sswap( 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 sgeqr2( m, ma, a, lda, tau, work, info )
223  IF( ma.LT.n ) THEN
224  CALL sorm2r( '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 ) = snrm2( 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 ) + isamax( n-i+1, work( i ), 1 )
246 *
247  IF( pvt.NE.i ) THEN
248  CALL sswap( 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 slarfg( m-i+1, a( i, i ), a( i+1, i ), 1, tau( i ) )
260  ELSE
261  CALL slarfg( 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 slarf( '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 ) = snrm2( 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 SGEQPF
305 *
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:56
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:126
subroutine sorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: sorm2r.f:161
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:108
subroutine sgeqr2(M, N, A, LDA, TAU, WORK, INFO)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: sgeqr2.f:123

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sgeqr2 ( integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 SGEQR2 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
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 sgeqr2.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  REAL a( lda, * ), tau( * ), work( * )
134 * ..
135 *
136 * =====================================================================
137 *
138 * .. Parameters ..
139  REAL one
140  parameter( one = 1.0e+0 )
141 * ..
142 * .. Local Scalars ..
143  INTEGER i, k
144  REAL aii
145 * ..
146 * .. External Subroutines ..
147  EXTERNAL slarf, slarfg, 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( 'SGEQR2', -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 slarfg( 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 slarf( '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 SGEQR2
191 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:126
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sgeqr2p ( integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 SGEQR2P 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
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 sgeqr2p.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  REAL a( lda, * ), tau( * ), work( * )
137 * ..
138 *
139 * =====================================================================
140 *
141 * .. Parameters ..
142  REAL one
143  parameter( one = 1.0e+0 )
144 * ..
145 * .. Local Scalars ..
146  INTEGER i, k
147  REAL aii
148 * ..
149 * .. External Subroutines ..
150  EXTERNAL slarf, slarfgp, 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( 'SGEQR2P', -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 slarfgp( 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 slarf( '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 SGEQR2P
194 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:126
subroutine slarfgp(N, ALPHA, X, INCX, TAU)
SLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: slarfgp.f:106

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGEQRF

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

Purpose:
 SGEQRF 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL 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 sgeqrf.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  REAL 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 sgeqr2, slarfb, slarft, 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, 'SGEQRF', ' ', 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( 'SGEQRF', -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, 'SGEQRF', ' ', 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, 'SGEQRF', ' ', 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 sgeqr2( 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 slarft( '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 slarfb( '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 sgeqr2( 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 SGEQRF
269 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgeqr2(M, N, A, LDA, TAU, WORK, INFO)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: sgeqr2.f:123
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine slarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: slarfb.f:197
subroutine slarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: slarft.f:165

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGEQRFP

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

Purpose:
 SGEQRFP 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL 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 sgeqrfp.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  REAL 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 sgeqr2p, slarfb, slarft, 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, 'SGEQRF', ' ', 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( 'SGEQRFP', -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, 'SGEQRF', ' ', 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, 'SGEQRF', ' ', 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 sgeqr2p( 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 slarft( '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 slarfb( '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 sgeqr2p( 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 SGEQRFP
272 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine slarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: slarfb.f:197
subroutine sgeqr2p(M, N, A, LDA, TAU, WORK, INFO)
SGEQR2P computes the QR factorization of a general rectangular matrix with non-negative diagonal elem...
Definition: sgeqr2p.f:126
subroutine slarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: slarft.f:165

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sgeqrt ( integer  M,
integer  N,
integer  NB,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldt, * )  T,
integer  LDT,
real, dimension( * )  WORK,
integer  INFO 
)

SGEQRT

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

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

Purpose:
 SGEQRT2 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 REAL 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 REAL 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 sgeqrt2.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  REAL a( lda, * ), t( ldt, * )
140 * ..
141 *
142 * =====================================================================
143 *
144 * .. Parameters ..
145  REAL one, zero
146  parameter( one = 1.0, zero = 0.0 )
147 * ..
148 * .. Local Scalars ..
149  INTEGER i, k
150  REAL aii, alpha
151 * ..
152 * .. External Subroutines ..
153  EXTERNAL slarfg, sgemv, sger, strmv, 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( 'SGEQRT2', -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 slarfg( 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 sgemv( '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 sger( 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 sgemv( '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 strmv( '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 SGEQRT2
226 *
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:132
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:108
subroutine strmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
STRMV
Definition: strmv.f:149

Here is the call graph for this function:

Here is the caller graph for this function:

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

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

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

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGERFS

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

Purpose:
 SGERFS 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 REAL 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 REAL array, dimension (LDAF,N)
          The factors L and U from the factorization A = P*L*U
          as computed by SGETRF.
[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 SGETRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[in]B
          B is REAL 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 REAL array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by SGETRS.
          On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]FERR
          FERR is REAL array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.
[out]BERR
          BERR is REAL array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).
[out]WORK
          WORK is REAL 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 sgerfs.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  REAL 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  REAL zero
209  parameter( zero = 0.0e+0 )
210  REAL one
211  parameter( one = 1.0e+0 )
212  REAL two
213  parameter( two = 2.0e+0 )
214  REAL three
215  parameter( three = 3.0e+0 )
216 * ..
217 * .. Local Scalars ..
218  LOGICAL notran
219  CHARACTER transt
220  INTEGER count, i, j, k, kase, nz
221  REAL eps, lstres, s, safe1, safe2, safmin, xk
222 * ..
223 * .. Local Arrays ..
224  INTEGER isave( 3 )
225 * ..
226 * .. External Subroutines ..
227  EXTERNAL saxpy, scopy, sgemv, sgetrs, slacn2, xerbla
228 * ..
229 * .. Intrinsic Functions ..
230  INTRINSIC abs, max
231 * ..
232 * .. External Functions ..
233  LOGICAL lsame
234  REAL slamch
235  EXTERNAL lsame, slamch
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( 'SGERFS', -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 = slamch( 'Epsilon' )
284  safmin = slamch( '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 scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
302  CALL sgemv( 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 sgetrs( trans, n, 1, af, ldaf, ipiv, work( n+1 ), n,
359  $ info )
360  CALL saxpy( 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 SLACN2 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 slacn2( 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 sgetrs( 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 sgetrs( 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 SGERFS
437 *
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
SGETRS
Definition: sgetrs.f:123
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: slacn2.f:138

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGERFSX

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

416 *
417 * -- LAPACK computational routine (version 3.4.0) --
418 * -- LAPACK is a software package provided by Univ. of Tennessee, --
419 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
420 * November 2011
421 *
422 * .. Scalar Arguments ..
423  CHARACTER trans, equed
424  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs, nparams,
425  $ n_err_bnds
426  REAL rcond
427 * ..
428 * .. Array Arguments ..
429  INTEGER ipiv( * ), iwork( * )
430  REAL a( lda, * ), af( ldaf, * ), b( ldb, * ),
431  $ x( ldx , * ), work( * )
432  REAL r( * ), c( * ), params( * ), berr( * ),
433  $ err_bnds_norm( nrhs, * ),
434  $ err_bnds_comp( nrhs, * )
435 * ..
436 *
437 * ==================================================================
438 *
439 * .. Parameters ..
440  REAL zero, one
441  parameter( zero = 0.0e+0, one = 1.0e+0 )
442  REAL itref_default, ithresh_default,
443  $ componentwise_default
444  REAL rthresh_default, dzthresh_default
445  parameter( itref_default = 1.0 )
446  parameter( ithresh_default = 10.0 )
447  parameter( componentwise_default = 1.0 )
448  parameter( rthresh_default = 0.5 )
449  parameter( dzthresh_default = 0.25 )
450  INTEGER la_linrx_itref_i, la_linrx_ithresh_i,
451  $ la_linrx_cwise_i
452  parameter( la_linrx_itref_i = 1,
453  $ la_linrx_ithresh_i = 2 )
454  parameter( la_linrx_cwise_i = 3 )
455  INTEGER la_linrx_trust_i, la_linrx_err_i,
456  $ la_linrx_rcond_i
457  parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
458  parameter( la_linrx_rcond_i = 3 )
459 * ..
460 * .. Local Scalars ..
461  CHARACTER(1) norm
462  LOGICAL rowequ, colequ, notran
463  INTEGER j, trans_type, prec_type, ref_type
464  INTEGER n_norms
465  REAL anorm, rcond_tmp
466  REAL illrcond_thresh, err_lbnd, cwise_wrong
467  LOGICAL ignore_cwise
468  INTEGER ithresh
469  REAL rthresh, unstable_thresh
470 * ..
471 * .. External Subroutines ..
473 * ..
474 * .. Intrinsic Functions ..
475  INTRINSIC max, sqrt
476 * ..
477 * .. External Functions ..
478  EXTERNAL lsame, blas_fpinfo_x, ilatrans, ilaprec
479  EXTERNAL slamch, slange, sla_gercond
480  REAL slamch, slange, sla_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.0 ) THEN
494  params( la_linrx_itref_i ) = itref_default
495  ELSE
496  ref_type = params( la_linrx_itref_i )
497  END IF
498  END IF
499 *
500 * Set default parameters.
501 *
502  illrcond_thresh = REAL( N ) * slamch( 'Epsilon' )
503  ithresh = int( ithresh_default )
504  rthresh = rthresh_default
505  unstable_thresh = dzthresh_default
506  ignore_cwise = componentwise_default .EQ. 0.0
507 *
508  IF ( nparams.GE.la_linrx_ithresh_i ) THEN
509  IF ( params( la_linrx_ithresh_i ).LT.0.0 ) THEN
510  params( la_linrx_ithresh_i ) = ithresh
511  ELSE
512  ithresh = int( params( la_linrx_ithresh_i ) )
513  END IF
514  END IF
515  IF ( nparams.GE.la_linrx_cwise_i ) THEN
516  IF ( params( la_linrx_cwise_i ).LT.0.0 ) THEN
517  IF ( ignore_cwise ) THEN
518  params( la_linrx_cwise_i ) = 0.0
519  ELSE
520  params( la_linrx_cwise_i ) = 1.0
521  END IF
522  ELSE
523  ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0
524  END IF
525  END IF
526  IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
527  n_norms = 0
528  ELSE IF ( ignore_cwise ) THEN
529  n_norms = 1
530  ELSE
531  n_norms = 2
532  END IF
533 *
534  notran = lsame( trans, 'N' )
535  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
536  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
537 *
538 * Test input parameters.
539 *
540  IF( trans_type.EQ.-1 ) THEN
541  info = -1
542  ELSE IF( .NOT.rowequ .AND. .NOT.colequ .AND.
543  $ .NOT.lsame( equed, 'N' ) ) THEN
544  info = -2
545  ELSE IF( n.LT.0 ) THEN
546  info = -3
547  ELSE IF( nrhs.LT.0 ) THEN
548  info = -4
549  ELSE IF( lda.LT.max( 1, n ) ) THEN
550  info = -6
551  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
552  info = -8
553  ELSE IF( ldb.LT.max( 1, n ) ) THEN
554  info = -13
555  ELSE IF( ldx.LT.max( 1, n ) ) THEN
556  info = -15
557  END IF
558  IF( info.NE.0 ) THEN
559  CALL xerbla( 'SGERFSX', -info )
560  RETURN
561  END IF
562 *
563 * Quick return if possible.
564 *
565  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
566  rcond = 1.0
567  DO j = 1, nrhs
568  berr( j ) = 0.0
569  IF ( n_err_bnds .GE. 1 ) THEN
570  err_bnds_norm( j, la_linrx_trust_i) = 1.0
571  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
572  END IF
573  IF ( n_err_bnds .GE. 2 ) THEN
574  err_bnds_norm( j, la_linrx_err_i) = 0.0
575  err_bnds_comp( j, la_linrx_err_i ) = 0.0
576  END IF
577  IF ( n_err_bnds .GE. 3 ) THEN
578  err_bnds_norm( j, la_linrx_rcond_i) = 1.0
579  err_bnds_comp( j, la_linrx_rcond_i ) = 1.0
580  END IF
581  END DO
582  RETURN
583  END IF
584 *
585 * Default to failure.
586 *
587  rcond = 0.0
588  DO j = 1, nrhs
589  berr( j ) = 1.0
590  IF ( n_err_bnds .GE. 1 ) THEN
591  err_bnds_norm( j, la_linrx_trust_i ) = 1.0
592  err_bnds_comp( j, la_linrx_trust_i ) = 1.0
593  END IF
594  IF ( n_err_bnds .GE. 2 ) THEN
595  err_bnds_norm( j, la_linrx_err_i ) = 1.0
596  err_bnds_comp( j, la_linrx_err_i ) = 1.0
597  END IF
598  IF ( n_err_bnds .GE. 3 ) THEN
599  err_bnds_norm( j, la_linrx_rcond_i ) = 0.0
600  err_bnds_comp( j, la_linrx_rcond_i ) = 0.0
601  END IF
602  END DO
603 *
604 * Compute the norm of A and the reciprocal of the condition
605 * number of A.
606 *
607  IF( notran ) THEN
608  norm = 'I'
609  ELSE
610  norm = '1'
611  END IF
612  anorm = slange( norm, n, n, a, lda, work )
613  CALL sgecon( 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( 'D' )
620 
621  IF ( notran ) THEN
622  CALL sla_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 sla_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.0, sqrt( REAL( N ) ) ) * slamch( '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 = sla_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 = sla_gercond( trans, n, a, lda, af, ldaf, ipiv,
648  $ -1, r, info, work, iwork )
649  ELSE
650  rcond_tmp = sla_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.0 )
659  $ err_bnds_norm( j, la_linrx_err_i ) = 1.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.0
665  err_bnds_norm( j, la_linrx_trust_i ) = 0.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.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( slamch( 'Epsilon' ) )
692  DO j = 1, nrhs
693  IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
694  $ THEN
695  rcond_tmp = sla_gercond( trans, n, a, lda, af, ldaf,
696  $ ipiv, 1, x(1,j), info, work, iwork )
697  ELSE
698  rcond_tmp = 0.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.0 )
705  $ err_bnds_comp( j, la_linrx_err_i ) = 1.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.0
711  err_bnds_comp( j, la_linrx_trust_i ) = 0.0
712  IF ( params( la_linrx_cwise_i ) .EQ. 1.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.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 SGERFSX
731 *
integer function ilatrans(TRANS)
ILATRANS
Definition: ilatrans.f:60
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
real function sla_gercond(TRANS, N, A, LDA, AF, LDAF, IPIV, CMODE, C, INFO, WORK, IWORK)
SLA_GERCOND estimates the Skeel condition number for a general matrix.
Definition: sla_gercond.f:152
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
integer function ilaprec(PREC)
ILAPREC
Definition: ilaprec.f:60
subroutine sgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
SGECON
Definition: sgecon.f:126
subroutine sla_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)
SLA_GERFSX_EXTENDED improves the computed solution to a system of linear equations for general matric...
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:116

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sgerq2 ( integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose:
 SGERQ2 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL 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 sgerq2.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  REAL a( lda, * ), tau( * ), work( * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  REAL one
142  parameter( one = 1.0e+0 )
143 * ..
144 * .. Local Scalars ..
145  INTEGER i, k
146  REAL aii
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL slarf, slarfg, 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( 'SGERQ2', -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 slarfg( 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 slarf( '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 SGERQ2
192 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition: slarf.f:126
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGERQF

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

Purpose:
 SGERQF 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 REAL 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 REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is REAL 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 sgerqf.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  REAL 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 sgerq2, slarfb, slarft, 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  ELSE IF( lwork.LT.max( 1, m ) .AND. .NOT.lquery ) THEN
183  info = -7
184  END IF
185 *
186  IF( info.EQ.0 ) THEN
187  k = min( m, n )
188  IF( k.EQ.0 ) THEN
189  lwkopt = 1
190  ELSE
191  nb = ilaenv( 1, 'SGERQF', ' ', m, n, -1, -1 )
192  lwkopt = m*nb
193  work( 1 ) = lwkopt
194  END IF
195  work( 1 ) = lwkopt
196 *
197  IF( lwork.LT.max( 1, m ) .AND. .NOT.lquery ) THEN
198  info = -7
199  END IF
200  END IF
201 *
202  IF( info.NE.0 ) THEN
203  CALL xerbla( 'SGERQF', -info )
204  RETURN
205  ELSE IF( lquery ) THEN
206  RETURN
207  END IF
208 *
209 * Quick return if possible
210 *
211  IF( k.EQ.0 ) THEN
212  RETURN
213  END IF
214 *
215  nbmin = 2
216  nx = 1
217  iws = m
218  IF( nb.GT.1 .AND. nb.LT.k ) THEN
219 *
220 * Determine when to cross over from blocked to unblocked code.
221 *
222  nx = max( 0, ilaenv( 3, 'SGERQF', ' ', m, n, -1, -1 ) )
223  IF( nx.LT.k ) THEN
224 *
225 * Determine if workspace is large enough for blocked code.
226 *
227  ldwork = m
228  iws = ldwork*nb
229  IF( lwork.LT.iws ) THEN
230 *
231 * Not enough workspace to use optimal NB: reduce NB and
232 * determine the minimum value of NB.
233 *
234  nb = lwork / ldwork
235  nbmin = max( 2, ilaenv( 2, 'SGERQF', ' ', m, n, -1,
236  $ -1 ) )
237  END IF
238  END IF
239  END IF
240 *
241  IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
242 *
243 * Use blocked code initially.
244 * The last kk rows are handled by the block method.
245 *
246  ki = ( ( k-nx-1 ) / nb )*nb
247  kk = min( k, ki+nb )
248 *
249  DO 10 i = k - kk + ki + 1, k - kk + 1, -nb
250  ib = min( k-i+1, nb )
251 *
252 * Compute the RQ factorization of the current block
253 * A(m-k+i:m-k+i+ib-1,1:n-k+i+ib-1)
254 *
255  CALL sgerq2( ib, n-k+i+ib-1, a( m-k+i, 1 ), lda, tau( i ),
256  $ work, iinfo )
257  IF( m-k+i.GT.1 ) THEN
258 *
259 * Form the triangular factor of the block reflector
260 * H = H(i+ib-1) . . . H(i+1) H(i)
261 *
262  CALL slarft( 'Backward', 'Rowwise', n-k+i+ib-1, ib,
263  $ a( m-k+i, 1 ), lda, tau( i ), work, ldwork )
264 *
265 * Apply H to A(1:m-k+i-1,1:n-k+i+ib-1) from the right
266 *
267  CALL slarfb( 'Right', 'No transpose', 'Backward',
268  $ 'Rowwise', m-k+i-1, n-k+i+ib-1, ib,
269  $ a( m-k+i, 1 ), lda, work, ldwork, a, lda,
270  $ work( ib+1 ), ldwork )
271  END IF
272  10 CONTINUE
273  mu = m - k + i + nb - 1
274  nu = n - k + i + nb - 1
275  ELSE
276  mu = m
277  nu = n
278  END IF
279 *
280 * Use unblocked code to factor the last or only block
281 *
282  IF( mu.GT.0 .AND. nu.GT.0 )
283  $ CALL sgerq2( mu, nu, a, lda, tau, work, iinfo )
284 *
285  work( 1 ) = iws
286  RETURN
287 *
288 * End of SGERQF
289 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgerq2(M, N, A, LDA, TAU, WORK, INFO)
SGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: sgerq2.f:125
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine slarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: slarfb.f:197
subroutine slarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: slarft.f:165

Here is the call graph for this function:

Here is the caller graph for this function:

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

SGESVJ

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

Purpose:
 SGESVJ 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.
 SGESVJ 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=SQRT(M), EPS=SLAMCH('E').
          = 'C': Analogous to JOBU='U', except that user can control the
                 level of numerical orthogonality of the computed left
                 singular vectors. TOL can be set to TOL = CTOL*EPS, where
                 CTOL is given on input in the array WORK.
                 No CTOL smaller than ONE is allowed. CTOL greater
                 than 1 / EPS is meaningless. The option 'C'
                 can be used if M*EPS is satisfactory orthogonality
                 of the computed left singular vectors, so CTOL=M could
                 save few sweeps of Jacobi rotations.
                 See the descriptions of A and WORK(1).
          = 'N': The matrix U is not computed. However, see the
                 description of A.
[in]JOBV
          JOBV is CHARACTER*1
          Specifies whether to compute the right singular vectors, that
          is, the matrix V:
          = 'V' : the matrix V is computed and returned in the array V
          = 'A' : the Jacobi rotations are applied to the MV-by-N
                  array V. In other words, the right singular vector
                  matrix V is not computed explicitly; instead it is
                  applied to an MV-by-N matrix initially stored in the
                  first MV rows of V.
          = 'N' : the matrix V is not computed and the array V is not
                  referenced
[in]M
          M is INTEGER
          The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0.
[in]N
          N is INTEGER
          The number of columns of the input matrix A.
          M >= N >= 0.
[in,out]A
          A is REAL array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit,
          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C':
                 If INFO .EQ. 0 :
                 RANKA orthonormal columns of U are returned in the
                 leading RANKA columns of the array A. Here RANKA <= N
                 is the number of computed singular values of A that are
                 above the underflow threshold SLAMCH('S'). The singular
                 vectors corresponding to underflowed or zero singular
                 values are not computed. The value of RANKA is returned
                 in the array 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=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
                 see the description of JOBU.
                 If INFO .GT. 0,
                 the procedure SGESVJ 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 SGESVJ did not converge in the given number
                 of iterations (sweeps).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]SVA
          SVA is REAL array, dimension (N)
          On exit,
          If INFO .EQ. 0 :
          depending on the value SCALE = 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 SGESVJ 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 SGESVJ
          is applied to the first MV rows of V. See the description of JOBV.
[in,out]V
          V is REAL 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 REAL 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=SLAMCH('E').
                    It is required that CTOL >= ONE, i.e. it is not
                    allowed to force the routine to obtain orthogonality
                    below EPSILON.
          On exit,
          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
                    are the computed singular vcalues 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 SGESVJ 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 : SGESVJ 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.nosp@m.@mat.nosp@m.h.hr. Thank you.

Definition at line 325 of file sgesvj.f.

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