LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ claqps()

subroutine claqps ( integer  M,
integer  N,
integer  OFFSET,
integer  NB,
integer  KB,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
complex, dimension( * )  TAU,
real, dimension( * )  VN1,
real, dimension( * )  VN2,
complex, dimension( * )  AUXV,
complex, dimension( ldf, * )  F,
integer  LDF 
)

CLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3.

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

Purpose:
 CLAQPS computes a step of QR factorization with column pivoting
 of a complex M-by-N matrix A by using Blas-3.  It tries to factorize
 NB columns from A starting from the row OFFSET+1, and updates all
 of the matrix with Blas-3 xGEMM.

 In some cases, due to catastrophic cancellations, it cannot
 factorize NB columns.  Hence, the actual number of factorized
 columns is returned in KB.

 Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
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]OFFSET
          OFFSET is INTEGER
          The number of rows of A that have been factorized in
          previous steps.
[in]NB
          NB is INTEGER
          The number of columns to factorize.
[out]KB
          KB is INTEGER
          The number of columns actually factorized.
[in,out]A
          A is COMPLEX array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, block A(OFFSET+1:M,1:KB) is the triangular
          factor obtained and block A(1:OFFSET,1:N) has been
          accordingly pivoted, but no factorized.
          The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has
          been updated.
[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)
          JPVT(I) = K <==> Column K of the full matrix A has been
          permuted into position I in AP.
[out]TAU
          TAU is COMPLEX array, dimension (KB)
          The scalar factors of the elementary reflectors.
[in,out]VN1
          VN1 is REAL array, dimension (N)
          The vector with the partial column norms.
[in,out]VN2
          VN2 is REAL array, dimension (N)
          The vector with the exact column norms.
[in,out]AUXV
          AUXV is COMPLEX array, dimension (NB)
          Auxiliary vector.
[in,out]F
          F is COMPLEX array, dimension (LDF,NB)
          Matrix  F**H = L * Y**H * A.
[in]LDF
          LDF is INTEGER
          The leading dimension of the array F. LDF >= max(1,N).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain X. Sun, Computer Science Dept., Duke University, USA


Partial column norm updating strategy modified on April 2011 Z. Drmac and Z. Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.

References:
LAPACK Working Note 176 [PDF]

Definition at line 176 of file claqps.f.

178 *
179 * -- LAPACK auxiliary routine --
180 * -- LAPACK is a software package provided by Univ. of Tennessee, --
181 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182 *
183 * .. Scalar Arguments ..
184  INTEGER KB, LDA, LDF, M, N, NB, OFFSET
185 * ..
186 * .. Array Arguments ..
187  INTEGER JPVT( * )
188  REAL VN1( * ), VN2( * )
189  COMPLEX A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  REAL ZERO, ONE
196  COMPLEX CZERO, CONE
197  parameter( zero = 0.0e+0, one = 1.0e+0,
198  $ czero = ( 0.0e+0, 0.0e+0 ),
199  $ cone = ( 1.0e+0, 0.0e+0 ) )
200 * ..
201 * .. Local Scalars ..
202  INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK
203  REAL TEMP, TEMP2, TOL3Z
204  COMPLEX AKK
205 * ..
206 * .. External Subroutines ..
207  EXTERNAL cgemm, cgemv, clarfg, cswap
208 * ..
209 * .. Intrinsic Functions ..
210  INTRINSIC abs, conjg, max, min, nint, real, sqrt
211 * ..
212 * .. External Functions ..
213  INTEGER ISAMAX
214  REAL SCNRM2, SLAMCH
215  EXTERNAL isamax, scnrm2, slamch
216 * ..
217 * .. Executable Statements ..
218 *
219  lastrk = min( m, n+offset )
220  lsticc = 0
221  k = 0
222  tol3z = sqrt(slamch('Epsilon'))
223 *
224 * Beginning of while loop.
225 *
226  10 CONTINUE
227  IF( ( k.LT.nb ) .AND. ( lsticc.EQ.0 ) ) THEN
228  k = k + 1
229  rk = offset + k
230 *
231 * Determine ith pivot column and swap if necessary
232 *
233  pvt = ( k-1 ) + isamax( n-k+1, vn1( k ), 1 )
234  IF( pvt.NE.k ) THEN
235  CALL cswap( m, a( 1, pvt ), 1, a( 1, k ), 1 )
236  CALL cswap( k-1, f( pvt, 1 ), ldf, f( k, 1 ), ldf )
237  itemp = jpvt( pvt )
238  jpvt( pvt ) = jpvt( k )
239  jpvt( k ) = itemp
240  vn1( pvt ) = vn1( k )
241  vn2( pvt ) = vn2( k )
242  END IF
243 *
244 * Apply previous Householder reflectors to column K:
245 * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**H.
246 *
247  IF( k.GT.1 ) THEN
248  DO 20 j = 1, k - 1
249  f( k, j ) = conjg( f( k, j ) )
250  20 CONTINUE
251  CALL cgemv( 'No transpose', m-rk+1, k-1, -cone, a( rk, 1 ),
252  $ lda, f( k, 1 ), ldf, cone, a( rk, k ), 1 )
253  DO 30 j = 1, k - 1
254  f( k, j ) = conjg( f( k, j ) )
255  30 CONTINUE
256  END IF
257 *
258 * Generate elementary reflector H(k).
259 *
260  IF( rk.LT.m ) THEN
261  CALL clarfg( m-rk+1, a( rk, k ), a( rk+1, k ), 1, tau( k ) )
262  ELSE
263  CALL clarfg( 1, a( rk, k ), a( rk, k ), 1, tau( k ) )
264  END IF
265 *
266  akk = a( rk, k )
267  a( rk, k ) = cone
268 *
269 * Compute Kth column of F:
270 *
271 * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K).
272 *
273  IF( k.LT.n ) THEN
274  CALL cgemv( 'Conjugate transpose', m-rk+1, n-k, tau( k ),
275  $ a( rk, k+1 ), lda, a( rk, k ), 1, czero,
276  $ f( k+1, k ), 1 )
277  END IF
278 *
279 * Padding F(1:K,K) with zeros.
280 *
281  DO 40 j = 1, k
282  f( j, k ) = czero
283  40 CONTINUE
284 *
285 * Incremental updating of F:
286 * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**H
287 * *A(RK:M,K).
288 *
289  IF( k.GT.1 ) THEN
290  CALL cgemv( 'Conjugate transpose', m-rk+1, k-1, -tau( k ),
291  $ a( rk, 1 ), lda, a( rk, k ), 1, czero,
292  $ auxv( 1 ), 1 )
293 *
294  CALL cgemv( 'No transpose', n, k-1, cone, f( 1, 1 ), ldf,
295  $ auxv( 1 ), 1, cone, f( 1, k ), 1 )
296  END IF
297 *
298 * Update the current row of A:
299 * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**H.
300 *
301  IF( k.LT.n ) THEN
302  CALL cgemm( 'No transpose', 'Conjugate transpose', 1, n-k,
303  $ k, -cone, a( rk, 1 ), lda, f( k+1, 1 ), ldf,
304  $ cone, a( rk, k+1 ), lda )
305  END IF
306 *
307 * Update partial column norms.
308 *
309  IF( rk.LT.lastrk ) THEN
310  DO 50 j = k + 1, n
311  IF( vn1( j ).NE.zero ) THEN
312 *
313 * NOTE: The following 4 lines follow from the analysis in
314 * Lapack Working Note 176.
315 *
316  temp = abs( a( rk, j ) ) / vn1( j )
317  temp = max( zero, ( one+temp )*( one-temp ) )
318  temp2 = temp*( vn1( j ) / vn2( j ) )**2
319  IF( temp2 .LE. tol3z ) THEN
320  vn2( j ) = real( lsticc )
321  lsticc = j
322  ELSE
323  vn1( j ) = vn1( j )*sqrt( temp )
324  END IF
325  END IF
326  50 CONTINUE
327  END IF
328 *
329  a( rk, k ) = akk
330 *
331 * End of while loop.
332 *
333  GO TO 10
334  END IF
335  kb = k
336  rk = offset + kb
337 *
338 * Apply the block reflector to the rest of the matrix:
339 * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -
340 * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**H.
341 *
342  IF( kb.LT.min( n, m-offset ) ) THEN
343  CALL cgemm( 'No transpose', 'Conjugate transpose', m-rk, n-kb,
344  $ kb, -cone, a( rk+1, 1 ), lda, f( kb+1, 1 ), ldf,
345  $ cone, a( rk+1, kb+1 ), lda )
346  END IF
347 *
348 * Recomputation of difficult columns.
349 *
350  60 CONTINUE
351  IF( lsticc.GT.0 ) THEN
352  itemp = nint( vn2( lsticc ) )
353  vn1( lsticc ) = scnrm2( m-rk, a( rk+1, lsticc ), 1 )
354 *
355 * NOTE: The computation of VN1( LSTICC ) relies on the fact that
356 * SNRM2 does not fail on vectors with norm below the value of
357 * SQRT(DLAMCH('S'))
358 *
359  vn2( lsticc ) = vn1( lsticc )
360  lsticc = itemp
361  GO TO 60
362  END IF
363 *
364  RETURN
365 *
366 * End of CLAQPS
367 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:106
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition: scnrm2.f90:90
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: