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

◆ slaqp2rk()

subroutine slaqp2rk ( integer  m,
integer  n,
integer  nrhs,
integer  ioffset,
integer  kmax,
real  abstol,
real  reltol,
integer  kp1,
real  maxc2nrm,
real, dimension( lda, * )  a,
integer  lda,
integer  k,
real  maxc2nrmk,
real  relmaxc2nrmk,
integer, dimension( * )  jpiv,
real, dimension( * )  tau,
real, dimension( * )  vn1,
real, dimension( * )  vn2,
real, dimension( * )  work,
integer  info 
)

SLAQP2RK computes truncated QR factorization with column pivoting of a real matrix block using Level 2 BLAS and overwrites a real m-by-nrhs matrix B with Q**T * B.

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

Purpose:
 SLAQP2RK computes a truncated (rank K) or full rank Householder QR
 factorization with column pivoting of a real matrix
 block A(IOFFSET+1:M,1:N) as

   A * P(K) = Q(K) * R(K).

 The routine uses Level 2 BLAS. The block A(1:IOFFSET,1:N)
 is accordingly pivoted, but not factorized.

 The routine also overwrites the right-hand-sides matrix block B
 stored in A(IOFFSET+1:M,N+1:N+NRHS) with Q(K)**T * B.
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]KMAX
          KMAX is INTEGER

          The first factorization stopping criterion. KMAX >= 0.

          The maximum number of columns of the matrix A to factorize,
          i.e. the maximum factorization rank.

          a) If KMAX >= min(M-IOFFSET,N), then this stopping
                criterion is not used, factorize columns
                depending on ABSTOL and RELTOL.

          b) If KMAX = 0, then this stopping criterion is
             satisfied on input and 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 DOUBLE PRECISION, cannot be NaN.

          The second factorization stopping criterion.

          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 KMAX 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 DOUBLE PRECISION, cannot be NaN.

          The third factorization stopping criterion.

          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 KMAX 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_mat.
[in]MAXC2NRM
          MAXC2NRM is DOUBLE PRECISION
          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:K) below
             the diagonal together with the array TAU represent
             the orthogonal matrix Q(K) as a product of elementary
             reflectors.
          2. The upper triangular block of the matrix A stored
             in A(IOFFSET+1:M,1:K) 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,K+1:N+NRHS).
             The left part A(IOFFSET+1:M,K+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(K)**T.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[out]K
          K 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 <= K <= min(M-IOFFSET,KMAX,N).

          K also represents the number of non-zero Householder
          vectors.
[out]MAXC2NRMK
          MAXC2NRMK is DOUBLE PRECISION
          The maximum column 2-norm of the residual matrix,
          when the factorization stopped at rank K. MAXC2NRMK >= 0.
[out]RELMAXC2NRMK
          RELMAXC2NRMK is DOUBLE PRECISION
          The ratio MAXC2NRMK / MAXC2NRM of the maximum column
          2-norm of the residual matrix (when the factorization
          stopped at rank K) to the maximum column 2-norm of the
          whole original matrix A. 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]WORK
          WORK is REAL array, dimension (N-1)
          Used in SLARF subroutine to apply an elementary
          reflector from the left.
[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 K+1 ( when K columns
             have been factorized ).

             On exit:
             K                  is set to the number of
                                   factorized columns without
                                   exception.
             MAXC2NRMK          is set to NaN.
             RELMAXC2NRMK       is set to NaN.
             TAU(K+1:min(M,N))  is not set and contains undefined
                                   elements. If j_1=K+1, TAU(K+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 factorization
             step K+1 ( when K 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 340 of file slaqp2rk.f.

344 IMPLICIT NONE
345*
346* -- LAPACK auxiliary routine --
347* -- LAPACK is a software package provided by Univ. of Tennessee, --
348* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
349*
350* .. Scalar Arguments ..
351 INTEGER INFO, IOFFSET, KP1, K, KMAX, LDA, M, N, NRHS
352 REAL ABSTOL, MAXC2NRM, MAXC2NRMK, RELMAXC2NRMK,
353 $ RELTOL
354* ..
355* .. Array Arguments ..
356 INTEGER JPIV( * )
357 REAL A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
358 $ WORK( * )
359* ..
360*
361* =====================================================================
362*
363* .. Parameters ..
364 REAL ZERO, ONE
365 parameter( zero = 0.0e+0, one = 1.0e+0 )
366* ..
367* .. Local Scalars ..
368 INTEGER I, ITEMP, J, JMAXC2NRM, KK, KP, MINMNFACT,
369 $ MINMNUPDT
370 REAL AIKK, HUGEVAL, TEMP, TEMP2, TOL3Z
371* ..
372* .. External Subroutines ..
373 EXTERNAL slarf, slarfg, sswap
374* ..
375* .. Intrinsic Functions ..
376 INTRINSIC abs, max, min, sqrt
377* ..
378* .. External Functions ..
379 LOGICAL SISNAN
380 INTEGER ISAMAX
381 REAL SLAMCH, SNRM2
382 EXTERNAL sisnan, slamch, isamax, snrm2
383* ..
384* .. Executable Statements ..
385*
386* Initialize INFO
387*
388 info = 0
389*
390* MINMNFACT in the smallest dimension of the submatrix
391* A(IOFFSET+1:M,1:N) to be factorized.
392*
393* MINMNUPDT is the smallest dimension
394* of the subarray A(IOFFSET+1:M,1:N+NRHS) to be udated, which
395* contains the submatrices A(IOFFSET+1:M,1:N) and
396* B(IOFFSET+1:M,1:NRHS) as column blocks.
397*
398 minmnfact = min( m-ioffset, n )
399 minmnupdt = min( m-ioffset, n+nrhs )
400 kmax = min( kmax, minmnfact )
401 tol3z = sqrt( slamch( 'Epsilon' ) )
402 hugeval = slamch( 'Overflow' )
403*
404* Compute the factorization, KK is the lomn loop index.
405*
406 DO kk = 1, kmax
407*
408 i = ioffset + kk
409*
410 IF( i.EQ.1 ) THEN
411*
412* ============================================================
413*
414* We are at the first column of the original whole matrix A,
415* therefore we use the computed KP1 and MAXC2NRM from the
416* main routine.
417*
418
419 kp = kp1
420*
421* ============================================================
422*
423 ELSE
424*
425* ============================================================
426*
427* Determine the pivot column in KK-th step, i.e. the index
428* of the column with the maximum 2-norm in the
429* submatrix A(I:M,K:N).
430*
431 kp = ( kk-1 ) + isamax( n-kk+1, vn1( kk ), 1 )
432*
433* Determine the maximum column 2-norm and the relative maximum
434* column 2-norm of the submatrix A(I:M,KK:N) in step KK.
435* RELMAXC2NRMK will be computed later, after somecondition
436* checks on MAXC2NRMK.
437*
438 maxc2nrmk = vn1( kp )
439*
440* ============================================================
441*
442* Check if the submatrix A(I:M,KK:N) contains NaN, and set
443* INFO parameter to the column number, where the first NaN
444* is found and return from the routine.
445* We need to check the condition only if the
446* column index (same as row index) of the original whole
447* matrix is larger than 1, since the condition for whole
448* original matrix is checked in the main routine.
449*
450 IF( sisnan( maxc2nrmk ) ) THEN
451*
452* Set K, the number of factorized columns.
453* that are not zero.
454*
455 k = kk - 1
456 info = k + kp
457*
458* Set RELMAXC2NRMK to NaN.
459*
460 relmaxc2nrmk = maxc2nrmk
461*
462* Array TAU(K+1:MINMNFACT) is not set and contains
463* undefined elements.
464*
465 RETURN
466 END IF
467*
468* ============================================================
469*
470* Quick return, if the submatrix A(I:M,KK:N) is
471* a zero matrix.
472* We need to check the condition only if the
473* column index (same as row index) of the original whole
474* matrix is larger than 1, since the condition for whole
475* original matrix is checked in the main routine.
476*
477 IF( maxc2nrmk.EQ.zero ) THEN
478*
479* Set K, the number of factorized columns.
480* that are not zero.
481*
482 k = kk - 1
483 relmaxc2nrmk = zero
484*
485* Set TAUs corresponding to the columns that were not
486* factorized to ZERO, i.e. set TAU(KK:MINMNFACT) to ZERO.
487*
488 DO j = kk, minmnfact
489 tau( j ) = zero
490 END DO
491*
492* Return from the routine.
493*
494 RETURN
495*
496 END IF
497*
498* ============================================================
499*
500* Check if the submatrix A(I:M,KK:N) contains Inf,
501* set INFO parameter to the column number, where
502* the first Inf is found plus N, and continue
503* the computation.
504* We need to check the condition only if the
505* column index (same as row index) of the original whole
506* matrix is larger than 1, since the condition for whole
507* original matrix is checked in the main routine.
508*
509 IF( info.EQ.0 .AND. maxc2nrmk.GT.hugeval ) THEN
510 info = n + kk - 1 + kp
511 END IF
512*
513* ============================================================
514*
515* Test for the second and third stopping criteria.
516* NOTE: There is no need to test for ABSTOL >= ZERO, since
517* MAXC2NRMK is non-negative. Similarly, there is no need
518* to test for RELTOL >= ZERO, since RELMAXC2NRMK is
519* non-negative.
520* We need to check the condition only if the
521* column index (same as row index) of the original whole
522* matrix is larger than 1, since the condition for whole
523* original matrix is checked in the main routine.
524
525 relmaxc2nrmk = maxc2nrmk / maxc2nrm
526*
527 IF( maxc2nrmk.LE.abstol .OR. relmaxc2nrmk.LE.reltol ) THEN
528*
529* Set K, the number of factorized columns.
530*
531 k = kk - 1
532*
533* Set TAUs corresponding to the columns that were not
534* factorized to ZERO, i.e. set TAU(KK:MINMNFACT) to ZERO.
535*
536 DO j = kk, minmnfact
537 tau( j ) = zero
538 END DO
539*
540* Return from the routine.
541*
542 RETURN
543*
544 END IF
545*
546* ============================================================
547*
548* End ELSE of IF(I.EQ.1)
549*
550 END IF
551*
552* ===============================================================
553*
554* If the pivot column is not the first column of the
555* subblock A(1:M,KK:N):
556* 1) swap the KK-th column and the KP-th pivot column
557* in A(1:M,1:N);
558* 2) copy the KK-th element into the KP-th element of the partial
559* and exact 2-norm vectors VN1 and VN2. ( Swap is not needed
560* for VN1 and VN2 since we use the element with the index
561* larger than KK in the next loop step.)
562* 3) Save the pivot interchange with the indices relative to the
563* the original matrix A, not the block A(1:M,1:N).
564*
565 IF( kp.NE.kk ) THEN
566 CALL sswap( m, a( 1, kp ), 1, a( 1, kk ), 1 )
567 vn1( kp ) = vn1( kk )
568 vn2( kp ) = vn2( kk )
569 itemp = jpiv( kp )
570 jpiv( kp ) = jpiv( kk )
571 jpiv( kk ) = itemp
572 END IF
573*
574* Generate elementary reflector H(KK) using the column A(I:M,KK),
575* if the column has more than one element, otherwise
576* the elementary reflector would be an identity matrix,
577* and TAU(KK) = ZERO.
578*
579 IF( i.LT.m ) THEN
580 CALL slarfg( m-i+1, a( i, kk ), a( i+1, kk ), 1,
581 $ tau( kk ) )
582 ELSE
583 tau( kk ) = zero
584 END IF
585*
586* Check if TAU(KK) contains NaN, set INFO parameter
587* to the column number where NaN is found and return from
588* the routine.
589* NOTE: There is no need to check TAU(KK) for Inf,
590* since SLARFG cannot produce TAU(KK) or Householder vector
591* below the diagonal containing Inf. Only BETA on the diagonal,
592* returned by SLARFG can contain Inf, which requires
593* TAU(KK) to contain NaN. Therefore, this case of generating Inf
594* by SLARFG is covered by checking TAU(KK) for NaN.
595*
596 IF( sisnan( tau(kk) ) ) THEN
597 k = kk - 1
598 info = kk
599*
600* Set MAXC2NRMK and RELMAXC2NRMK to NaN.
601*
602 maxc2nrmk = tau( kk )
603 relmaxc2nrmk = tau( kk )
604*
605* Array TAU(KK:MINMNFACT) is not set and contains
606* undefined elements, except the first element TAU(KK) = NaN.
607*
608 RETURN
609 END IF
610*
611* Apply H(KK)**T to A(I:M,KK+1:N+NRHS) from the left.
612* ( If M >= N, then at KK = N there is no residual matrix,
613* i.e. no columns of A to update, only columns of B.
614* If M < N, then at KK = M-IOFFSET, I = M and we have a
615* one-row residual matrix in A and the elementary
616* reflector is a unit matrix, TAU(KK) = ZERO, i.e. no update
617* is needed for the residual matrix in A and the
618* right-hand-side-matrix in B.
619* Therefore, we update only if
620* KK < MINMNUPDT = min(M-IOFFSET, N+NRHS)
621* condition is satisfied, not only KK < N+NRHS )
622*
623 IF( kk.LT.minmnupdt ) THEN
624 aikk = a( i, kk )
625 a( i, kk ) = one
626 CALL slarf( 'Left', m-i+1, n+nrhs-kk, a( i, kk ), 1,
627 $ tau( kk ), a( i, kk+1 ), lda, work( 1 ) )
628 a( i, kk ) = aikk
629 END IF
630*
631 IF( kk.LT.minmnfact ) THEN
632*
633* Update the partial column 2-norms for the residual matrix,
634* only if the residual matrix A(I+1:M,KK+1:N) exists, i.e.
635* when KK < min(M-IOFFSET, N).
636*
637 DO j = kk + 1, n
638 IF( vn1( j ).NE.zero ) THEN
639*
640* NOTE: The following lines follow from the analysis in
641* Lapack Working Note 176.
642*
643 temp = one - ( abs( a( i, j ) ) / vn1( j ) )**2
644 temp = max( temp, zero )
645 temp2 = temp*( vn1( j ) / vn2( j ) )**2
646 IF( temp2 .LE. tol3z ) THEN
647*
648* Compute the column 2-norm for the partial
649* column A(I+1:M,J) by explicitly computing it,
650* and store it in both partial 2-norm vector VN1
651* and exact column 2-norm vector VN2.
652*
653 vn1( j ) = snrm2( m-i, a( i+1, j ), 1 )
654 vn2( j ) = vn1( j )
655*
656 ELSE
657*
658* Update the column 2-norm for the partial
659* column A(I+1:M,J) by removing one
660* element A(I,J) and store it in partial
661* 2-norm vector VN1.
662*
663 vn1( j ) = vn1( j )*sqrt( temp )
664*
665 END IF
666 END IF
667 END DO
668*
669 END IF
670*
671* End factorization loop
672*
673 END DO
674*
675* If we reached this point, all colunms have been factorized,
676* i.e. no condition was triggered to exit the routine.
677* Set the number of factorized columns.
678*
679 k = kmax
680*
681* We reached the end of the loop, i.e. all KMAX columns were
682* factorized, we need to set MAXC2NRMK and RELMAXC2NRMK before
683* we return.
684*
685 IF( k.LT.minmnfact ) THEN
686*
687 jmaxc2nrm = k + isamax( n-k, vn1( k+1 ), 1 )
688 maxc2nrmk = vn1( jmaxc2nrm )
689*
690 IF( k.EQ.0 ) THEN
691 relmaxc2nrmk = one
692 ELSE
693 relmaxc2nrmk = maxc2nrmk / maxc2nrm
694 END IF
695*
696 ELSE
697 maxc2nrmk = zero
698 relmaxc2nrmk = zero
699 END IF
700*
701* We reached the end of the loop, i.e. all KMAX columns were
702* factorized, set TAUs corresponding to the columns that were
703* not factorized to ZERO, i.e. TAU(K+1:MINMNFACT) set to ZERO.
704*
705 DO j = k + 1, minmnfact
706 tau( j ) = zero
707 END DO
708*
709 RETURN
710*
711* End of SLAQP2RK
712*
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:59
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine slarf(side, m, n, v, incv, tau, c, ldc, work)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition slarf.f:124
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
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: