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

◆ zlaqp3rk()

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

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

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

Purpose:
 ZLAQP3RK computes a step of truncated QR factorization with column
 pivoting of a complex 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)**H * 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 ZGEQP3RK. 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 ZGEQP3RK.
          MAXC2NRM >= 0.
[in,out]A
          A is COMPLEX*16 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 unitary 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)**H.
[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 COMPLEX*16 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 COMPLEX*16 array, dimension (NB)
          Auxiliary vector.
[out]F
          F is COMPLEX*16 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+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 392 of file zlaqp3rk.f.

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