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

◆ dlaqp3rk()

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

DLAQP3RK 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 DLAQP3RK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLAQP3RK 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 DOUBLE PRECISION, 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 DOUBLE PRECISION, 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 DGEQP3RK. 1 <= KP1 <= N_orig.
[in]MAXC2NRM
          MAXC2NRM is DOUBLE PRECISION
          The maximum column 2-norm of the whole original
          matrix A_orig computed in the main routine DGEQP3RK.
          MAXC2NRM >= 0.
[in,out]A
          A is DOUBLE PRECISION 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]

verbatim 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.

Parameters
[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 DOUBLE PRECISION
          The maximum column 2-norm of the residual matrix,
          when the factorization stopped at rank KB. 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 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 DOUBLE PRECISION array, dimension (min(M-IOFFSET,N))
          The scalar factors of the elementary reflectors.
[in,out]VN1
          VN1 is DOUBLE PRECISION array, dimension (N)
          The vector with the partial column norms.
[in,out]VN2
          VN2 is DOUBLE PRECISION array, dimension (N)
          The vector with the exact column norms.
[out]AUXV
          AUXV is DOUBLE PRECISION array, dimension (NB)
          Auxiliary vector.
[out]F
          F is DOUBLE PRECISION 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 "bad" 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 398 of file dlaqp3rk.f.

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