LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slaqp3rk()

subroutine slaqp3rk ( integer m,
integer n,
integer nrhs,
integer ioffset,
integer nb,
real abstol,
real reltol,
integer kp1,
real maxc2nrm,
real, dimension( lda, * ) a,
integer lda,
logical done,
integer kb,
real maxc2nrmk,
real relmaxc2nrmk,
integer, dimension( * ) jpiv,
real, dimension( * ) tau,
real, dimension( * ) vn1,
real, dimension( * ) vn2,
real, dimension( * ) auxv,
real, dimension( ldf, * ) f,
integer ldf,
integer, dimension( * ) iwork,
integer info )

SLAQP3RK computes a step of truncated QR factorization with column pivoting of a real m-by-n matrix A using Level 3 BLAS and overwrites a real m-by-nrhs matrix B with Q**T * B.

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

Purpose:
!>
!> SLAQP3RK computes a step of truncated QR factorization with column
!> pivoting of a real M-by-N matrix A block A(IOFFSET+1:M,1:N)
!> by using Level 3 BLAS as
!>
!>   A * P(KB) = Q(KB) * R(KB).
!>
!> The routine tries to factorize NB columns from A starting from
!> the row IOFFSET+1 and updates the residual matrix with BLAS 3
!> xGEMM. The number of actually factorized columns is returned
!> is smaller than NB.
!>
!> Block A(1:IOFFSET,1:N) is accordingly pivoted, but not factorized.
!>
!> The routine also overwrites the right-hand-sides B matrix stored
!> in A(IOFFSET+1:M,1:N+1:N+NRHS) with Q(KB)**T * B.
!>
!> Cases when the number of factorized columns KB < NB:
!>
!> (1) In some cases, due to catastrophic cancellations, it cannot
!> factorize all NB columns and need to update the residual matrix.
!> Hence, the actual number of factorized columns in the block returned
!> in KB is smaller than NB. The logical DONE is returned as FALSE.
!> The factorization of the whole original matrix A_orig must proceed
!> with the next block.
!>
!> (2) Whenever the stopping criterion ABSTOL or RELTOL is satisfied,
!> the factorization of the whole original matrix A_orig is stopped,
!> the logical DONE is returned as TRUE. The number of factorized
!> columns which is smaller than NB is returned in KB.
!>
!> (3) In case both stopping criteria ABSTOL or RELTOL are not used,
!> and when the residual matrix is a zero matrix in some factorization
!> step KB, the factorization of the whole original matrix A_orig is
!> stopped, the logical DONE is returned as TRUE. The number of
!> factorized columns which is smaller than NB is returned in KB.
!>
!> (4) Whenever NaN is detected in the matrix A or in the array TAU,
!> the factorization of the whole original matrix A_orig is stopped,
!> the logical DONE is returned as TRUE. The number of factorized
!> columns which is smaller than NB is returned in KB. The INFO
!> parameter is set to the column index of the first NaN occurrence.
!>
!>
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]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of
!>          columns of the matrix B. NRHS >= 0.
!>
[in]IOFFSET
!>          IOFFSET is INTEGER
!>          The number of rows of the matrix A that must be pivoted
!>          but not factorized. IOFFSET >= 0.
!>
!>          IOFFSET also represents the number of columns of the whole
!>          original matrix A_orig that have been factorized
!>          in the previous steps.
!>
[in]NB
!>          NB is INTEGER
!>          Factorization block size, i.e the number of columns
!>          to factorize in the matrix A. 0 <= NB
!>
!>          If NB = 0, then the routine exits immediately.
!>             This means that the factorization is not performed,
!>             the matrices A and B and the arrays TAU, IPIV
!>             are not modified.
!>
[in]ABSTOL
!>          ABSTOL is REAL, cannot be NaN.
!>
!>          The absolute tolerance (stopping threshold) for
!>          maximum column 2-norm of the residual matrix.
!>          The algorithm converges (stops the factorization) when
!>          the maximum column 2-norm of the residual matrix
!>          is less than or equal to ABSTOL.
!>
!>          a) If ABSTOL < 0.0, then this stopping criterion is not
!>                used, the routine factorizes columns depending
!>                on NB and RELTOL.
!>                This includes the case ABSTOL = -Inf.
!>
!>          b) If 0.0 <= ABSTOL then the input value
!>                of ABSTOL is used.
!>
[in]RELTOL
!>          RELTOL is REAL, cannot be NaN.
!>
!>          The tolerance (stopping threshold) for the ratio of the
!>          maximum column 2-norm of the residual matrix to the maximum
!>          column 2-norm of the original matrix A_orig. The algorithm
!>          converges (stops the factorization), when this ratio is
!>          less than or equal to RELTOL.
!>
!>          a) If RELTOL < 0.0, then this stopping criterion is not
!>                used, the routine factorizes columns depending
!>                on NB and ABSTOL.
!>                This includes the case RELTOL = -Inf.
!>
!>          d) If 0.0 <= RELTOL then the input value of RELTOL
!>                is used.
!>
[in]KP1
!>          KP1 is INTEGER
!>          The index of the column with the maximum 2-norm in
!>          the whole original matrix A_orig determined in the
!>          main routine SGEQP3RK. 1 <= KP1 <= N_orig.
!>
[in]MAXC2NRM
!>          MAXC2NRM is REAL
!>          The maximum column 2-norm of the whole original
!>          matrix A_orig computed in the main routine SGEQP3RK.
!>          MAXC2NRM >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA,N+NRHS)
!>          On entry:
!>              the M-by-N matrix A and M-by-NRHS matrix B, as in
!>
!>                                  N     NRHS
!>              array_A   =   M  [ mat_A, mat_B ]
!>
!>          On exit:
!>          1. The elements in block A(IOFFSET+1:M,1:KB) below
!>             the diagonal together with the array TAU represent
!>             the orthogonal matrix Q(KB) as a product of elementary
!>             reflectors.
!>          2. The upper triangular block of the matrix A stored
!>             in A(IOFFSET+1:M,1:KB) is the triangular factor obtained.
!>          3. The block of the matrix A stored in A(1:IOFFSET,1:N)
!>             has been accordingly pivoted, but not factorized.
!>          4. The rest of the array A, block A(IOFFSET+1:M,KB+1:N+NRHS).
!>             The left part A(IOFFSET+1:M,KB+1:N) of this block
!>             contains the residual of the matrix A, and,
!>             if NRHS > 0, the right part of the block
!>             A(IOFFSET+1:M,N+1:N+NRHS) contains the block of
!>             the right-hand-side matrix B. Both these blocks have been
!>             updated by multiplication from the left by Q(KB)**T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,M).
!>
[out]DONE
!>          DONE is LOGICAL
!>          TRUE: a) if the factorization completed before processing
!>                   all min(M-IOFFSET,NB,N) columns due to ABSTOL
!>                   or RELTOL criterion,
!>                b) if the factorization completed before processing
!>                   all min(M-IOFFSET,NB,N) columns due to the
!>                   residual matrix being a ZERO matrix.
!>                c) when NaN was detected in the matrix A
!>                   or in the array TAU.
!>          FALSE: otherwise.
!>
[out]KB
!>          KB is INTEGER
!>          Factorization rank of the matrix A, i.e. the rank of
!>          the factor R, which is the same as the number of non-zero
!>          rows of the factor R.  0 <= KB <= min(M-IOFFSET,NB,N).
!>
!>          KB also represents the number of non-zero Householder
!>          vectors.
!>
[out]MAXC2NRMK
!>          MAXC2NRMK is REAL
!>          The maximum column 2-norm of the residual matrix,
!>          when the factorization stopped at rank KB. MAXC2NRMK >= 0.
!>
[out]RELMAXC2NRMK
!>          RELMAXC2NRMK is REAL
!>          The ratio MAXC2NRMK / MAXC2NRM of the maximum column
!>          2-norm of the residual matrix (when the factorization
!>          stopped at rank KB) to the maximum column 2-norm of the
!>          original matrix A_orig. RELMAXC2NRMK >= 0.
!>
[out]JPIV
!>          JPIV is INTEGER array, dimension (N)
!>          Column pivot indices, for 1 <= j <= N, column j
!>          of the matrix A was interchanged with column JPIV(j).
!>
[out]TAU
!>          TAU is REAL array, dimension (min(M-IOFFSET,N))
!>          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.
!>
[out]AUXV
!>          AUXV is REAL array, dimension (NB)
!>          Auxiliary vector.
!>
[out]F
!>          F is REAL array, dimension (LDF,NB)
!>          Matrix F**T = L*(Y**T)*A.
!>
[in]LDF
!>          LDF is INTEGER
!>          The leading dimension of the array F. LDF >= max(1,N+NRHS).
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N-1).
!>          Is a work array. ( IWORK is used to store indices
!>          of  columns for norm downdating in the residual
!>          matrix ).
!>
[out]INFO
!>          INFO is INTEGER
!>          1) INFO = 0: successful exit.
!>          2) If INFO = j_1, where 1 <= j_1 <= N, then NaN was
!>             detected and the routine stops the computation.
!>             The j_1-th column of the matrix A or the j_1-th
!>             element of array TAU contains the first occurrence
!>             of NaN in the factorization step KB+1 ( when KB columns
!>             have been factorized ).
!>
!>             On exit:
!>             KB                  is set to the number of
!>                                    factorized columns without
!>                                    exception.
!>             MAXC2NRMK           is set to NaN.
!>             RELMAXC2NRMK        is set to NaN.
!>             TAU(KB+1:min(M,N))     is not set and contains undefined
!>                                    elements. If j_1=KB+1, TAU(KB+1)
!>                                    may contain NaN.
!>          3) If INFO = j_2, where N+1 <= j_2 <= 2*N, then no NaN
!>             was detected, but +Inf (or -Inf) was detected and
!>             the routine continues the computation until completion.
!>             The (j_2-N)-th column of the matrix A contains the first
!>             occurrence of +Inf (or -Inf) in the actorization
!>             step KB+1 ( when KB columns have been factorized ).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
References:
[1] A Level 3 BLAS QR factorization algorithm with column pivoting developed in 1996. G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain. X. Sun, Computer Science Dept., Duke University, USA. C. H. Bischof, Math. and Comp. Sci. Div., Argonne National Lab, USA. A BLAS-3 version of the QR factorization with column pivoting. LAPACK Working Note 114 https://www.netlib.org/lapack/lawnspdf/lawn114.pdf and in SIAM J. Sci. Comput., 19(5):1486-1494, Sept. 1998. https://doi.org/10.1137/S1064827595296732

[2] A partial column norm updating strategy developed in 2006. Z. Drmac and Z. Bujanovic, Dept. of Math., University of Zagreb, Croatia. On the failure of rank revealing QR factorization software – a case study. LAPACK Working Note 176. http://www.netlib.org/lapack/lawnspdf/lawn176.pdf and in ACM Trans. Math. Softw. 35, 2, Article 12 (July 2008), 28 pages. https://doi.org/10.1145/1377612.1377616

Contributors:
!>
!>  November  2023, Igor Kozachenko, James Demmel,
!>                  EECS Department,
!>                  University of California, Berkeley, USA.
!>
!>

Definition at line 388 of file slaqp3rk.f.

392 IMPLICIT NONE
393*
394* -- LAPACK auxiliary routine --
395* -- LAPACK is a software package provided by Univ. of Tennessee, --
396* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
397*
398* .. Scalar Arguments ..
399 LOGICAL DONE
400 INTEGER INFO, IOFFSET, KB, KP1, LDA, LDF, M, N,
401 $ NB, NRHS
402 REAL ABSTOL, MAXC2NRM, MAXC2NRMK, RELMAXC2NRMK,
403 $ RELTOL
404* ..
405* .. Array Arguments ..
406 INTEGER IWORK( * ), JPIV( * )
407 REAL A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ),
408 $ VN1( * ), VN2( * )
409* ..
410*
411* =====================================================================
412*
413* .. Parameters ..
414 REAL ZERO, ONE
415 parameter( zero = 0.0e+0, one = 1.0e+0 )
416* ..
417* .. Local Scalars ..
418 INTEGER ITEMP, J, K, MINMNFACT, MINMNUPDT,
419 $ LSTICC, KP, I, IF
420 REAL AIK, HUGEVAL, TEMP, TEMP2, TOL3Z
421* ..
422* .. External Subroutines ..
423 EXTERNAL sgemm, sgemv, slarfg, sswap
424* ..
425* .. Intrinsic Functions ..
426 INTRINSIC abs, max, min, sqrt
427* ..
428* .. External Functions ..
429 LOGICAL SISNAN
430 INTEGER ISAMAX
431 REAL SLAMCH, SNRM2
432 EXTERNAL sisnan, slamch, isamax, snrm2
433* ..
434* .. Executable Statements ..
435*
436* Initialize INFO
437*
438 info = 0
439*
440* MINMNFACT in the smallest dimension of the submatrix
441* A(IOFFSET+1:M,1:N) to be factorized.
442*
443 minmnfact = min( m-ioffset, n )
444 minmnupdt = min( m-ioffset, n+nrhs )
445 nb = min( nb, minmnfact )
446 tol3z = sqrt( slamch( 'Epsilon' ) )
447 hugeval = slamch( 'Overflow' )
448*
449* Compute factorization in a while loop over NB columns,
450* K is the column index in the block A(1:M,1:N).
451*
452 k = 0
453 lsticc = 0
454 done = .false.
455*
456 DO WHILE ( k.LT.nb .AND. lsticc.EQ.0 )
457 k = k + 1
458 i = ioffset + k
459*
460 IF( i.EQ.1 ) THEN
461*
462* We are at the first column of the original whole matrix A_orig,
463* therefore we use the computed KP1 and MAXC2NRM from the
464* main routine.
465*
466 kp = kp1
467*
468 ELSE
469*
470* Determine the pivot column in K-th step, i.e. the index
471* of the column with the maximum 2-norm in the
472* submatrix A(I:M,K:N).
473*
474 kp = ( k-1 ) + isamax( n-k+1, vn1( k ), 1 )
475*
476* Determine the maximum column 2-norm and the relative maximum
477* column 2-norm of the submatrix A(I:M,K:N) in step K.
478*
479 maxc2nrmk = vn1( kp )
480*
481* ============================================================
482*
483* Check if the submatrix A(I:M,K:N) contains NaN, set
484* INFO parameter to the column number, where the first NaN
485* is found and return from the routine.
486* We need to check the condition only if the
487* column index (same as row index) of the original whole
488* matrix is larger than 1, since the condition for whole
489* original matrix is checked in the main routine.
490*
491 IF( sisnan( maxc2nrmk ) ) THEN
492*
493 done = .true.
494*
495* Set KB, the number of factorized partial columns
496* that are non-zero in each step in the block,
497* i.e. the rank of the factor R.
498* Set IF, the number of processed rows in the block, which
499* is the same as the number of processed rows in
500* the original whole matrix A_orig.
501*
502 kb = k - 1
503 IF = i - 1
504 info = kb + kp
505*
506* Set RELMAXC2NRMK to NaN.
507*
508 relmaxc2nrmk = maxc2nrmk
509*
510* There is no need to apply the block reflector to the
511* residual of the matrix A stored in A(KB+1:M,KB+1:N),
512* since the submatrix contains NaN and we stop
513* the computation.
514* But, we need to apply the block reflector to the residual
515* right hand sides stored in A(KB+1:M,N+1:N+NRHS), if the
516* residual right hand sides exist. This occurs
517* when ( NRHS != 0 AND KB <= (M-IOFFSET) ):
518*
519* A(I+1:M,N+1:N+NRHS) := A(I+1:M,N+1:N+NRHS) -
520* A(I+1:M,1:KB) * F(N+1:N+NRHS,1:KB)**T.
521
522 IF( nrhs.GT.0 .AND. kb.LT.(m-ioffset) ) THEN
523 CALL sgemm( 'No transpose', 'Transpose',
524 $ m-IF, nrhs, kb, -one, a( if+1, 1 ), lda,
525 $ f( n+1, 1 ), ldf, one, a( if+1, n+1 ), lda )
526 END IF
527*
528* There is no need to recompute the 2-norm of the
529* difficult columns, since we stop the factorization.
530*
531* Array TAU(KF+1:MINMNFACT) is not set and contains
532* undefined elements.
533*
534* Return from the routine.
535*
536 RETURN
537 END IF
538*
539* Quick return, if the submatrix A(I:M,K:N) is
540* a zero matrix. We need to check it only if the column index
541* (same as row index) is larger than 1, since the condition
542* for the whole original matrix A_orig is checked in the main
543* routine.
544*
545 IF( maxc2nrmk.EQ.zero ) THEN
546*
547 done = .true.
548*
549* Set KB, the number of factorized partial columns
550* that are non-zero in each step in the block,
551* i.e. the rank of the factor R.
552* Set IF, the number of processed rows in the block, which
553* is the same as the number of processed rows in
554* the original whole matrix A_orig.
555*
556 kb = k - 1
557 IF = i - 1
558 relmaxc2nrmk = zero
559*
560* There is no need to apply the block reflector to the
561* residual of the matrix A stored in A(KB+1:M,KB+1:N),
562* since the submatrix is zero and we stop the computation.
563* But, we need to apply the block reflector to the residual
564* right hand sides stored in A(KB+1:M,N+1:N+NRHS), if the
565* residual right hand sides exist. This occurs
566* when ( NRHS != 0 AND KB <= (M-IOFFSET) ):
567*
568* A(I+1:M,N+1:N+NRHS) := A(I+1:M,N+1:N+NRHS) -
569* A(I+1:M,1:KB) * F(N+1:N+NRHS,1:KB)**T.
570*
571 IF( nrhs.GT.0 .AND. kb.LT.(m-ioffset) ) THEN
572 CALL sgemm( 'No transpose', 'Transpose',
573 $ m-IF, nrhs, kb, -one, a( if+1, 1 ), lda,
574 $ f( n+1, 1 ), ldf, one, a( if+1, n+1 ), lda )
575 END IF
576*
577* There is no need to recompute the 2-norm of the
578* difficult columns, since we stop the factorization.
579*
580* Set TAUs corresponding to the columns that were not
581* factorized to ZERO, i.e. set TAU(KB+1:MINMNFACT) = ZERO,
582* which is equivalent to seting TAU(K:MINMNFACT) = ZERO.
583*
584 DO j = k, minmnfact
585 tau( j ) = zero
586 END DO
587*
588* Return from the routine.
589*
590 RETURN
591*
592 END IF
593*
594* ============================================================
595*
596* Check if the submatrix A(I:M,K:N) contains Inf,
597* set INFO parameter to the column number, where
598* the first Inf is found plus N, and continue
599* the computation.
600* We need to check the condition only if the
601* column index (same as row index) of the original whole
602* matrix is larger than 1, since the condition for whole
603* original matrix is checked in the main routine.
604*
605 IF( info.EQ.0 .AND. maxc2nrmk.GT.hugeval ) THEN
606 info = n + k - 1 + kp
607 END IF
608*
609* ============================================================
610*
611* Test for the second and third tolerance stopping criteria.
612* NOTE: There is no need to test for ABSTOL.GE.ZERO, since
613* MAXC2NRMK is non-negative. Similarly, there is no need
614* to test for RELTOL.GE.ZERO, since RELMAXC2NRMK is
615* non-negative.
616* We need to check the condition only if the
617* column index (same as row index) of the original whole
618* matrix is larger than 1, since the condition for whole
619* original matrix is checked in the main routine.
620*
621 relmaxc2nrmk = maxc2nrmk / maxc2nrm
622*
623 IF( maxc2nrmk.LE.abstol .OR. relmaxc2nrmk.LE.reltol ) THEN
624*
625 done = .true.
626*
627* Set KB, the number of factorized partial columns
628* that are non-zero in each step in the block,
629* i.e. the rank of the factor R.
630* Set IF, the number of processed rows in the block, which
631* is the same as the number of processed rows in
632* the original whole matrix A_orig;
633*
634 kb = k - 1
635 IF = i - 1
636*
637* Apply the block reflector to the residual of the
638* matrix A and the residual of the right hand sides B, if
639* the residual matrix and and/or the residual of the right
640* hand sides exist, i.e. if the submatrix
641* A(I+1:M,KB+1:N+NRHS) exists. This occurs when
642* KB < MINMNUPDT = min( M-IOFFSET, N+NRHS ):
643*
644* A(IF+1:M,K+1:N+NRHS) := A(IF+1:M,KB+1:N+NRHS) -
645* A(IF+1:M,1:KB) * F(KB+1:N+NRHS,1:KB)**T.
646*
647 IF( kb.LT.minmnupdt ) THEN
648 CALL sgemm( 'No transpose', 'Transpose',
649 $ m-IF, n+nrhs-kb, kb,-one, a( if+1, 1 ), lda,
650 $ f( kb+1, 1 ), ldf, one, a( if+1, kb+1 ), lda )
651 END IF
652*
653* There is no need to recompute the 2-norm of the
654* difficult columns, since we stop the factorization.
655*
656* Set TAUs corresponding to the columns that were not
657* factorized to ZERO, i.e. set TAU(KB+1:MINMNFACT) = ZERO,
658* which is equivalent to seting TAU(K:MINMNFACT) = ZERO.
659*
660 DO j = k, minmnfact
661 tau( j ) = zero
662 END DO
663*
664* Return from the routine.
665*
666 RETURN
667*
668 END IF
669*
670* ============================================================
671*
672* End ELSE of IF(I.EQ.1)
673*
674 END IF
675*
676* ===============================================================
677*
678* If the pivot column is not the first column of the
679* subblock A(1:M,K:N):
680* 1) swap the K-th column and the KP-th pivot column
681* in A(1:M,1:N);
682* 2) swap the K-th row and the KP-th row in F(1:N,1:K-1)
683* 3) copy the K-th element into the KP-th element of the partial
684* and exact 2-norm vectors VN1 and VN2. (Swap is not needed
685* for VN1 and VN2 since we use the element with the index
686* larger than K in the next loop step.)
687* 4) Save the pivot interchange with the indices relative to the
688* the original matrix A_orig, not the block A(1:M,1:N).
689*
690 IF( kp.NE.k ) THEN
691 CALL sswap( m, a( 1, kp ), 1, a( 1, k ), 1 )
692 CALL sswap( k-1, f( kp, 1 ), ldf, f( k, 1 ), ldf )
693 vn1( kp ) = vn1( k )
694 vn2( kp ) = vn2( k )
695 itemp = jpiv( kp )
696 jpiv( kp ) = jpiv( k )
697 jpiv( k ) = itemp
698 END IF
699*
700* Apply previous Householder reflectors to column K:
701* A(I:M,K) := A(I:M,K) - A(I:M,1:K-1)*F(K,1:K-1)**T.
702*
703 IF( k.GT.1 ) THEN
704 CALL sgemv( 'No transpose', m-i+1, k-1, -one, a( i, 1 ),
705 $ lda, f( k, 1 ), ldf, one, a( i, k ), 1 )
706 END IF
707*
708* Generate elementary reflector H(k) using the column A(I:M,K).
709*
710 IF( i.LT.m ) THEN
711 CALL slarfg( m-i+1, a( i, k ), a( i+1, k ), 1, tau( k ) )
712 ELSE
713 tau( k ) = zero
714 END IF
715*
716* Check if TAU(K) contains NaN, set INFO parameter
717* to the column number where NaN is found and return from
718* the routine.
719* NOTE: There is no need to check TAU(K) for Inf,
720* since SLARFG cannot produce TAU(K) or Householder vector
721* below the diagonal containing Inf. Only BETA on the diagonal,
722* returned by SLARFG can contain Inf, which requires
723* TAU(K) to contain NaN. Therefore, this case of generating Inf
724* by SLARFG is covered by checking TAU(K) for NaN.
725*
726 IF( sisnan( tau(k) ) ) THEN
727*
728 done = .true.
729*
730* Set KB, the number of factorized partial columns
731* that are non-zero in each step in the block,
732* i.e. the rank of the factor R.
733* Set IF, the number of processed rows in the block, which
734* is the same as the number of processed rows in
735* the original whole matrix A_orig.
736*
737 kb = k - 1
738 IF = i - 1
739 info = k
740*
741* Set MAXC2NRMK and RELMAXC2NRMK to NaN.
742*
743 maxc2nrmk = tau( k )
744 relmaxc2nrmk = tau( k )
745*
746* There is no need to apply the block reflector to the
747* residual of the matrix A stored in A(KB+1:M,KB+1:N),
748* since the submatrix contains NaN and we stop
749* the computation.
750* But, we need to apply the block reflector to the residual
751* right hand sides stored in A(KB+1:M,N+1:N+NRHS), if the
752* residual right hand sides exist. This occurs
753* when ( NRHS != 0 AND KB <= (M-IOFFSET) ):
754*
755* A(I+1:M,N+1:N+NRHS) := A(I+1:M,N+1:N+NRHS) -
756* A(I+1:M,1:KB) * F(N+1:N+NRHS,1:KB)**T.
757*
758 IF( nrhs.GT.0 .AND. kb.LT.(m-ioffset) ) THEN
759 CALL sgemm( 'No transpose', 'Transpose',
760 $ m-IF, nrhs, kb, -one, a( if+1, 1 ), lda,
761 $ f( n+1, 1 ), ldf, one, a( if+1, n+1 ), lda )
762 END IF
763*
764* There is no need to recompute the 2-norm of the
765* difficult columns, since we stop the factorization.
766*
767* Array TAU(KF+1:MINMNFACT) is not set and contains
768* undefined elements.
769*
770* Return from the routine.
771*
772 RETURN
773 END IF
774*
775* ===============================================================
776*
777 aik = a( i, k )
778 a( i, k ) = one
779*
780* ===============================================================
781*
782* Compute the current K-th column of F:
783* 1) F(K+1:N,K) := tau(K) * A(I:M,K+1:N)**T * A(I:M,K).
784*
785 IF( k.LT.n+nrhs ) THEN
786 CALL sgemv( 'Transpose', m-i+1, n+nrhs-k,
787 $ tau( k ), a( i, k+1 ), lda, a( i, k ), 1,
788 $ zero, f( k+1, k ), 1 )
789 END IF
790*
791* 2) Zero out elements above and on the diagonal of the
792* column K in matrix F, i.e elements F(1:K,K).
793*
794 DO j = 1, k
795 f( j, k ) = zero
796 END DO
797*
798* 3) Incremental updating of the K-th column of F:
799* F(1:N,K) := F(1:N,K) - tau(K) * F(1:N,1:K-1) * A(I:M,1:K-1)**T
800* * A(I:M,K).
801*
802 IF( k.GT.1 ) THEN
803 CALL sgemv( 'Transpose', m-i+1, k-1, -tau( k ),
804 $ a( i, 1 ), lda, a( i, k ), 1, zero,
805 $ auxv( 1 ), 1 )
806*
807 CALL sgemv( 'No transpose', n+nrhs, k-1, one,
808 $ f( 1, 1 ), ldf, auxv( 1 ), 1, one,
809 $ f( 1, k ), 1 )
810 END IF
811*
812* ===============================================================
813*
814* Update the current I-th row of A:
815* A(I,K+1:N+NRHS) := A(I,K+1:N+NRHS)
816* - A(I,1:K)*F(K+1:N+NRHS,1:K)**T.
817*
818 IF( k.LT.n+nrhs ) THEN
819 CALL sgemv( 'No transpose', n+nrhs-k, k, -one,
820 $ f( k+1, 1 ), ldf, a( i, 1 ), lda, one,
821 $ a( i, k+1 ), lda )
822 END IF
823*
824 a( i, k ) = aik
825*
826* Update the partial column 2-norms for the residual matrix,
827* only if the residual matrix A(I+1:M,K+1:N) exists, i.e.
828* when K < MINMNFACT = min( M-IOFFSET, N ).
829*
830 IF( k.LT.minmnfact ) THEN
831*
832 DO j = k + 1, n
833 IF( vn1( j ).NE.zero ) THEN
834*
835* NOTE: The following lines follow from the analysis in
836* Lapack Working Note 176.
837*
838 temp = abs( a( i, j ) ) / vn1( j )
839 temp = max( zero, ( one+temp )*( one-temp ) )
840 temp2 = temp*( vn1( j ) / vn2( j ) )**2
841 IF( temp2.LE.tol3z ) THEN
842*
843* At J-index, we have a difficult column for the
844* update of the 2-norm. Save the index of the previous
845* difficult column in IWORK(J-1).
846* NOTE: ILSTCC > 1, threfore we can use IWORK only
847* with N-1 elements, where the elements are
848* shifted by 1 to the left.
849*
850 iwork( j-1 ) = lsticc
851*
852* Set the index of the last difficult column LSTICC.
853*
854 lsticc = j
855*
856 ELSE
857 vn1( j ) = vn1( j )*sqrt( temp )
858 END IF
859 END IF
860 END DO
861*
862 END IF
863*
864* End of while loop.
865*
866 END DO
867*
868* Now, afler the loop:
869* Set KB, the number of factorized columns in the block;
870* Set IF, the number of processed rows in the block, which
871* is the same as the number of processed rows in
872* the original whole matrix A_orig, IF = IOFFSET + KB.
873*
874 kb = k
875 IF = i
876*
877* Apply the block reflector to the residual of the matrix A
878* and the residual of the right hand sides B, if the residual
879* matrix and and/or the residual of the right hand sides
880* exist, i.e. if the submatrix A(I+1:M,KB+1:N+NRHS) exists.
881* This occurs when KB < MINMNUPDT = min( M-IOFFSET, N+NRHS ):
882*
883* A(IF+1:M,K+1:N+NRHS) := A(IF+1:M,KB+1:N+NRHS) -
884* A(IF+1:M,1:KB) * F(KB+1:N+NRHS,1:KB)**T.
885*
886 IF( kb.LT.minmnupdt ) THEN
887 CALL sgemm( 'No transpose', 'Transpose',
888 $ m-IF, n+nrhs-kb, kb, -one, a( if+1, 1 ), lda,
889 $ f( kb+1, 1 ), ldf, one, a( if+1, kb+1 ), lda )
890 END IF
891*
892* Recompute the 2-norm of the difficult columns.
893* Loop over the index of the difficult columns from the largest
894* to the smallest index.
895*
896 DO WHILE( lsticc.GT.0 )
897*
898* LSTICC is the index of the last difficult column is greater
899* than 1.
900* ITEMP is the index of the previous difficult column.
901*
902 itemp = iwork( lsticc-1 )
903*
904* Compute the 2-norm explicilty for the last difficult column and
905* save it in the partial and exact 2-norm vectors VN1 and VN2.
906*
907* NOTE: The computation of VN1( LSTICC ) relies on the fact that
908* SNRM2 does not fail on vectors with norm below the value of
909* SQRT(SLAMCH('S'))
910*
911 vn1( lsticc ) = snrm2( m-IF, a( if+1, lsticc ), 1 )
912 vn2( lsticc ) = vn1( lsticc )
913*
914* Downdate the index of the last difficult column to
915* the index of the previous difficult column.
916*
917 lsticc = itemp
918*
919 END DO
920*
921 RETURN
922*
923* End of SLAQP3RK
924*
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:57
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slarfg
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:104
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: