LAPACK 3.3.1
Linear Algebra PACKage

claqp2.f

Go to the documentation of this file.
00001       SUBROUTINE CLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
00002      $                   WORK )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.3.1) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *  -- April 2011                                                      --
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            LDA, M, N, OFFSET
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            JPVT( * )
00014       REAL               VN1( * ), VN2( * )
00015       COMPLEX            A( LDA, * ), TAU( * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CLAQP2 computes a QR factorization with column pivoting of
00022 *  the block A(OFFSET+1:M,1:N).
00023 *  The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
00024 *
00025 *  Arguments
00026 *  =========
00027 *
00028 *  M       (input) INTEGER
00029 *          The number of rows of the matrix A. M >= 0.
00030 *
00031 *  N       (input) INTEGER
00032 *          The number of columns of the matrix A. N >= 0.
00033 *
00034 *  OFFSET  (input) INTEGER
00035 *          The number of rows of the matrix A that must be pivoted
00036 *          but no factorized. OFFSET >= 0.
00037 *
00038 *  A       (input/output) COMPLEX array, dimension (LDA,N)
00039 *          On entry, the M-by-N matrix A.
00040 *          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is 
00041 *          the triangular factor obtained; the elements in block
00042 *          A(OFFSET+1:M,1:N) below the diagonal, together with the
00043 *          array TAU, represent the orthogonal matrix Q as a product of
00044 *          elementary reflectors. Block A(1:OFFSET,1:N) has been
00045 *          accordingly pivoted, but no factorized.
00046 *
00047 *  LDA     (input) INTEGER
00048 *          The leading dimension of the array A. LDA >= max(1,M).
00049 *
00050 *  JPVT    (input/output) INTEGER array, dimension (N)
00051 *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
00052 *          to the front of A*P (a leading column); if JPVT(i) = 0,
00053 *          the i-th column of A is a free column.
00054 *          On exit, if JPVT(i) = k, then the i-th column of A*P
00055 *          was the k-th column of A.
00056 *
00057 *  TAU     (output) COMPLEX array, dimension (min(M,N))
00058 *          The scalar factors of the elementary reflectors.
00059 *
00060 *  VN1     (input/output) REAL array, dimension (N)
00061 *          The vector with the partial column norms.
00062 *
00063 *  VN2     (input/output) REAL array, dimension (N)
00064 *          The vector with the exact column norms.
00065 *
00066 *  WORK    (workspace) COMPLEX array, dimension (N)
00067 *
00068 *  Further Details
00069 *  ===============
00070 *
00071 *  Based on contributions by
00072 *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
00073 *    X. Sun, Computer Science Dept., Duke University, USA
00074 *
00075 *  Partial column norm updating strategy modified by
00076 *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
00077 *    University of Zagreb, Croatia.
00078 *  -- April 2011                                                      --
00079 *  For more details see LAPACK Working Note 176.
00080 *  =====================================================================
00081 *
00082 *     .. Parameters ..
00083       REAL               ZERO, ONE
00084       COMPLEX            CONE
00085       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0,
00086      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00087 *     ..
00088 *     .. Local Scalars ..
00089       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
00090       REAL               TEMP, TEMP2, TOL3Z
00091       COMPLEX            AII
00092 *     ..
00093 *     .. External Subroutines ..
00094       EXTERNAL           CLARF, CLARFG, CSWAP
00095 *     ..
00096 *     .. Intrinsic Functions ..
00097       INTRINSIC          ABS, CONJG, MAX, MIN, SQRT
00098 *     ..
00099 *     .. External Functions ..
00100       INTEGER            ISAMAX
00101       REAL               SCNRM2, SLAMCH
00102       EXTERNAL           ISAMAX, SCNRM2, SLAMCH
00103 *     ..
00104 *     .. Executable Statements ..
00105 *
00106       MN = MIN( M-OFFSET, N )
00107       TOL3Z = SQRT(SLAMCH('Epsilon'))
00108 *
00109 *     Compute factorization.
00110 *
00111       DO 20 I = 1, MN
00112 *
00113          OFFPI = OFFSET + I
00114 *
00115 *        Determine ith pivot column and swap if necessary.
00116 *
00117          PVT = ( I-1 ) + ISAMAX( N-I+1, VN1( I ), 1 )
00118 *
00119          IF( PVT.NE.I ) THEN
00120             CALL CSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
00121             ITEMP = JPVT( PVT )
00122             JPVT( PVT ) = JPVT( I )
00123             JPVT( I ) = ITEMP
00124             VN1( PVT ) = VN1( I )
00125             VN2( PVT ) = VN2( I )
00126          END IF
00127 *
00128 *        Generate elementary reflector H(i).
00129 *
00130          IF( OFFPI.LT.M ) THEN
00131             CALL CLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
00132      $                   TAU( I ) )
00133          ELSE
00134             CALL CLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
00135          END IF
00136 *
00137          IF( I.LT.N ) THEN
00138 *
00139 *           Apply H(i)**H to A(offset+i:m,i+1:n) from the left.
00140 *
00141             AII = A( OFFPI, I )
00142             A( OFFPI, I ) = CONE
00143             CALL CLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
00144      $                  CONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA,
00145      $                  WORK( 1 ) )
00146             A( OFFPI, I ) = AII
00147          END IF
00148 *
00149 *        Update partial column norms.
00150 *
00151          DO 10 J = I + 1, N
00152             IF( VN1( J ).NE.ZERO ) THEN
00153 *
00154 *              NOTE: The following 4 lines follow from the analysis in
00155 *              Lapack Working Note 176.
00156 *
00157                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
00158                TEMP = MAX( TEMP, ZERO )
00159                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
00160                IF( TEMP2 .LE. TOL3Z ) THEN
00161                   IF( OFFPI.LT.M ) THEN
00162                      VN1( J ) = SCNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
00163                      VN2( J ) = VN1( J )
00164                   ELSE
00165                      VN1( J ) = ZERO
00166                      VN2( J ) = ZERO
00167                   END IF
00168                ELSE
00169                   VN1( J ) = VN1( J )*SQRT( TEMP )
00170                END IF
00171             END IF
00172    10    CONTINUE
00173 *
00174    20 CONTINUE
00175 *
00176       RETURN
00177 *
00178 *     End of CLAQP2
00179 *
00180       END
 All Files Functions