1 SUBROUTINE PSLABRD( M, N, NB, A, IA, JA, DESCA, D, E, TAUQ, TAUP, 2 $ X, IX, JX, DESCX, Y, IY, JY, DESCY, WORK ) 3* 4* -- ScaLAPACK auxiliary routine (version 1.7) -- 5* University of Tennessee, Knoxville, Oak Ridge National Laboratory, 6* and University of California, Berkeley. 7* May 1, 1997 8* 9* .. Scalar Arguments .. 10 INTEGER IA, IX, IY, JA, JX, JY, M, N, NB 11* .. 12* .. Array Arguments .. 13 INTEGER DESCA( * ), DESCX( * ), DESCY( * ) 14 REAL A( * ), D( * ), E( * ), TAUP( * ), 15 $ TAUQ( * ), X( * ), Y( * ), WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* PSLABRD reduces the first NB rows and columns of a real general 22* M-by-N distributed matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1) to upper 23* or lower bidiagonal form by an orthogonal transformation Q' * A * P, 24* and returns the matrices X and Y which are needed to apply the 25* transformation to the unreduced part of sub( A ). 26* 27* If M >= N, sub( A ) is reduced to upper bidiagonal form; if M < N, to 28* lower bidiagonal form. 29* 30* This is an auxiliary routine called by PSGEBRD. 31* 32* Notes 33* ===== 34* 35* Each global data object is described by an associated description 36* vector. This vector stores the information required to establish 37* the mapping between an object element and its corresponding process 38* and memory location. 39* 40* Let A be a generic term for any 2D block cyclicly distributed array. 41* Such a global array has an associated description vector DESCA. 42* In the following comments, the character _ should be read as 43* "of the global array". 44* 45* NOTATION STORED IN EXPLANATION 46* --------------- -------------- -------------------------------------- 47* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 48* DTYPE_A = 1. 49* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 50* the BLACS process grid A is distribu- 51* ted over. The context itself is glo- 52* bal, but the handle (the integer 53* value) may vary. 54* M_A (global) DESCA( M_ ) The number of rows in the global 55* array A. 56* N_A (global) DESCA( N_ ) The number of columns in the global 57* array A. 58* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 59* the rows of the array. 60* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 61* the columns of the array. 62* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 63* row of the array A is distributed. 64* CSRC_A (global) DESCA( CSRC_ ) The process column over which the 65* first column of the array A is 66* distributed. 67* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 68* array. LLD_A >= MAX(1,LOCr(M_A)). 69* 70* Let K be the number of rows or columns of a distributed matrix, 71* and assume that its process grid has dimension p x q. 72* LOCr( K ) denotes the number of elements of K that a process 73* would receive if K were distributed over the p processes of its 74* process column. 75* Similarly, LOCc( K ) denotes the number of elements of K that a 76* process would receive if K were distributed over the q processes of 77* its process row. 78* The values of LOCr() and LOCc() may be determined via a call to the 79* ScaLAPACK tool function, NUMROC: 80* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 81* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 82* An upper bound for these quantities may be computed by: 83* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 84* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 85* 86* Arguments 87* ========= 88* 89* M (global input) INTEGER 90* The number of rows to be operated on, i.e. the number of rows 91* of the distributed submatrix sub( A ). M >= 0. 92* 93* N (global input) INTEGER 94* The number of columns to be operated on, i.e. the number of 95* columns of the distributed submatrix sub( A ). N >= 0. 96* 97* NB (global input) INTEGER 98* The number of leading rows and columns of sub( A ) to be 99* reduced. 100* 101* A (local input/local output) REAL pointer into the 102* local memory to an array of dimension (LLD_A,LOCc(JA+N-1)). 103* On entry, this array contains the local pieces of the 104* general distributed matrix sub( A ) to be reduced. On exit, 105* the first NB rows and columns of the matrix are overwritten; 106* the rest of the distributed matrix sub( A ) is unchanged. 107* If m >= n, elements on and below the diagonal in the first NB 108* columns, with the array TAUQ, represent the orthogonal 109* matrix Q as a product of elementary reflectors; and 110* elements above the diagonal in the first NB rows, with the 111* array TAUP, represent the orthogonal matrix P as a product 112* of elementary reflectors. 113* If m < n, elements below the diagonal in the first NB 114* columns, with the array TAUQ, represent the orthogonal 115* matrix Q as a product of elementary reflectors, and 116* elements on and above the diagonal in the first NB rows, 117* with the array TAUP, represent the orthogonal matrix P as 118* a product of elementary reflectors. 119* See Further Details. 120* 121* IA (global input) INTEGER 122* The row index in the global array A indicating the first 123* row of sub( A ). 124* 125* JA (global input) INTEGER 126* The column index in the global array A indicating the 127* first column of sub( A ). 128* 129* DESCA (global and local input) INTEGER array of dimension DLEN_. 130* The array descriptor for the distributed matrix A. 131* 132* D (local output) REAL array, dimension 133* LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-1) otherwise. 134* The distributed diagonal elements of the bidiagonal matrix 135* B: D(i) = A(ia+i-1,ja+i-1). D is tied to the distributed 136* matrix A. 137* 138* E (local output) REAL array, dimension 139* LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-2) otherwise. 140* The distributed off-diagonal elements of the bidiagonal 141* distributed matrix B: 142* if m >= n, E(i) = A(ia+i-1,ja+i) for i = 1,2,...,n-1; 143* if m < n, E(i) = A(ia+i,ja+i-1) for i = 1,2,...,m-1. 144* E is tied to the distributed matrix A. 145* 146* TAUQ (local output) REAL array dimension 147* LOCc(JA+MIN(M,N)-1). The scalar factors of the elementary 148* reflectors which represent the orthogonal matrix Q. TAUQ 149* is tied to the distributed matrix A. See Further Details. 150* 151* TAUP (local output) REAL array, dimension 152* LOCr(IA+MIN(M,N)-1). The scalar factors of the elementary 153* reflectors which represent the orthogonal matrix P. TAUP 154* is tied to the distributed matrix A. See Further Details. 155* 156* X (local output) REAL pointer into the local memory 157* to an array of dimension (LLD_X,NB). On exit, the local 158* pieces of the distributed M-by-NB matrix 159* X(IX:IX+M-1,JX:JX+NB-1) required to update the unreduced 160* part of sub( A ). 161* 162* IX (global input) INTEGER 163* The row index in the global array X indicating the first 164* row of sub( X ). 165* 166* JX (global input) INTEGER 167* The column index in the global array X indicating the 168* first column of sub( X ). 169* 170* DESCX (global and local input) INTEGER array of dimension DLEN_. 171* The array descriptor for the distributed matrix X. 172* 173* Y (local output) REAL pointer into the local memory 174* to an array of dimension (LLD_Y,NB). On exit, the local 175* pieces of the distributed N-by-NB matrix 176* Y(IY:IY+N-1,JY:JY+NB-1) required to update the unreduced 177* part of sub( A ). 178* 179* IY (global input) INTEGER 180* The row index in the global array Y indicating the first 181* row of sub( Y ). 182* 183* JY (global input) INTEGER 184* The column index in the global array Y indicating the 185* first column of sub( Y ). 186* 187* DESCY (global and local input) INTEGER array of dimension DLEN_. 188* The array descriptor for the distributed matrix Y. 189* 190* WORK (local workspace) REAL array, dimension (LWORK) 191* LWORK >= NB_A + NQ, with 192* 193* NQ = NUMROC( N+MOD( IA-1, NB_Y ), NB_Y, MYCOL, IACOL, NPCOL ) 194* IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ) 195* 196* INDXG2P and NUMROC are ScaLAPACK tool functions; 197* MYROW, MYCOL, NPROW and NPCOL can be determined by calling 198* the subroutine BLACS_GRIDINFO. 199* 200* Further Details 201* =============== 202* 203* The matrices Q and P are represented as products of elementary 204* reflectors: 205* 206* Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb) 207* 208* Each H(i) and G(i) has the form: 209* 210* H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' 211* 212* where tauq and taup are real scalars, and v and u are real vectors. 213* 214* If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in 215* A(ia+i-1:ia+m-1,ja+i-1); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is 216* stored on exit in A(ia+i-1,ja+i:ja+n-1); tauq is stored in 217* TAUQ(ja+i-1) and taup in TAUP(ia+i-1). 218* 219* If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in 220* A(ia+i+1:ia+m-1,ja+i-1); u(1:i-1) = 0, u(i) = 1, and u(i:n) is 221* stored on exit in A(ia+i-1,ja+i:ja+n-1); tauq is stored in 222* TAUQ(ja+i-1) and taup in TAUP(ia+i-1). 223* 224* The elements of the vectors v and u together form the m-by-nb matrix 225* V and the nb-by-n matrix U' which are needed, with X and Y, to apply 226* the transformation to the unreduced part of the matrix, using a block 227* update of the form: sub( A ) := sub( A ) - V*Y' - X*U'. 228* 229* The contents of sub( A ) on exit are illustrated by the following 230* examples with nb = 2: 231* 232* m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): 233* 234* ( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 ) 235* ( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 ) 236* ( v1 v2 a a a ) ( v1 1 a a a a ) 237* ( v1 v2 a a a ) ( v1 v2 a a a a ) 238* ( v1 v2 a a a ) ( v1 v2 a a a a ) 239* ( v1 v2 a a a ) 240* 241* where a denotes an element of the original matrix which is unchanged, 242* vi denotes an element of the vector defining H(i), and ui an element 243* of the vector defining G(i). 244* 245* ===================================================================== 246* 247* .. Parameters .. 248 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 249 $ LLD_, MB_, M_, NB_, N_, RSRC_ 250 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 251 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 252 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 253 REAL ONE, ZERO 254 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 255* .. 256* .. Local Scalars .. 257 INTEGER I, IACOL, IAROW, ICTXT, II, IPY, IW, J, JJ, 258 $ JWY, K, MYCOL, MYROW, NPCOL, NPROW 259 REAL ALPHA, TAU 260 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ), 261 $ DESCTP( DLEN_ ), DESCTQ( DLEN_ ), 262 $ DESCW( DLEN_ ), DESCWY( DLEN_ ) 263* .. 264* .. External Subroutines .. 265 EXTERNAL BLACS_GRIDINFO, DESCSET, INFOG2L, PSCOPY, 266 $ PSELGET, PSELSET, PSGEMV, PSLARFG, 267 $ PSSCAL 268* .. 269* .. Intrinsic Functions .. 270 INTRINSIC MIN, MOD 271* .. 272* .. Executable Statements .. 273* 274* Quick return if possible 275* 276 IF( M.LE.0 .OR. N.LE.0 ) 277 $ RETURN 278* 279 ICTXT = DESCA( CTXT_ ) 280 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 281 CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, II, JJ, 282 $ IAROW, IACOL ) 283 IPY = DESCA( MB_ ) + 1 284 IW = MOD( IA-1, DESCA( NB_ ) ) + 1 285 ALPHA = ZERO 286* 287 CALL DESCSET( DESCWY, 1, N+MOD( IA-1, DESCY( NB_ ) ), 1, 288 $ DESCA( NB_ ), IAROW, IACOL, ICTXT, 1 ) 289 CALL DESCSET( DESCW, DESCA( MB_ ), 1, DESCA( MB_ ), 1, IAROW, 290 $ IACOL, ICTXT, DESCA( MB_ ) ) 291 CALL DESCSET( DESCTQ, 1, JA+MIN(M,N)-1, 1, DESCA( NB_ ), IAROW, 292 $ DESCA( CSRC_ ), DESCA( CTXT_ ), 1 ) 293 CALL DESCSET( DESCTP, IA+MIN(M,N)-1, 1, DESCA( MB_ ), 1, 294 $ DESCA( RSRC_ ), IACOL, DESCA( CTXT_ ), 295 $ DESCA( LLD_ ) ) 296* 297 IF( M.GE.N ) THEN 298* 299* Reduce to upper bidiagonal form 300* 301 CALL DESCSET( DESCD, 1, JA+MIN(M,N)-1, 1, DESCA( NB_ ), MYROW, 302 $ DESCA( CSRC_ ), DESCA( CTXT_ ), 1 ) 303 CALL DESCSET( DESCE, IA+MIN(M,N)-1, 1, DESCA( MB_ ), 1, 304 $ DESCA( RSRC_ ), MYCOL, DESCA( CTXT_ ), 305 $ DESCA( LLD_ ) ) 306 DO 10 K = 1, NB 307 I = IA + K - 1 308 J = JA + K - 1 309 JWY = IW + K 310* 311* Update A(i:ia+m-1,j) 312* 313 IF( K.GT.1 ) THEN 314 CALL PSGEMV( 'No transpose', M-K+1, K-1, -ONE, A, I, JA, 315 $ DESCA, Y, IY, JY+K-1, DESCY, 1, ONE, A, I, 316 $ J, DESCA, 1 ) 317 CALL PSGEMV( 'No transpose', M-K+1, K-1, -ONE, X, IX+K-1, 318 $ JX, DESCX, A, IA, J, DESCA, 1, ONE, A, I, J, 319 $ DESCA, 1 ) 320 CALL PSELSET( A, I-1, J, DESCA, ALPHA ) 321 END IF 322* 323* Generate reflection Q(i) to annihilate A(i+1:ia+m-1,j) 324* 325 CALL PSLARFG( M-K+1, ALPHA, I, J, A, I+1, J, DESCA, 1, 326 $ TAUQ ) 327 CALL PSELSET( D, 1, J, DESCD, ALPHA ) 328 CALL PSELSET( A, I, J, DESCA, ONE ) 329* 330* Compute Y(IA+I:IA+N-1,J) 331* 332 CALL PSGEMV( 'Transpose', M-K+1, N-K, ONE, A, I, J+1, DESCA, 333 $ A, I, J, DESCA, 1, ZERO, WORK( IPY ), 1, JWY, 334 $ DESCWY, DESCWY( M_ ) ) 335 CALL PSGEMV( 'Transpose', M-K+1, K-1, ONE, A, I, JA, DESCA, 336 $ A, I, J, DESCA, 1, ZERO, WORK, IW, 1, DESCW, 337 $ 1 ) 338 CALL PSGEMV( 'Transpose', K-1, N-K, -ONE, Y, IY, JY+K, 339 $ DESCY, WORK, IW, 1, DESCW, 1, ONE, WORK( IPY ), 340 $ 1, JWY, DESCWY, DESCWY( M_ ) ) 341 CALL PSGEMV( 'Transpose', M-K+1, K-1, ONE, X, IX+K-1, JX, 342 $ DESCX, A, I, J, DESCA, 1, ZERO, WORK, IW, 1, 343 $ DESCW, 1 ) 344 CALL PSGEMV( 'Transpose', K-1, N-K, -ONE, A, IA, J+1, DESCA, 345 $ WORK, IW, 1, DESCW, 1, ONE, WORK( IPY ), 1, 346 $ JWY, DESCWY, DESCWY( M_ ) ) 347* 348 CALL PSELGET( 'Rowwise', ' ', TAU, TAUQ, 1, J, DESCTQ ) 349 CALL PSSCAL( N-K, TAU, WORK( IPY ), 1, JWY, DESCWY, 350 $ DESCWY( M_ ) ) 351 CALL PSCOPY( N-K, WORK( IPY ), 1, JWY, DESCWY, DESCWY( M_ ), 352 $ Y, IY+K-1, JY+K, DESCY, DESCY( M_ ) ) 353* 354* Update A(i,j+1:ja+n-1) 355* 356 CALL PSGEMV( 'Transpose', K, N-K, -ONE, Y, IY, JY+K, DESCY, 357 $ A, I, JA, DESCA, DESCA( M_ ), ONE, A, I, J+1, 358 $ DESCA, DESCA( M_ ) ) 359 CALL PSGEMV( 'Transpose', K-1, N-K, -ONE, A, IA, J+1, DESCA, 360 $ X, IX+K-1, JX, DESCX, DESCX( M_ ), ONE, A, I, 361 $ J+1, DESCA, DESCA( M_ ) ) 362 CALL PSELSET( A, I, J, DESCA, ALPHA ) 363* 364* Generate reflection P(i) to annihilate A(i,j+2:ja+n-1) 365* 366 CALL PSLARFG( N-K, ALPHA, I, J+1, A, I, 367 $ MIN( J+2, N+JA-1 ), DESCA, DESCA( M_ ), TAUP ) 368 CALL PSELSET( E, I, 1, DESCE, ALPHA ) 369 CALL PSELSET( A, I, J+1, DESCA, ONE ) 370* 371* Compute X(I+1:IA+M-1,J) 372* 373 CALL PSGEMV( 'No transpose', M-K, N-K, ONE, A, I+1, J+1, 374 $ DESCA, A, I, J+1, DESCA, DESCA( M_ ), ZERO, X, 375 $ IX+K, JX+K-1, DESCX, 1 ) 376 CALL PSGEMV( 'No transpose', K, N-K, ONE, Y, IY, JY+K, 377 $ DESCY, A, I, J+1, DESCA, DESCA( M_ ), ZERO, 378 $ WORK, IW, 1, DESCW, 1 ) 379 CALL PSGEMV( 'No transpose', M-K, K, -ONE, A, I+1, JA, 380 $ DESCA, WORK, IW, 1, DESCW, 1, ONE, X, IX+K, 381 $ JX+K-1, DESCX, 1 ) 382 CALL PSGEMV( 'No transpose', K-1, N-K, ONE, A, IA, J+1, 383 $ DESCA, A, I, J+1, DESCA, DESCA( M_ ), ZERO, 384 $ WORK, IW, 1, DESCW, 1 ) 385 CALL PSGEMV( 'No transpose', M-K, K-1, -ONE, X, IX+K, JX, 386 $ DESCX, WORK, IW, 1, DESCW, 1, ONE, X, IX+K, 387 $ JX+K-1, DESCX, 1 ) 388* 389 CALL PSELGET( 'Columnwise', ' ', TAU, TAUP, I, 1, DESCTP ) 390 CALL PSSCAL( M-K, TAU, X, IX+K, JX+K-1, DESCX, 1 ) 391 10 CONTINUE 392* 393 ELSE 394* 395* Reduce to lower bidiagonal form 396* 397 CALL DESCSET( DESCD, IA+MIN(M,N)-1, 1, DESCA( MB_ ), 1, 398 $ DESCA( RSRC_ ), MYCOL, DESCA( CTXT_ ), 399 $ DESCA( LLD_ ) ) 400 CALL DESCSET( DESCE, 1, JA+MIN(M,N)-1, 1, DESCA( NB_ ), MYROW, 401 $ DESCA( CSRC_ ), DESCA( CTXT_ ), 1 ) 402 DO 20 K = 1, NB 403 I = IA + K - 1 404 J = JA + K - 1 405 JWY = IW + K 406* 407* Update A(i,j:ja+n-1) 408* 409 IF( K.GT.1 ) THEN 410 CALL PSGEMV( 'Transpose', K-1, N-K+1, -ONE, Y, IY, 411 $ JY+K-1, DESCY, A, I, JA, DESCA, DESCA( M_ ), 412 $ ONE, A, I, J, DESCA, DESCA( M_ ) ) 413 CALL PSGEMV( 'Transpose', K-1, N-K+1, -ONE, A, IA, J, 414 $ DESCA, X, IX+K-1, JX, DESCX, DESCX( M_ ), 415 $ ONE, A, I, J, DESCA, DESCA( M_ ) ) 416 CALL PSELSET( A, I, J-1, DESCA, ALPHA ) 417 END IF 418* 419* Generate reflection P(i) to annihilate A(i,j+1:ja+n-1) 420* 421 CALL PSLARFG( N-K+1, ALPHA, I, J, A, I, J+1, DESCA, 422 $ DESCA( M_ ), TAUP ) 423 CALL PSELSET( D, I, 1, DESCD, ALPHA ) 424 CALL PSELSET( A, I, J, DESCA, ONE ) 425* 426* Compute X(i+1:ia+m-1,j) 427* 428 CALL PSGEMV( 'No transpose', M-K, N-K+1, ONE, A, I+1, J, 429 $ DESCA, A, I, J, DESCA, DESCA( M_ ), ZERO, X, 430 $ IX+K, JX+K-1, DESCX, 1 ) 431 CALL PSGEMV( 'No transpose', K-1, N-K+1, ONE, Y, IY, JY+K-1, 432 $ DESCY, A, I, J, DESCA, DESCA( M_ ), ZERO, 433 $ WORK, IW, 1, DESCW, 1 ) 434 CALL PSGEMV( 'No transpose', M-K, K-1, -ONE, A, I+1, JA, 435 $ DESCA, WORK, IW, 1, DESCW, 1, ONE, X, IX+K, 436 $ JX+K-1, DESCX, 1 ) 437 CALL PSGEMV( 'No transpose', K-1, N-K+1, ONE, A, IA, J, 438 $ DESCA, A, I, J, DESCA, DESCA( M_ ), ZERO, 439 $ WORK, IW, 1, DESCW, 1 ) 440 CALL PSGEMV( 'No transpose', M-K, K-1, -ONE, X, IX+K, JX, 441 $ DESCX, WORK, IW, 1, DESCW, 1, ONE, X, IX+K, 442 $ JX+K-1, DESCX, 1 ) 443* 444 CALL PSELGET( 'Columnwise', ' ', TAU, TAUP, I, 1, DESCTP ) 445 CALL PSSCAL( M-K, TAU, X, IX+K, JX+K-1, DESCX, 1 ) 446* 447* Update A(i+1:ia+m-1,j) 448* 449 CALL PSGEMV( 'No transpose', M-K, K-1, -ONE, A, I+1, JA, 450 $ DESCA, Y, IY, JY+K-1, DESCY, 1, ONE, A, I+1, J, 451 $ DESCA, 1 ) 452 CALL PSGEMV( 'No transpose', M-K, K, -ONE, X, IX+K, JX, 453 $ DESCX, A, IA, J, DESCA, 1, ONE, A, I+1, J, 454 $ DESCA, 1 ) 455 CALL PSELSET( A, I, J, DESCA, ALPHA ) 456* 457* Generate reflection Q(i) to annihilate A(i+2:ia+m-1,j) 458* 459 CALL PSLARFG( M-K, ALPHA, I+1, J, A, MIN( I+2, M+IA-1 ), 460 $ J, DESCA, 1, TAUQ ) 461 CALL PSELSET( E, 1, J, DESCE, ALPHA ) 462 CALL PSELSET( A, I+1, J, DESCA, ONE ) 463* 464* Compute Y(ia+i:ia+n-1,j) 465* 466 CALL PSGEMV( 'Transpose', M-K, N-K, ONE, A, I+1, J+1, DESCA, 467 $ A, I+1, J, DESCA, 1, ZERO, WORK( IPY ), 1, 468 $ JWY, DESCWY, DESCWY( M_ ) ) 469 CALL PSGEMV( 'Transpose', M-K, K-1, ONE, A, I+1, JA, DESCA, 470 $ A, I+1, J, DESCA, 1, ZERO, WORK, IW, 1, DESCW, 471 $ 1 ) 472 CALL PSGEMV( 'Transpose', K-1, N-K, -ONE, Y, IY, JY+K, 473 $ DESCY, WORK, IW, 1, DESCW, 1, ONE, WORK( IPY ), 474 $ 1, JWY, DESCWY, DESCWY( M_ ) ) 475 CALL PSGEMV( 'Transpose', M-K, K, ONE, X, IX+K, JX, DESCX, 476 $ A, I+1, J, DESCA, 1, ZERO, WORK, IW, 1, DESCW, 477 $ 1 ) 478 CALL PSGEMV( 'Transpose', K, N-K, -ONE, A, IA, J+1, DESCA, 479 $ WORK, IW, 1, DESCW, 1, ONE, WORK( IPY ), 1, 480 $ JWY, DESCWY, DESCWY( M_ ) ) 481* 482 CALL PSELGET( 'Rowwise', ' ', TAU, TAUQ, 1, J, DESCTQ ) 483 CALL PSSCAL( N-K, TAU, WORK( IPY ), 1, JWY, DESCWY, 484 $ DESCWY( M_ ) ) 485 CALL PSCOPY( N-K, WORK( IPY ), 1, JWY, DESCWY, DESCWY( M_ ), 486 $ Y, IY+K-1, JY+K, DESCY, DESCY( M_ ) ) 487 20 CONTINUE 488 END IF 489* 490 RETURN 491* 492* End of PSLABRD 493* 494 END 495