1 SUBROUTINE PDGEHDRV( N, ILO, IHI, A, IA, JA, DESCA, TAU, WORK ) 2* 3* -- ScaLAPACK routine (version 1.7) -- 4* University of Tennessee, Knoxville, Oak Ridge National Laboratory, 5* and University of California, Berkeley. 6* May 28, 2001 7* 8* .. Scalar Arguments .. 9 INTEGER IA, IHI, ILO, JA, N 10* .. 11* .. Array Arguments .. 12 INTEGER DESCA( * ) 13 DOUBLE PRECISION A( * ), TAU( * ), WORK( * ) 14* .. 15* 16* Purpose 17* ======= 18* 19* PDGEHDRV computes sub( A ) = A(IA:IA+N-1,JA:JA+N-1) from the 20* orthogonal matrix Q, the Hessenberg matrix, and the array TAU 21* returned by PDGEHRD: 22* sub( A ) := Q * H * Q' 23* 24* Arguments 25* ========= 26* 27* N (global input) INTEGER 28* The number of rows and columns to be operated on, i.e. the 29* order of the distributed submatrix sub( A ). N >= 0. 30* 31* ILO (global input) INTEGER 32* IHI (global input) INTEGER 33* It is assumed that sub( A ) is already upper triangular in 34* rows and columns 1:ILO-1 and IHI+1:N. If N > 0, 35* 1 <= ILO <= IHI <= N; otherwise set ILO = 1, IHI = N. 36* 37* A (local input/local output) DOUBLE PRECISION pointer into the 38* local memory to an array of dimension (LLD_A,LOCc(JA+N-1)). 39* On entry, this array contains the local pieces of the N-by-N 40* general distributed matrix sub( A ) reduced to Hessenberg 41* form by PDGEHRD. The upper triangle and the first sub- 42* diagonal of sub( A ) contain the upper Hessenberg matrix H, 43* and the elements below the first subdiagonal, with the array 44* TAU, represent the orthogonal matrix Q as a product of 45* elementary reflectors. On exit, the original distributed 46* N-by-N matrix sub( A ) is recovered. 47* 48* IA (global input) INTEGER 49* The row index in the global array A indicating the first 50* row of sub( A ). 51* 52* JA (global input) INTEGER 53* The column index in the global array A indicating the 54* first column of sub( A ). 55* 56* DESCA (global and local input) INTEGER array of dimension DLEN_. 57* The array descriptor for the distributed matrix A. 58* 59* TAU (local input) DOUBLE PRECISION array, dimension LOCc(JA+N-2) 60* The scalar factors of the elementary reflectors returned by 61* PDGEHRD. TAU is tied to the distributed matrix A. 62* 63* WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK). 64* LWORK >= NB*NB + NB*IHLP + MAX[ NB*( IHLP+INLQ ), 65* NB*( IHLQ + MAX[ IHIP, 66* IHLP+NUMROC( NUMROC( IHI-ILO+LOFF+1, NB, 0, 0, 67* NPCOL ), NB, 0, 0, LCMQ ) ] ) ] 68* 69* where NB = MB_A = NB_A, 70* LCM is the least common multiple of NPROW and NPCOL, 71* LCM = ILCM( NPROW, NPCOL ), LCMQ = LCM / NPCOL, 72* 73* IROFFA = MOD( IA-1, NB ), 74* IAROW = INDXG2P( IA, NB, MYROW, RSRC_A, NPROW ), 75* IHIP = NUMROC( IHI+IROFFA, NB, MYROW, IAROW, NPROW ), 76* 77* ILROW = INDXG2P( IA+ILO-1, NB, MYROW, RSRC_A, NPROW ), 78* ILCOL = INDXG2P( JA+ILO-1, NB, MYCOL, CSRC_A, NPCOL ), 79* IHLP = NUMROC( IHI-ILO+IROFFA+1, NB, MYROW, ILROW, NPROW ), 80* IHLQ = NUMROC( IHI-ILO+IROFFA+1, NB, MYCOL, ILCOL, NPCOL ), 81* INLQ = NUMROC( N-ILO+IROFFA+1, NB, MYCOL, ILCOL, NPCOL ). 82* 83* ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions; 84* MYROW, MYCOL, NPROW and NPCOL can be determined by calling 85* the subroutine BLACS_GRIDINFO. 86* 87* ===================================================================== 88* 89* .. Parameters .. 90 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 91 $ LLD_, MB_, M_, NB_, N_, RSRC_ 92 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 93 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 94 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 95 DOUBLE PRECISION ZERO 96 PARAMETER ( ZERO = 0.0D+0 ) 97* .. 98* .. Local Scalars .. 99 INTEGER I, IACOL, IAROW, ICTXT, IHLP, II, IOFF, IPT, 100 $ IPV, IPW, IV, J, JB, JJ, JL, K, MYCOL, MYROW, 101 $ NB, NPCOL, NPROW 102* .. 103* .. Local Arrays .. 104 INTEGER DESCV( DLEN_ ) 105* .. 106* .. External Functions .. 107 INTEGER INDXG2P, NUMROC 108 EXTERNAL INDXG2P, NUMROC 109* .. 110* .. External Subroutines .. 111 EXTERNAL BLACS_GRIDINFO, DESCSET, INFOG2L, PDLARFB, 112 $ PDLARFT, PDLACPY, PDLASET 113* .. 114* .. Intrinsic Functions .. 115 INTRINSIC MAX, MIN, MOD 116* .. 117* .. Executable statements .. 118* 119* Get grid parameters 120* 121 ICTXT = DESCA( CTXT_ ) 122 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 123* 124* Quick return if possible 125* 126 IF( IHI-ILO.LE.0 ) 127 $ RETURN 128* 129 NB = DESCA( MB_ ) 130 IOFF = MOD( IA+ILO-2, NB ) 131 CALL INFOG2L( IA+ILO-1, JA+ILO-1, DESCA, NPROW, NPCOL, MYROW, 132 $ MYCOL, II, JJ, IAROW, IACOL ) 133 IHLP = NUMROC( IHI-ILO+IOFF+1, NB, MYROW, IAROW, NPROW ) 134* 135 IPT = 1 136 IPV = IPT + NB * NB 137 IPW = IPV + IHLP * NB 138 JL = MAX( ( ( JA+IHI-2 ) / NB ) * NB + 1, JA + ILO - 1 ) 139 CALL DESCSET( DESCV, IHI-ILO+IOFF+1, NB, NB, NB, IAROW, 140 $ INDXG2P( JL, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ), 141 $ NPCOL ), ICTXT, MAX( 1, IHLP ) ) 142* 143 DO 10 J = JL, ILO+JA+NB-IOFF-1, -NB 144 JB = MIN( JA+IHI-J-1, NB ) 145 I = IA + J - JA 146 K = I - IA + 1 147 IV = K - ILO + IOFF + 1 148* 149* Compute upper triangular matrix T from TAU. 150* 151 CALL PDLARFT( 'Forward', 'Columnwise', IHI-K, JB, A, I+1, J, 152 $ DESCA, TAU, WORK( IPT ), WORK( IPW ) ) 153* 154* Copy Householder vectors into workspace. 155* 156 CALL PDLACPY( 'All', IHI-K, JB, A, I+1, J, DESCA, WORK( IPV ), 157 $ IV+1, 1, DESCV ) 158* 159* Zero out the strict lower triangular part of A. 160* 161 CALL PDLASET( 'Lower', IHI-K-1, JB, ZERO, ZERO, A, I+2, J, 162 $ DESCA ) 163* 164* Apply block Householder transformation from Left. 165* 166 CALL PDLARFB( 'Left', 'No transpose', 'Forward', 'Columnwise', 167 $ IHI-K, N-K+1, JB, WORK( IPV ), IV+1, 1, DESCV, 168 $ WORK( IPT ), A, I+1, J, DESCA, WORK( IPW ) ) 169* 170* Apply block Householder transformation from Right. 171* 172 CALL PDLARFB( 'Right', 'Transpose', 'Forward', 'Columnwise', 173 $ IHI, IHI-K, JB, WORK( IPV ), IV+1, 1, DESCV, 174 $ WORK( IPT ), A, IA, J+1, DESCA, WORK( IPW ) ) 175* 176 DESCV( CSRC_ ) = MOD( DESCV( CSRC_ ) + NPCOL - 1, NPCOL ) 177* 178 10 CONTINUE 179* 180* Handle the first block separately 181* 182 IV = IOFF + 1 183 I = IA + ILO - 1 184 J = JA + ILO - 1 185 JB = MIN( NB-IOFF, JA+IHI-J-1 ) 186* 187* Compute upper triangular matrix T from TAU. 188* 189 CALL PDLARFT( 'Forward', 'Columnwise', IHI-ILO, JB, A, I+1, J, 190 $ DESCA, TAU, WORK( IPT ), WORK( IPW ) ) 191* 192* Copy Householder vectors into workspace. 193* 194 CALL PDLACPY( 'All', IHI-ILO, JB, A, I+1, J, DESCA, WORK( IPV ), 195 $ IV+1, 1, DESCV ) 196* 197* Zero out the strict lower triangular part of A. 198* 199 IF( IHI-ILO.GT.0 ) 200 $ CALL PDLASET( 'Lower', IHI-ILO-1, JB, ZERO, ZERO, A, I+2, J, 201 $ DESCA ) 202* 203* Apply block Householder transformation from Left. 204* 205 CALL PDLARFB( 'Left', 'No transpose', 'Forward', 'Columnwise', 206 $ IHI-ILO, N-ILO+1, JB, WORK( IPV ), IV+1, 1, DESCV, 207 $ WORK( IPT ), A, I+1, J, DESCA, WORK( IPW ) ) 208* 209* Apply block Householder transformation from Right. 210* 211 CALL PDLARFB( 'Right', 'Transpose', 'Forward', 'Columnwise', IHI, 212 $ IHI-ILO, JB, WORK( IPV ), IV+1, 1, DESCV, 213 $ WORK( IPT ), A, IA, J+1, DESCA, WORK( IPW ) ) 214* 215 RETURN 216* 217* End of PDGEHDRV 218* 219 END 220