|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE PZGEBDRV( M, N, A, IA, JA, DESCA, D, E, TAUQ, TAUP, 00002 $ WORK, INFO ) 00003 * 00004 * -- ScaLAPACK routine (version 1.7) -- 00005 * University of Tennessee, Knoxville, Oak Ridge National Laboratory, 00006 * and University of California, Berkeley. 00007 * May 1, 1997 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, IA, JA, M, N 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER DESCA( * ) 00014 DOUBLE PRECISION D( * ), E( * ) 00015 COMPLEX*16 A( * ), TAUP( * ), TAUQ( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * PZGEBDRV computes sub( A ) = A(IA:IA+M-1,JA:JA+N-1) from sub( A ), 00022 * Q, P returned by PZGEBRD: 00023 * 00024 * sub( A ) := Q * B * P'. 00025 * 00026 * Notes 00027 * ===== 00028 * 00029 * Each global data object is described by an associated description 00030 * vector. This vector stores the information required to establish 00031 * the mapping between an object element and its corresponding process 00032 * and memory location. 00033 * 00034 * Let A be a generic term for any 2D block cyclicly distributed array. 00035 * Such a global array has an associated description vector DESCA. 00036 * In the following comments, the character _ should be read as 00037 * "of the global array". 00038 * 00039 * NOTATION STORED IN EXPLANATION 00040 * --------------- -------------- -------------------------------------- 00041 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 00042 * DTYPE_A = 1. 00043 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 00044 * the BLACS process grid A is distribu- 00045 * ted over. The context itself is glo- 00046 * bal, but the handle (the integer 00047 * value) may vary. 00048 * M_A (global) DESCA( M_ ) The number of rows in the global 00049 * array A. 00050 * N_A (global) DESCA( N_ ) The number of columns in the global 00051 * array A. 00052 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 00053 * the rows of the array. 00054 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 00055 * the columns of the array. 00056 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 00057 * row of the array A is distributed. 00058 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the 00059 * first column of the array A is 00060 * distributed. 00061 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 00062 * array. LLD_A >= MAX(1,LOCr(M_A)). 00063 * 00064 * Let K be the number of rows or columns of a distributed matrix, 00065 * and assume that its process grid has dimension p x q. 00066 * LOCr( K ) denotes the number of elements of K that a process 00067 * would receive if K were distributed over the p processes of its 00068 * process column. 00069 * Similarly, LOCc( K ) denotes the number of elements of K that a 00070 * process would receive if K were distributed over the q processes of 00071 * its process row. 00072 * The values of LOCr() and LOCc() may be determined via a call to the 00073 * ScaLAPACK tool function, NUMROC: 00074 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 00075 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 00076 * An upper bound for these quantities may be computed by: 00077 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 00078 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 00079 * 00080 * Arguments 00081 * ========= 00082 * 00083 * M (global input) INTEGER 00084 * The number of rows to be operated on, i.e. the number of rows 00085 * of the distributed submatrix sub( A ). M >= 0. 00086 * 00087 * N (global input) INTEGER 00088 * The number of columns to be operated on, i.e. the number of 00089 * columns of the distributed submatrix sub( A ). N >= 0. 00090 * 00091 * A (local input/local output) COMPLEX*16 pointer into the 00092 * local memory to an array of dimension (LLD_A, LOCc(JA+N-1)). 00093 * On entry, this array contains the local pieces of sub( A ) 00094 * as returned by PZGEBRD. On exit, the original distribu- 00095 * ted matrix sub( A ) is restored. 00096 * 00097 * IA (global input) INTEGER 00098 * The row index in the global array A indicating the first 00099 * row of sub( A ). 00100 * 00101 * JA (global input) INTEGER 00102 * The column index in the global array A indicating the 00103 * first column of sub( A ). 00104 * 00105 * DESCA (global and local input) INTEGER array of dimension DLEN_. 00106 * The array descriptor for the distributed matrix A. 00107 * 00108 * D (local input) DOUBLE PRECISION array, dimension 00109 * LOCc(JA+MIN(M,N)-1) if M >= N; LOCr(IA+MIN(M,N)-1) otherwise. 00110 * The distributed diagonal elements of the bidiagonal matrix 00111 * B: D(i) = A(i,i). D is tied to the distributed matrix A. 00112 * 00113 * E (local input) DOUBLE PRECISION array, dimension 00114 * LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-2) otherwise. 00115 * The distributed off-diagonal elements of the bidiagonal 00116 * distributed matrix B: 00117 * if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1; 00118 * if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1. 00119 * E is tied to the distributed matrix A. 00120 * 00121 * TAUQ (local input) COMPLEX*16 array dimension 00122 * LOCc(JA+MIN(M,N)-1). The scalar factors of the elementary 00123 * reflectors which represent the unitary matrix Q. TAUQ is 00124 * tied to the distributed matrix A. See Further Details. 00125 * 00126 * TAUP (local input) COMPLEX*16 array, dimension 00127 * LOCr(IA+MIN(M,N)-1). The scalar factors of the elementary 00128 * reflectors which represent the unitary matrix P. TAUP is 00129 * tied to the distributed matrix A. See Further Details. 00130 * 00131 * WORK (local workspace) COMPLEX*16 array, dimension (LWORK) 00132 * LWORK >= 2*NB*( MP + NQ + NB ) 00133 * 00134 * where NB = MB_A = NB_A, 00135 * IROFFA = MOD( IA-1, NB ), ICOFFA = MOD( JA-1, NB ), 00136 * IAROW = INDXG2P( IA, NB, MYROW, RSRC_A, NPROW ), 00137 * IACOL = INDXG2P( JA, NB, MYCOL, CSRC_A, NPCOL ), 00138 * MP = NUMROC( M+IROFFA, NB, MYROW, IAROW, NPROW ), 00139 * NQ = NUMROC( N+ICOFFA, NB, MYCOL, IACOL, NPCOL ). 00140 * 00141 * INDXG2P and NUMROC are ScaLAPACK tool functions; 00142 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling 00143 * the subroutine BLACS_GRIDINFO. 00144 * 00145 * INFO (global output) INTEGER 00146 * On exit, if INFO <> 0, a discrepancy has been found between 00147 * the diagonal and off-diagonal elements of A and the copies 00148 * contained in the arrays D and E. 00149 * 00150 * ===================================================================== 00151 * 00152 * .. Parameters .. 00153 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 00154 $ LLD_, MB_, M_, NB_, N_, RSRC_ 00155 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00156 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00157 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00158 DOUBLE PRECISION REIGHT, RZERO 00159 PARAMETER ( REIGHT = 8.0D+0, RZERO = 0.0D+0 ) 00160 COMPLEX*16 ONE, ZERO 00161 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 00162 $ ZERO = ( 0.0D+0, 0.0D+0 ) ) 00163 * .. 00164 * .. Local Scalars .. 00165 INTEGER I, IACOL, IAROW, ICTXT, IIA, IL, IPTP, IPTQ, 00166 $ IPV, IPW, IPWK, IOFF, IV, J, JB, JJA, JL, JV, 00167 $ K, MN, MP, MYCOL, MYROW, NB, NPCOL, NPROW, NQ 00168 DOUBLE PRECISION ADDBND, D2, E2 00169 COMPLEX*16 D1, E1 00170 * .. 00171 * .. Local Arrays .. 00172 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ), DESCV( DLEN_ ), 00173 $ DESCW( DLEN_ ) 00174 * .. 00175 * .. External Subroutines .. 00176 EXTERNAL BLACS_GRIDINFO, DESCSET, IGSUM2D, INFOG2L, 00177 $ PDELGET, PZLACPY, PZLARFB, PZLARFT, 00178 $ PZLASET, PZELGET 00179 * .. 00180 * .. External Functions .. 00181 INTEGER INDXG2P, NUMROC 00182 DOUBLE PRECISION PDLAMCH 00183 EXTERNAL INDXG2P, NUMROC, PDLAMCH 00184 * .. 00185 * .. Intrinsic Functions .. 00186 INTRINSIC ABS, DCMPLX, MAX, MIN, MOD 00187 * .. 00188 * .. Executable Statements .. 00189 * 00190 * Get grid parameters 00191 * 00192 ICTXT = DESCA( CTXT_ ) 00193 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 00194 * 00195 INFO = 0 00196 NB = DESCA( MB_ ) 00197 IOFF = MOD( IA-1, NB ) 00198 CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, IIA, JJA, 00199 $ IAROW, IACOL ) 00200 MP = NUMROC( M+IOFF, NB, MYROW, IAROW, NPROW ) 00201 NQ = NUMROC( N+IOFF, NB, MYCOL, IACOL, NPCOL ) 00202 IPV = 1 00203 IPW = IPV + MP*NB 00204 IPTP = IPW + NQ*NB 00205 IPTQ = IPTP + NB*NB 00206 IPWK = IPTQ + NB*NB 00207 * 00208 IV = 1 00209 JV = 1 00210 MN = MIN( M, N ) 00211 IL = MAX( ( (IA+MN-2) / NB )*NB + 1, IA ) 00212 JL = MAX( ( (JA+MN-2) / NB )*NB + 1, JA ) 00213 IAROW = INDXG2P( IL, NB, MYROW, DESCA( RSRC_ ), NPROW ) 00214 IACOL = INDXG2P( JL, NB, MYCOL, DESCA( CSRC_ ), NPCOL ) 00215 CALL DESCSET( DESCV, IA+M-IL, NB, NB, NB, IAROW, IACOL, ICTXT, 00216 $ MAX( 1, MP ) ) 00217 CALL DESCSET( DESCW, NB, JA+N-JL, NB, NB, IAROW, IACOL, ICTXT, 00218 $ NB ) 00219 * 00220 ADDBND = REIGHT * PDLAMCH( ICTXT, 'eps' ) 00221 * 00222 * When A is an upper bidiagonal form 00223 * 00224 IF( M.GE.N ) THEN 00225 * 00226 CALL DESCSET( DESCD, 1, JA+MN-1, 1, DESCA( NB_ ), MYROW, 00227 $ DESCA( CSRC_ ), DESCA( CTXT_ ), 1 ) 00228 CALL DESCSET( DESCE, IA+MN-1, 1, DESCA( MB_ ), 1, 00229 $ DESCA( RSRC_ ), MYCOL, DESCA( CTXT_ ), 00230 $ DESCA( LLD_ ) ) 00231 * 00232 DO 10 J = 0, MN-1 00233 D1 = ZERO 00234 E1 = ZERO 00235 D2 = RZERO 00236 E2 = RZERO 00237 CALL PDELGET( ' ', ' ', D2, D, 1, JA+J, DESCD ) 00238 CALL PZELGET( 'Columnwise', ' ', D1, A, IA+J, JA+J, DESCA ) 00239 IF( J.LT.(MN-1) ) THEN 00240 CALL PDELGET( ' ', ' ', E2, E, IA+J, 1, DESCE ) 00241 CALL PZELGET( 'Rowwise', ' ', E1, A, IA+J, JA+J+1, 00242 $ DESCA ) 00243 END IF 00244 * 00245 IF( ( ABS( D1-DCMPLX( D2 ) ).GT.( ABS( D2 )*ADDBND ) ) .OR. 00246 $ ( ABS( E1-DCMPLX( E2 ) ).GT.( ABS( E2 )*ADDBND ) ) ) 00247 $ INFO = INFO + 1 00248 10 CONTINUE 00249 * 00250 DO 20 J = JL, JA+NB-IOFF, -NB 00251 JB = MIN( JA+N-J, NB ) 00252 I = IA + J - JA 00253 K = I - IA + 1 00254 * 00255 * Compute upper triangular matrix TQ from TAUQ. 00256 * 00257 CALL PZLARFT( 'Forward', 'Columnwise', M-K+1, JB, A, I, J, 00258 $ DESCA, TAUQ, WORK( IPTQ ), WORK( IPWK ) ) 00259 * 00260 * Copy Householder vectors into workspace. 00261 * 00262 CALL PZLACPY( 'Lower', M-K+1, JB, A, I, J, DESCA, 00263 $ WORK( IPV ), IV, JV, DESCV ) 00264 CALL PZLASET( 'Upper', M-K+1, JB, ZERO, ONE, WORK( IPV ), 00265 $ IV, JV, DESCV ) 00266 * 00267 * Zero out the strict lower triangular part of A. 00268 * 00269 CALL PZLASET( 'Lower', M-K, JB, ZERO, ZERO, A, I+1, J, 00270 $ DESCA ) 00271 * 00272 * Compute upper triangular matrix TP from TAUP. 00273 * 00274 CALL PZLARFT( 'Forward', 'Rowwise', N-K, JB, A, I, J+1, 00275 $ DESCA, TAUP, WORK( IPTP ), WORK( IPWK ) ) 00276 * 00277 * Copy Householder vectors into workspace. 00278 * 00279 CALL PZLACPY( 'Upper', JB, N-K, A, I, J+1, DESCA, 00280 $ WORK( IPW ), IV, JV+1, DESCW ) 00281 CALL PZLASET( 'Lower', JB, N-K, ZERO, ONE, WORK( IPW ), IV, 00282 $ JV+1, DESCW ) 00283 * 00284 * Zero out the strict+1 upper triangular part of A. 00285 * 00286 CALL PZLASET( 'Upper', JB, N-K-1, ZERO, ZERO, A, I, J+2, 00287 $ DESCA ) 00288 * 00289 * Apply block Householder transformation from Left. 00290 * 00291 CALL PZLARFB( 'Left', 'No transpose', 'Forward', 00292 $ 'Columnwise', M-K+1, N-K+1, JB, WORK( IPV ), 00293 $ IV, JV, DESCV, WORK( IPTQ ), A, I, J, DESCA, 00294 $ WORK( IPWK ) ) 00295 * 00296 * Apply block Householder transformation from Right. 00297 * 00298 CALL PZLARFB( 'Right', 'Conjugate transpose', 'Forward', 00299 $ 'Rowwise', M-K+1, N-K, JB, WORK( IPW ), IV, 00300 $ JV+1, DESCW, WORK( IPTP ), A, I, J+1, DESCA, 00301 $ WORK( IPWK ) ) 00302 * 00303 DESCV( M_ ) = DESCV( M_ ) + NB 00304 DESCV( RSRC_ ) = MOD( DESCV( RSRC_ ) + NPROW - 1, NPROW ) 00305 DESCV( CSRC_ ) = MOD( DESCV( CSRC_ ) + NPCOL - 1, NPCOL ) 00306 DESCW( N_ ) = DESCW( N_ ) + NB 00307 DESCW( RSRC_ ) = DESCV( RSRC_ ) 00308 DESCW( CSRC_ ) = DESCV( CSRC_ ) 00309 * 00310 20 CONTINUE 00311 * 00312 * Handle first block separately 00313 * 00314 JB = MIN( N, NB - IOFF ) 00315 IV = IOFF + 1 00316 JV = IOFF + 1 00317 * 00318 * Compute upper triangular matrix TQ from TAUQ. 00319 * 00320 CALL PZLARFT( 'Forward', 'Columnwise', M, JB, A, IA, JA, DESCA, 00321 $ TAUQ, WORK( IPTQ ), WORK( IPWK ) ) 00322 * 00323 * Copy Householder vectors into workspace. 00324 * 00325 CALL PZLACPY( 'Lower', M, JB, A, IA, JA, DESCA, WORK( IPV ), 00326 $ IV, JV, DESCV ) 00327 CALL PZLASET( 'Upper', M, JB, ZERO, ONE, WORK( IPV ), IV, JV, 00328 $ DESCV ) 00329 * 00330 * Zero out the strict lower triangular part of A. 00331 * 00332 CALL PZLASET( 'Lower', M-1, JB, ZERO, ZERO, A, IA+1, JA, 00333 $ DESCA ) 00334 * 00335 * Compute upper triangular matrix TP from TAUP. 00336 * 00337 CALL PZLARFT( 'Forward', 'Rowwise', N-1, JB, A, IA, JA+1, 00338 $ DESCA, TAUP, WORK( IPTP ), WORK( IPWK ) ) 00339 * 00340 * Copy Householder vectors into workspace. 00341 * 00342 CALL PZLACPY( 'Upper', JB, N-1, A, IA, JA+1, DESCA, 00343 $ WORK( IPW ), IV, JV+1, DESCW ) 00344 CALL PZLASET( 'Lower', JB, N-1, ZERO, ONE, WORK( IPW ), IV, 00345 $ JV+1, DESCW ) 00346 * 00347 * Zero out the strict+1 upper triangular part of A. 00348 * 00349 CALL PZLASET( 'Upper', JB, N-2, ZERO, ZERO, A, IA, JA+2, 00350 $ DESCA ) 00351 * 00352 * Apply block Householder transformation from left. 00353 * 00354 CALL PZLARFB( 'Left', 'No transpose', 'Forward', 'Columnwise', 00355 $ M, N, JB, WORK( IPV ), IV, JV, DESCV, 00356 $ WORK( IPTQ ), A, IA, JA, DESCA, WORK( IPWK ) ) 00357 * 00358 * Apply block Householder transformation from right. 00359 * 00360 CALL PZLARFB( 'Right', 'Conjugate transpose', 'Forward', 00361 $ 'Rowwise', M, N-1, JB, WORK( IPW ), IV, JV+1, 00362 $ DESCW, WORK( IPTP ), A, IA, JA+1, DESCA, 00363 $ WORK( IPWK ) ) 00364 * 00365 ELSE 00366 * 00367 CALL DESCSET( DESCD, IA+MN-1, 1, DESCA( MB_ ), 1, 00368 $ DESCA( RSRC_ ), MYCOL, DESCA( CTXT_ ), 00369 $ DESCA( LLD_ ) ) 00370 CALL DESCSET( DESCE, 1, JA+MN-2, 1, DESCA( NB_ ), MYROW, 00371 $ DESCA( CSRC_ ), DESCA( CTXT_ ), 1 ) 00372 * 00373 DO 30 J = 0, MN-1 00374 D1 = ZERO 00375 E1 = ZERO 00376 D2 = RZERO 00377 E2 = RZERO 00378 CALL PDELGET( ' ', ' ', D2, D, IA+J, 1, DESCD ) 00379 CALL PZELGET( 'Rowwise', ' ', D1, A, IA+J, JA+J, DESCA ) 00380 IF( J.LT.(MN-1) ) THEN 00381 CALL PDELGET( ' ', ' ', E2, E, 1, JA+J, DESCE ) 00382 CALL PZELGET( 'Columnwise', ' ', E1, A, IA+J+1, JA+J, 00383 $ DESCA ) 00384 END IF 00385 * 00386 IF( ( ABS( D1-DCMPLX( D2 ) ).GT.( ABS( D2 )*ADDBND ) ) .OR. 00387 $ ( ABS( E1-DCMPLX( E2 ) ).GT.( ABS( E2 )*ADDBND ) ) ) 00388 $ INFO = INFO + 1 00389 30 CONTINUE 00390 * 00391 DO 40 I = IL, IA+NB-IOFF, -NB 00392 JB = MIN( IA+M-I, NB ) 00393 J = JA + I - IA 00394 K = J - JA + 1 00395 * 00396 * Compute upper triangular matrix TQ from TAUQ. 00397 * 00398 CALL PZLARFT( 'Forward', 'Columnwise', M-K, JB, A, I+1, J, 00399 $ DESCA, TAUQ, WORK( IPTQ ), WORK( IPWK ) ) 00400 * 00401 * Copy Householder vectors into workspace. 00402 * 00403 CALL PZLACPY( 'Lower', M-K, JB, A, I+1, J, DESCA, 00404 $ WORK( IPV ), IV+1, JV, DESCV ) 00405 CALL PZLASET( 'Upper', M-K, JB, ZERO, ONE, WORK( IPV ), 00406 $ IV+1, JV, DESCV ) 00407 * 00408 * Zero out the strict lower triangular part of A. 00409 * 00410 CALL PZLASET( 'Lower', M-K-1, JB, ZERO, ZERO, A, I+2, J, 00411 $ DESCA ) 00412 * 00413 * Compute upper triangular matrix TP from TAUP. 00414 * 00415 CALL PZLARFT( 'Forward', 'Rowwise', N-K+1, JB, A, I, J, 00416 $ DESCA, TAUP, WORK( IPTP ), WORK( IPWK ) ) 00417 * 00418 * Copy Householder vectors into workspace. 00419 * 00420 CALL PZLACPY( 'Upper', JB, N-K+1, A, I, J, DESCA, 00421 $ WORK( IPW ), IV, JV, DESCW ) 00422 CALL PZLASET( 'Lower', JB, N-K+1, ZERO, ONE, WORK( IPW ), 00423 $ IV, JV, DESCW ) 00424 * 00425 * Zero out the strict+1 upper triangular part of A. 00426 * 00427 CALL PZLASET( 'Upper', JB, N-K, ZERO, ZERO, A, I, J+1, 00428 $ DESCA ) 00429 * 00430 * Apply block Householder transformation from Left. 00431 * 00432 CALL PZLARFB( 'Left', 'No transpose', 'Forward', 00433 $ 'Columnwise', M-K, N-K+1, JB, WORK( IPV ), 00434 $ IV+1, JV, DESCV, WORK( IPTQ ), A, I+1, J, 00435 $ DESCA, WORK( IPWK ) ) 00436 * 00437 * Apply block Householder transformation from Right. 00438 * 00439 CALL PZLARFB( 'Right', 'Conjugate transpose', 'Forward', 00440 $ 'Rowwise', M-K+1, N-K+1, JB, WORK( IPW ), IV, 00441 $ JV, DESCW, WORK( IPTP ), A, I, J, DESCA, 00442 $ WORK( IPWK ) ) 00443 * 00444 DESCV( M_ ) = DESCV( M_ ) + NB 00445 DESCV( RSRC_ ) = MOD( DESCV( RSRC_ ) + NPROW - 1, NPROW ) 00446 DESCV( CSRC_ ) = MOD( DESCV( CSRC_ ) + NPCOL - 1, NPCOL ) 00447 DESCW( N_ ) = DESCW( N_ ) + NB 00448 DESCW( RSRC_ ) = DESCV( RSRC_ ) 00449 DESCW( CSRC_ ) = DESCV( CSRC_ ) 00450 * 00451 40 CONTINUE 00452 * 00453 * Handle first block separately 00454 * 00455 JB = MIN( M, NB - IOFF ) 00456 IV = IOFF + 1 00457 JV = IOFF + 1 00458 * 00459 * Compute upper triangular matrix TQ from TAUQ. 00460 * 00461 CALL PZLARFT( 'Forward', 'Columnwise', M-1, JB, A, IA+1, JA, 00462 $ DESCA, TAUQ, WORK( IPTQ ), WORK( IPWK ) ) 00463 * 00464 * Copy Householder vectors into workspace. 00465 * 00466 CALL PZLACPY( 'Lower', M-1, JB, A, IA+1, JA, DESCA, 00467 $ WORK( IPV ), IV+1, JV, DESCV ) 00468 CALL PZLASET( 'Upper', M-1, JB, ZERO, ONE, WORK( IPV ), IV+1, 00469 $ JV, DESCV ) 00470 * 00471 * Zero out the strict lower triangular part of A. 00472 * 00473 CALL PZLASET( 'Lower', M-2, JB, ZERO, ZERO, A, IA+2, JA, 00474 $ DESCA ) 00475 * 00476 * Compute upper triangular matrix TP from TAUP. 00477 * 00478 CALL PZLARFT( 'Forward', 'Rowwise', N, JB, A, IA, JA, DESCA, 00479 $ TAUP, WORK( IPTP ), WORK( IPWK ) ) 00480 * 00481 * Copy Householder vectors into workspace. 00482 * 00483 CALL PZLACPY( 'Upper', JB, N, A, IA, JA, DESCA, WORK( IPW ), 00484 $ IV, JV, DESCW ) 00485 CALL PZLASET( 'Lower', JB, N, ZERO, ONE, WORK( IPW ), IV, JV, 00486 $ DESCW ) 00487 * 00488 * Zero out the strict+1 upper triangular part of A. 00489 * 00490 CALL PZLASET( 'Upper', JB, N-1, ZERO, ZERO, A, IA, JA+1, 00491 $ DESCA ) 00492 * 00493 * Apply block Householder transformation from left 00494 * 00495 CALL PZLARFB( 'Left', 'No transpose', 'Forward', 'Columnwise', 00496 $ M-1, N, JB, WORK( IPV ), IV+1, JV, DESCV, 00497 $ WORK( IPTQ ), A, IA+1, JA, DESCA, WORK( IPWK ) ) 00498 * 00499 * Apply block Householder transformation from right 00500 * 00501 CALL PZLARFB( 'Right', 'Conjugate transpose', 'Forward', 00502 $ 'Rowwise', M, N, JB, WORK( IPW ), IV, JV, DESCW, 00503 $ WORK( IPTP ), A, IA, JA, DESCA, WORK( IPWK ) ) 00504 END IF 00505 * 00506 CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, INFO, 1, -1, 0 ) 00507 * 00508 RETURN 00509 * 00510 * End of PZGEBDRV 00511 * 00512 END
