LAPACK 3.3.0
|
00001 SUBROUTINE ZUNGHR( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER IHI, ILO, INFO, LDA, LWORK, N 00010 * .. 00011 * .. Array Arguments .. 00012 COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * ZUNGHR generates a complex unitary matrix Q which is defined as the 00019 * product of IHI-ILO elementary reflectors of order N, as returned by 00020 * ZGEHRD: 00021 * 00022 * Q = H(ilo) H(ilo+1) . . . H(ihi-1). 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * N (input) INTEGER 00028 * The order of the matrix Q. N >= 0. 00029 * 00030 * ILO (input) INTEGER 00031 * IHI (input) INTEGER 00032 * ILO and IHI must have the same values as in the previous call 00033 * of ZGEHRD. Q is equal to the unit matrix except in the 00034 * submatrix Q(ilo+1:ihi,ilo+1:ihi). 00035 * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 00036 * 00037 * A (input/output) COMPLEX*16 array, dimension (LDA,N) 00038 * On entry, the vectors which define the elementary reflectors, 00039 * as returned by ZGEHRD. 00040 * On exit, the N-by-N unitary matrix Q. 00041 * 00042 * LDA (input) INTEGER 00043 * The leading dimension of the array A. LDA >= max(1,N). 00044 * 00045 * TAU (input) COMPLEX*16 array, dimension (N-1) 00046 * TAU(i) must contain the scalar factor of the elementary 00047 * reflector H(i), as returned by ZGEHRD. 00048 * 00049 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00050 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00051 * 00052 * LWORK (input) INTEGER 00053 * The dimension of the array WORK. LWORK >= IHI-ILO. 00054 * For optimum performance LWORK >= (IHI-ILO)*NB, where NB is 00055 * the optimal blocksize. 00056 * 00057 * If LWORK = -1, then a workspace query is assumed; the routine 00058 * only calculates the optimal size of the WORK array, returns 00059 * this value as the first entry of the WORK array, and no error 00060 * message related to LWORK is issued by XERBLA. 00061 * 00062 * INFO (output) INTEGER 00063 * = 0: successful exit 00064 * < 0: if INFO = -i, the i-th argument had an illegal value 00065 * 00066 * ===================================================================== 00067 * 00068 * .. Parameters .. 00069 COMPLEX*16 ZERO, ONE 00070 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 00071 $ ONE = ( 1.0D+0, 0.0D+0 ) ) 00072 * .. 00073 * .. Local Scalars .. 00074 LOGICAL LQUERY 00075 INTEGER I, IINFO, J, LWKOPT, NB, NH 00076 * .. 00077 * .. External Subroutines .. 00078 EXTERNAL XERBLA, ZUNGQR 00079 * .. 00080 * .. External Functions .. 00081 INTEGER ILAENV 00082 EXTERNAL ILAENV 00083 * .. 00084 * .. Intrinsic Functions .. 00085 INTRINSIC MAX, MIN 00086 * .. 00087 * .. Executable Statements .. 00088 * 00089 * Test the input arguments 00090 * 00091 INFO = 0 00092 NH = IHI - ILO 00093 LQUERY = ( LWORK.EQ.-1 ) 00094 IF( N.LT.0 ) THEN 00095 INFO = -1 00096 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN 00097 INFO = -2 00098 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN 00099 INFO = -3 00100 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00101 INFO = -5 00102 ELSE IF( LWORK.LT.MAX( 1, NH ) .AND. .NOT.LQUERY ) THEN 00103 INFO = -8 00104 END IF 00105 * 00106 IF( INFO.EQ.0 ) THEN 00107 NB = ILAENV( 1, 'ZUNGQR', ' ', NH, NH, NH, -1 ) 00108 LWKOPT = MAX( 1, NH )*NB 00109 WORK( 1 ) = LWKOPT 00110 END IF 00111 * 00112 IF( INFO.NE.0 ) THEN 00113 CALL XERBLA( 'ZUNGHR', -INFO ) 00114 RETURN 00115 ELSE IF( LQUERY ) THEN 00116 RETURN 00117 END IF 00118 * 00119 * Quick return if possible 00120 * 00121 IF( N.EQ.0 ) THEN 00122 WORK( 1 ) = 1 00123 RETURN 00124 END IF 00125 * 00126 * Shift the vectors which define the elementary reflectors one 00127 * column to the right, and set the first ilo and the last n-ihi 00128 * rows and columns to those of the unit matrix 00129 * 00130 DO 40 J = IHI, ILO + 1, -1 00131 DO 10 I = 1, J - 1 00132 A( I, J ) = ZERO 00133 10 CONTINUE 00134 DO 20 I = J + 1, IHI 00135 A( I, J ) = A( I, J-1 ) 00136 20 CONTINUE 00137 DO 30 I = IHI + 1, N 00138 A( I, J ) = ZERO 00139 30 CONTINUE 00140 40 CONTINUE 00141 DO 60 J = 1, ILO 00142 DO 50 I = 1, N 00143 A( I, J ) = ZERO 00144 50 CONTINUE 00145 A( J, J ) = ONE 00146 60 CONTINUE 00147 DO 80 J = IHI + 1, N 00148 DO 70 I = 1, N 00149 A( I, J ) = ZERO 00150 70 CONTINUE 00151 A( J, J ) = ONE 00152 80 CONTINUE 00153 * 00154 IF( NH.GT.0 ) THEN 00155 * 00156 * Generate Q(ilo+1:ihi,ilo+1:ihi) 00157 * 00158 CALL ZUNGQR( NH, NH, NH, A( ILO+1, ILO+1 ), LDA, TAU( ILO ), 00159 $ WORK, LWORK, IINFO ) 00160 END IF 00161 WORK( 1 ) = LWKOPT 00162 RETURN 00163 * 00164 * End of ZUNGHR 00165 * 00166 END