LAPACK 3.12.1
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]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 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  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 382 of file zlaqp3rk.f.

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