LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dtpqrt2.f
Go to the documentation of this file.
00001 *> \brief \b DTPQRT2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DTPQRT2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtpqrt2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtpqrt2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtpqrt2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER   INFO, LDA, LDB, LDT, N, M, L
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), T( LDT, * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *> DTPQRT2 computes a QR factorization of a real "triangular-pentagonal"
00037 *> matrix C, which is composed of a triangular block A and pentagonal block B, 
00038 *> using the compact WY representation for Q.
00039 *> \endverbatim
00040 *
00041 *  Arguments:
00042 *  ==========
00043 *
00044 *> \param[in] M
00045 *> \verbatim
00046 *>          M is INTEGER
00047 *>          The total number of rows of the matrix B.  
00048 *>          M >= 0.
00049 *> \endverbatim
00050 *>
00051 *> \param[in] N
00052 *> \verbatim
00053 *>          N is INTEGER
00054 *>          The number of columns of the matrix B, and the order of
00055 *>          the triangular matrix A.
00056 *>          N >= 0.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] L
00060 *> \verbatim
00061 *>          L is INTEGER
00062 *>          The number of rows of the upper trapezoidal part of B.  
00063 *>          MIN(M,N) >= L >= 0.  See Further Details.
00064 *> \endverbatim
00065 *>
00066 *> \param[in,out] A
00067 *> \verbatim
00068 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00069 *>          On entry, the upper triangular N-by-N matrix A.
00070 *>          On exit, the elements on and above the diagonal of the array
00071 *>          contain the upper triangular matrix R.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] LDA
00075 *> \verbatim
00076 *>          LDA is INTEGER
00077 *>          The leading dimension of the array A.  LDA >= max(1,N).
00078 *> \endverbatim
00079 *>
00080 *> \param[in,out] B
00081 *> \verbatim
00082 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
00083 *>          On entry, the pentagonal M-by-N matrix B.  The first M-L rows 
00084 *>          are rectangular, and the last L rows are upper trapezoidal.
00085 *>          On exit, B contains the pentagonal matrix V.  See Further Details.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] LDB
00089 *> \verbatim
00090 *>          LDB is INTEGER
00091 *>          The leading dimension of the array B.  LDB >= max(1,M).
00092 *> \endverbatim
00093 *>
00094 *> \param[out] T
00095 *> \verbatim
00096 *>          T is DOUBLE PRECISION array, dimension (LDT,N)
00097 *>          The N-by-N upper triangular factor T of the block reflector.
00098 *>          See Further Details.
00099 *> \endverbatim
00100 *>
00101 *> \param[in] LDT
00102 *> \verbatim
00103 *>          LDT is INTEGER
00104 *>          The leading dimension of the array T.  LDT >= max(1,N)
00105 *> \endverbatim
00106 *>
00107 *> \param[out] INFO
00108 *> \verbatim
00109 *>          INFO is INTEGER
00110 *>          = 0: successful exit
00111 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00112 *> \endverbatim
00113 *
00114 *  Authors:
00115 *  ========
00116 *
00117 *> \author Univ. of Tennessee 
00118 *> \author Univ. of California Berkeley 
00119 *> \author Univ. of Colorado Denver 
00120 *> \author NAG Ltd. 
00121 *
00122 *> \date April 2012
00123 *
00124 *> \ingroup doubleOTHERcomputational
00125 *
00126 *> \par Further Details:
00127 *  =====================
00128 *>
00129 *> \verbatim
00130 *>
00131 *>  The input matrix C is a (N+M)-by-N matrix  
00132 *>
00133 *>               C = [ A ]
00134 *>                   [ B ]        
00135 *>
00136 *>  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
00137 *>  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
00138 *>  upper trapezoidal matrix B2:
00139 *>
00140 *>               B = [ B1 ]  <- (M-L)-by-N rectangular
00141 *>                   [ B2 ]  <-     L-by-N upper trapezoidal.
00142 *>
00143 *>  The upper trapezoidal matrix B2 consists of the first L rows of a
00144 *>  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0, 
00145 *>  B is rectangular M-by-N; if M=L=N, B is upper triangular.  
00146 *>
00147 *>  The matrix W stores the elementary reflectors H(i) in the i-th column
00148 *>  below the diagonal (of A) in the (N+M)-by-N input matrix C
00149 *>
00150 *>               C = [ A ]  <- upper triangular N-by-N
00151 *>                   [ B ]  <- M-by-N pentagonal
00152 *>
00153 *>  so that W can be represented as
00154 *>
00155 *>               W = [ I ]  <- identity, N-by-N
00156 *>                   [ V ]  <- M-by-N, same form as B.
00157 *>
00158 *>  Thus, all of information needed for W is contained on exit in B, which
00159 *>  we call V above.  Note that V has the same form as B; that is, 
00160 *>
00161 *>               V = [ V1 ] <- (M-L)-by-N rectangular
00162 *>                   [ V2 ] <-     L-by-N upper trapezoidal.
00163 *>
00164 *>  The columns of V represent the vectors which define the H(i)'s.  
00165 *>  The (M+N)-by-(M+N) block reflector H is then given by
00166 *>
00167 *>               H = I - W * T * W**T
00168 *>
00169 *>  where W^H is the conjugate transpose of W and T is the upper triangular
00170 *>  factor of the block reflector.
00171 *> \endverbatim
00172 *>
00173 *  =====================================================================
00174       SUBROUTINE DTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
00175 *
00176 *  -- LAPACK computational routine (version 3.4.1) --
00177 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00178 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00179 *     April 2012
00180 *
00181 *     .. Scalar Arguments ..
00182       INTEGER   INFO, LDA, LDB, LDT, N, M, L
00183 *     ..
00184 *     .. Array Arguments ..
00185       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), T( LDT, * )
00186 *     ..
00187 *
00188 *  =====================================================================
00189 *
00190 *     .. Parameters ..
00191       DOUBLE PRECISION  ONE, ZERO
00192       PARAMETER( ONE = 1.0, ZERO = 0.0 )
00193 *     ..
00194 *     .. Local Scalars ..
00195       INTEGER   I, J, P, MP, NP
00196       DOUBLE PRECISION   ALPHA
00197 *     ..
00198 *     .. External Subroutines ..
00199       EXTERNAL  DLARFG, DGEMV, DGER, DTRMV, XERBLA
00200 *     ..
00201 *     .. Intrinsic Functions ..
00202       INTRINSIC MAX, MIN
00203 *     ..
00204 *     .. Executable Statements ..
00205 *
00206 *     Test the input arguments
00207 *
00208       INFO = 0
00209       IF( M.LT.0 ) THEN
00210          INFO = -1
00211       ELSE IF( N.LT.0 ) THEN
00212          INFO = -2
00213       ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
00214          INFO = -3
00215       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00216          INFO = -5
00217       ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
00218          INFO = -7
00219       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00220          INFO = -9
00221       END IF
00222       IF( INFO.NE.0 ) THEN
00223          CALL XERBLA( 'DTPQRT2', -INFO )
00224          RETURN
00225       END IF
00226 *
00227 *     Quick return if possible
00228 *
00229       IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
00230 *      
00231       DO I = 1, N
00232 *
00233 *        Generate elementary reflector H(I) to annihilate B(:,I)
00234 *
00235          P = M-L+MIN( L, I )
00236          CALL DLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
00237          IF( I.LT.N ) THEN
00238 *
00239 *           W(1:N-I) := C(I:M,I+1:N)^H * C(I:M,I) [use W = T(:,N)]
00240 *
00241             DO J = 1, N-I
00242                T( J, N ) = (A( I, I+J ))
00243             END DO
00244             CALL DGEMV( 'T', P, N-I, ONE, B( 1, I+1 ), LDB, 
00245      $                  B( 1, I ), 1, ONE, T( 1, N ), 1 )
00246 *
00247 *           C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)^H
00248 *
00249             ALPHA = -(T( I, 1 ))            
00250             DO J = 1, N-I
00251                A( I, I+J ) = A( I, I+J ) + ALPHA*(T( J, N ))
00252             END DO
00253             CALL DGER( P, N-I, ALPHA, B( 1, I ), 1, 
00254      $           T( 1, N ), 1, B( 1, I+1 ), LDB )
00255          END IF
00256       END DO
00257 *
00258       DO I = 2, N
00259 *
00260 *        T(1:I-1,I) := C(I:M,1:I-1)^H * (alpha * C(I:M,I))
00261 *
00262          ALPHA = -T( I, 1 )
00263 
00264          DO J = 1, I-1
00265             T( J, I ) = ZERO
00266          END DO
00267          P = MIN( I-1, L )
00268          MP = MIN( M-L+1, M )
00269          NP = MIN( P+1, N )
00270 *
00271 *        Triangular part of B2
00272 *
00273          DO J = 1, P
00274             T( J, I ) = ALPHA*B( M-L+J, I )
00275          END DO
00276          CALL DTRMV( 'U', 'T', 'N', P, B( MP, 1 ), LDB,
00277      $               T( 1, I ), 1 )
00278 *
00279 *        Rectangular part of B2
00280 *
00281          CALL DGEMV( 'T', L, I-1-P, ALPHA, B( MP, NP ), LDB, 
00282      $               B( MP, I ), 1, ZERO, T( NP, I ), 1 )
00283 *
00284 *        B1
00285 *
00286          CALL DGEMV( 'T', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, 
00287      $               ONE, T( 1, I ), 1 )         
00288 *
00289 *        T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
00290 *
00291          CALL DTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
00292 *
00293 *        T(I,I) = tau(I)
00294 *
00295          T( I, I ) = T( I, 1 )
00296          T( I, 1 ) = ZERO
00297       END DO
00298    
00299 *
00300 *     End of DTPQRT2
00301 *
00302       END
 All Files Functions