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

◆ claqp2rk()

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

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

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

Purpose:
 CLAQP2RK computes a truncated (rank K) or full rank Householder QR
 factorization with column pivoting of the complex 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)**H * 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 REAL, 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 REAL, 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 CGEQP3RK. 1 <= KP1 <= N_orig_mat.
[in]MAXC2NRM
          MAXC2NRM is REAL
          The maximum column 2-norm of the whole original
          matrix A_orig computed in the main routine CGEQP3RK.
          MAXC2NRM >= 0.
[in,out]A
          A is COMPLEX 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 unitary 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)**H.
[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 REAL
          The maximum column 2-norm of the residual matrix,
          when the factorization stopped at rank K. MAXC2NRMK >= 0.
[out]RELMAXC2NRMK
          RELMAXC2NRMK is REAL
          The ratio MAXC2NRMK / MAXC2NRM of the maximum column
          2-norm of the residual matrix (when the factorization
          stopped at rank 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 COMPLEX 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 COMPLEX array, dimension (N-1)
          Used in CLARF 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 341 of file claqp2rk.f.

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