|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE PSLABRD( M, N, NB, A, IA, JA, DESCA, D, E, TAUQ, TAUP, 00002 $ X, IX, JX, DESCX, Y, IY, JY, DESCY, WORK ) 00003 * 00004 * -- ScaLAPACK auxiliary 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 IA, IX, IY, JA, JX, JY, M, N, NB 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER DESCA( * ), DESCX( * ), DESCY( * ) 00014 REAL A( * ), D( * ), E( * ), TAUP( * ), 00015 $ TAUQ( * ), X( * ), Y( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * PSLABRD reduces the first NB rows and columns of a real general 00022 * M-by-N distributed matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1) to upper 00023 * or lower bidiagonal form by an orthogonal transformation Q' * A * P, 00024 * and returns the matrices X and Y which are needed to apply the 00025 * transformation to the unreduced part of sub( A ). 00026 * 00027 * If M >= N, sub( A ) is reduced to upper bidiagonal form; if M < N, to 00028 * lower bidiagonal form. 00029 * 00030 * This is an auxiliary routine called by PSGEBRD. 00031 * 00032 * Notes 00033 * ===== 00034 * 00035 * Each global data object is described by an associated description 00036 * vector. This vector stores the information required to establish 00037 * the mapping between an object element and its corresponding process 00038 * and memory location. 00039 * 00040 * Let A be a generic term for any 2D block cyclicly distributed array. 00041 * Such a global array has an associated description vector DESCA. 00042 * In the following comments, the character _ should be read as 00043 * "of the global array". 00044 * 00045 * NOTATION STORED IN EXPLANATION 00046 * --------------- -------------- -------------------------------------- 00047 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 00048 * DTYPE_A = 1. 00049 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 00050 * the BLACS process grid A is distribu- 00051 * ted over. The context itself is glo- 00052 * bal, but the handle (the integer 00053 * value) may vary. 00054 * M_A (global) DESCA( M_ ) The number of rows in the global 00055 * array A. 00056 * N_A (global) DESCA( N_ ) The number of columns in the global 00057 * array A. 00058 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 00059 * the rows of the array. 00060 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 00061 * the columns of the array. 00062 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 00063 * row of the array A is distributed. 00064 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the 00065 * first column of the array A is 00066 * distributed. 00067 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 00068 * array. LLD_A >= MAX(1,LOCr(M_A)). 00069 * 00070 * Let K be the number of rows or columns of a distributed matrix, 00071 * and assume that its process grid has dimension p x q. 00072 * LOCr( K ) denotes the number of elements of K that a process 00073 * would receive if K were distributed over the p processes of its 00074 * process column. 00075 * Similarly, LOCc( K ) denotes the number of elements of K that a 00076 * process would receive if K were distributed over the q processes of 00077 * its process row. 00078 * The values of LOCr() and LOCc() may be determined via a call to the 00079 * ScaLAPACK tool function, NUMROC: 00080 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 00081 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 00082 * An upper bound for these quantities may be computed by: 00083 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 00084 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 00085 * 00086 * Arguments 00087 * ========= 00088 * 00089 * M (global input) INTEGER 00090 * The number of rows to be operated on, i.e. the number of rows 00091 * of the distributed submatrix sub( A ). M >= 0. 00092 * 00093 * N (global input) INTEGER 00094 * The number of columns to be operated on, i.e. the number of 00095 * columns of the distributed submatrix sub( A ). N >= 0. 00096 * 00097 * NB (global input) INTEGER 00098 * The number of leading rows and columns of sub( A ) to be 00099 * reduced. 00100 * 00101 * A (local input/local output) REAL pointer into the 00102 * local memory to an array of dimension (LLD_A,LOCc(JA+N-1)). 00103 * On entry, this array contains the local pieces of the 00104 * general distributed matrix sub( A ) to be reduced. On exit, 00105 * the first NB rows and columns of the matrix are overwritten; 00106 * the rest of the distributed matrix sub( A ) is unchanged. 00107 * If m >= n, elements on and below the diagonal in the first NB 00108 * columns, with the array TAUQ, represent the orthogonal 00109 * matrix Q as a product of elementary reflectors; and 00110 * elements above the diagonal in the first NB rows, with the 00111 * array TAUP, represent the orthogonal matrix P as a product 00112 * of elementary reflectors. 00113 * If m < n, elements below the diagonal in the first NB 00114 * columns, with the array TAUQ, represent the orthogonal 00115 * matrix Q as a product of elementary reflectors, and 00116 * elements on and above the diagonal in the first NB rows, 00117 * with the array TAUP, represent the orthogonal matrix P as 00118 * a product of elementary reflectors. 00119 * See Further Details. 00120 * 00121 * IA (global input) INTEGER 00122 * The row index in the global array A indicating the first 00123 * row of sub( A ). 00124 * 00125 * JA (global input) INTEGER 00126 * The column index in the global array A indicating the 00127 * first column of sub( A ). 00128 * 00129 * DESCA (global and local input) INTEGER array of dimension DLEN_. 00130 * The array descriptor for the distributed matrix A. 00131 * 00132 * D (local output) REAL array, dimension 00133 * LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-1) otherwise. 00134 * The distributed diagonal elements of the bidiagonal matrix 00135 * B: D(i) = A(ia+i-1,ja+i-1). D is tied to the distributed 00136 * matrix A. 00137 * 00138 * E (local output) REAL array, dimension 00139 * LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-2) otherwise. 00140 * The distributed off-diagonal elements of the bidiagonal 00141 * distributed matrix B: 00142 * if m >= n, E(i) = A(ia+i-1,ja+i) for i = 1,2,...,n-1; 00143 * if m < n, E(i) = A(ia+i,ja+i-1) for i = 1,2,...,m-1. 00144 * E is tied to the distributed matrix A. 00145 * 00146 * TAUQ (local output) REAL array dimension 00147 * LOCc(JA+MIN(M,N)-1). The scalar factors of the elementary 00148 * reflectors which represent the orthogonal matrix Q. TAUQ 00149 * is tied to the distributed matrix A. See Further Details. 00150 * 00151 * TAUP (local output) REAL array, dimension 00152 * LOCr(IA+MIN(M,N)-1). The scalar factors of the elementary 00153 * reflectors which represent the orthogonal matrix P. TAUP 00154 * is tied to the distributed matrix A. See Further Details. 00155 * 00156 * X (local output) REAL pointer into the local memory 00157 * to an array of dimension (LLD_X,NB). On exit, the local 00158 * pieces of the distributed M-by-NB matrix 00159 * X(IX:IX+M-1,JX:JX+NB-1) required to update the unreduced 00160 * part of sub( A ). 00161 * 00162 * IX (global input) INTEGER 00163 * The row index in the global array X indicating the first 00164 * row of sub( X ). 00165 * 00166 * JX (global input) INTEGER 00167 * The column index in the global array X indicating the 00168 * first column of sub( X ). 00169 * 00170 * DESCX (global and local input) INTEGER array of dimension DLEN_. 00171 * The array descriptor for the distributed matrix X. 00172 * 00173 * Y (local output) REAL pointer into the local memory 00174 * to an array of dimension (LLD_Y,NB). On exit, the local 00175 * pieces of the distributed N-by-NB matrix 00176 * Y(IY:IY+N-1,JY:JY+NB-1) required to update the unreduced 00177 * part of sub( A ). 00178 * 00179 * IY (global input) INTEGER 00180 * The row index in the global array Y indicating the first 00181 * row of sub( Y ). 00182 * 00183 * JY (global input) INTEGER 00184 * The column index in the global array Y indicating the 00185 * first column of sub( Y ). 00186 * 00187 * DESCY (global and local input) INTEGER array of dimension DLEN_. 00188 * The array descriptor for the distributed matrix Y. 00189 * 00190 * WORK (local workspace) REAL array, dimension (LWORK) 00191 * LWORK >= NB_A + NQ, with 00192 * 00193 * NQ = NUMROC( N+MOD( IA-1, NB_Y ), NB_Y, MYCOL, IACOL, NPCOL ) 00194 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ) 00195 * 00196 * INDXG2P and NUMROC are ScaLAPACK tool functions; 00197 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling 00198 * the subroutine BLACS_GRIDINFO. 00199 * 00200 * Further Details 00201 * =============== 00202 * 00203 * The matrices Q and P are represented as products of elementary 00204 * reflectors: 00205 * 00206 * Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb) 00207 * 00208 * Each H(i) and G(i) has the form: 00209 * 00210 * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' 00211 * 00212 * where tauq and taup are real scalars, and v and u are real vectors. 00213 * 00214 * If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in 00215 * A(ia+i-1:ia+m-1,ja+i-1); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is 00216 * stored on exit in A(ia+i-1,ja+i:ja+n-1); tauq is stored in 00217 * TAUQ(ja+i-1) and taup in TAUP(ia+i-1). 00218 * 00219 * If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in 00220 * A(ia+i+1:ia+m-1,ja+i-1); u(1:i-1) = 0, u(i) = 1, and u(i:n) is 00221 * stored on exit in A(ia+i-1,ja+i:ja+n-1); tauq is stored in 00222 * TAUQ(ja+i-1) and taup in TAUP(ia+i-1). 00223 * 00224 * The elements of the vectors v and u together form the m-by-nb matrix 00225 * V and the nb-by-n matrix U' which are needed, with X and Y, to apply 00226 * the transformation to the unreduced part of the matrix, using a block 00227 * update of the form: sub( A ) := sub( A ) - V*Y' - X*U'. 00228 * 00229 * The contents of sub( A ) on exit are illustrated by the following 00230 * examples with nb = 2: 00231 * 00232 * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): 00233 * 00234 * ( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 ) 00235 * ( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 ) 00236 * ( v1 v2 a a a ) ( v1 1 a a a a ) 00237 * ( v1 v2 a a a ) ( v1 v2 a a a a ) 00238 * ( v1 v2 a a a ) ( v1 v2 a a a a ) 00239 * ( v1 v2 a a a ) 00240 * 00241 * where a denotes an element of the original matrix which is unchanged, 00242 * vi denotes an element of the vector defining H(i), and ui an element 00243 * of the vector defining G(i). 00244 * 00245 * ===================================================================== 00246 * 00247 * .. Parameters .. 00248 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 00249 $ LLD_, MB_, M_, NB_, N_, RSRC_ 00250 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00251 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00252 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00253 REAL ONE, ZERO 00254 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00255 * .. 00256 * .. Local Scalars .. 00257 INTEGER I, IACOL, IAROW, ICTXT, II, IPY, IW, J, JJ, 00258 $ JWY, K, MYCOL, MYROW, NPCOL, NPROW 00259 REAL ALPHA, TAU 00260 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ), 00261 $ DESCTP( DLEN_ ), DESCTQ( DLEN_ ), 00262 $ DESCW( DLEN_ ), DESCWY( DLEN_ ) 00263 * .. 00264 * .. External Subroutines .. 00265 EXTERNAL BLACS_GRIDINFO, DESCSET, INFOG2L, PSCOPY, 00266 $ PSELGET, PSELSET, PSGEMV, PSLARFG, 00267 $ PSSCAL 00268 * .. 00269 * .. Intrinsic Functions .. 00270 INTRINSIC MIN, MOD 00271 * .. 00272 * .. Executable Statements .. 00273 * 00274 * Quick return if possible 00275 * 00276 IF( M.LE.0 .OR. N.LE.0 ) 00277 $ RETURN 00278 * 00279 ICTXT = DESCA( CTXT_ ) 00280 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 00281 CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, II, JJ, 00282 $ IAROW, IACOL ) 00283 IPY = DESCA( MB_ ) + 1 00284 IW = MOD( IA-1, DESCA( NB_ ) ) + 1 00285 ALPHA = ZERO 00286 * 00287 CALL DESCSET( DESCWY, 1, N+MOD( IA-1, DESCY( NB_ ) ), 1, 00288 $ DESCA( NB_ ), IAROW, IACOL, ICTXT, 1 ) 00289 CALL DESCSET( DESCW, DESCA( MB_ ), 1, DESCA( MB_ ), 1, IAROW, 00290 $ IACOL, ICTXT, DESCA( MB_ ) ) 00291 CALL DESCSET( DESCTQ, 1, JA+MIN(M,N)-1, 1, DESCA( NB_ ), IAROW, 00292 $ DESCA( CSRC_ ), DESCA( CTXT_ ), 1 ) 00293 CALL DESCSET( DESCTP, IA+MIN(M,N)-1, 1, DESCA( MB_ ), 1, 00294 $ DESCA( RSRC_ ), IACOL, DESCA( CTXT_ ), 00295 $ DESCA( LLD_ ) ) 00296 * 00297 IF( M.GE.N ) THEN 00298 * 00299 * Reduce to upper bidiagonal form 00300 * 00301 CALL DESCSET( DESCD, 1, JA+MIN(M,N)-1, 1, DESCA( NB_ ), MYROW, 00302 $ DESCA( CSRC_ ), DESCA( CTXT_ ), 1 ) 00303 CALL DESCSET( DESCE, IA+MIN(M,N)-1, 1, DESCA( MB_ ), 1, 00304 $ DESCA( RSRC_ ), MYCOL, DESCA( CTXT_ ), 00305 $ DESCA( LLD_ ) ) 00306 DO 10 K = 1, NB 00307 I = IA + K - 1 00308 J = JA + K - 1 00309 JWY = IW + K 00310 * 00311 * Update A(i:ia+m-1,j) 00312 * 00313 IF( K.GT.1 ) THEN 00314 CALL PSGEMV( 'No transpose', M-K+1, K-1, -ONE, A, I, JA, 00315 $ DESCA, Y, IY, JY+K-1, DESCY, 1, ONE, A, I, 00316 $ J, DESCA, 1 ) 00317 CALL PSGEMV( 'No transpose', M-K+1, K-1, -ONE, X, IX+K-1, 00318 $ JX, DESCX, A, IA, J, DESCA, 1, ONE, A, I, J, 00319 $ DESCA, 1 ) 00320 CALL PSELSET( A, I-1, J, DESCA, ALPHA ) 00321 END IF 00322 * 00323 * Generate reflection Q(i) to annihilate A(i+1:ia+m-1,j) 00324 * 00325 CALL PSLARFG( M-K+1, ALPHA, I, J, A, I+1, J, DESCA, 1, 00326 $ TAUQ ) 00327 CALL PSELSET( D, 1, J, DESCD, ALPHA ) 00328 CALL PSELSET( A, I, J, DESCA, ONE ) 00329 * 00330 * Compute Y(IA+I:IA+N-1,J) 00331 * 00332 CALL PSGEMV( 'Transpose', M-K+1, N-K, ONE, A, I, J+1, DESCA, 00333 $ A, I, J, DESCA, 1, ZERO, WORK( IPY ), 1, JWY, 00334 $ DESCWY, DESCWY( M_ ) ) 00335 CALL PSGEMV( 'Transpose', M-K+1, K-1, ONE, A, I, JA, DESCA, 00336 $ A, I, J, DESCA, 1, ZERO, WORK, IW, 1, DESCW, 00337 $ 1 ) 00338 CALL PSGEMV( 'Transpose', K-1, N-K, -ONE, Y, IY, JY+K, 00339 $ DESCY, WORK, IW, 1, DESCW, 1, ONE, WORK( IPY ), 00340 $ 1, JWY, DESCWY, DESCWY( M_ ) ) 00341 CALL PSGEMV( 'Transpose', M-K+1, K-1, ONE, X, IX+K-1, JX, 00342 $ DESCX, A, I, J, DESCA, 1, ZERO, WORK, IW, 1, 00343 $ DESCW, 1 ) 00344 CALL PSGEMV( 'Transpose', K-1, N-K, -ONE, A, IA, J+1, DESCA, 00345 $ WORK, IW, 1, DESCW, 1, ONE, WORK( IPY ), 1, 00346 $ JWY, DESCWY, DESCWY( M_ ) ) 00347 * 00348 CALL PSELGET( 'Rowwise', ' ', TAU, TAUQ, 1, J, DESCTQ ) 00349 CALL PSSCAL( N-K, TAU, WORK( IPY ), 1, JWY, DESCWY, 00350 $ DESCWY( M_ ) ) 00351 CALL PSCOPY( N-K, WORK( IPY ), 1, JWY, DESCWY, DESCWY( M_ ), 00352 $ Y, IY+K-1, JY+K, DESCY, DESCY( M_ ) ) 00353 * 00354 * Update A(i,j+1:ja+n-1) 00355 * 00356 CALL PSGEMV( 'Transpose', K, N-K, -ONE, Y, IY, JY+K, DESCY, 00357 $ A, I, JA, DESCA, DESCA( M_ ), ONE, A, I, J+1, 00358 $ DESCA, DESCA( M_ ) ) 00359 CALL PSGEMV( 'Transpose', K-1, N-K, -ONE, A, IA, J+1, DESCA, 00360 $ X, IX+K-1, JX, DESCX, DESCX( M_ ), ONE, A, I, 00361 $ J+1, DESCA, DESCA( M_ ) ) 00362 CALL PSELSET( A, I, J, DESCA, ALPHA ) 00363 * 00364 * Generate reflection P(i) to annihilate A(i,j+2:ja+n-1) 00365 * 00366 CALL PSLARFG( N-K, ALPHA, I, J+1, A, I, 00367 $ MIN( J+2, N+JA-1 ), DESCA, DESCA( M_ ), TAUP ) 00368 CALL PSELSET( E, I, 1, DESCE, ALPHA ) 00369 CALL PSELSET( A, I, J+1, DESCA, ONE ) 00370 * 00371 * Compute X(I+1:IA+M-1,J) 00372 * 00373 CALL PSGEMV( 'No transpose', M-K, N-K, ONE, A, I+1, J+1, 00374 $ DESCA, A, I, J+1, DESCA, DESCA( M_ ), ZERO, X, 00375 $ IX+K, JX+K-1, DESCX, 1 ) 00376 CALL PSGEMV( 'No transpose', K, N-K, ONE, Y, IY, JY+K, 00377 $ DESCY, A, I, J+1, DESCA, DESCA( M_ ), ZERO, 00378 $ WORK, IW, 1, DESCW, 1 ) 00379 CALL PSGEMV( 'No transpose', M-K, K, -ONE, A, I+1, JA, 00380 $ DESCA, WORK, IW, 1, DESCW, 1, ONE, X, IX+K, 00381 $ JX+K-1, DESCX, 1 ) 00382 CALL PSGEMV( 'No transpose', K-1, N-K, ONE, A, IA, J+1, 00383 $ DESCA, A, I, J+1, DESCA, DESCA( M_ ), ZERO, 00384 $ WORK, IW, 1, DESCW, 1 ) 00385 CALL PSGEMV( 'No transpose', M-K, K-1, -ONE, X, IX+K, JX, 00386 $ DESCX, WORK, IW, 1, DESCW, 1, ONE, X, IX+K, 00387 $ JX+K-1, DESCX, 1 ) 00388 * 00389 CALL PSELGET( 'Columnwise', ' ', TAU, TAUP, I, 1, DESCTP ) 00390 CALL PSSCAL( M-K, TAU, X, IX+K, JX+K-1, DESCX, 1 ) 00391 10 CONTINUE 00392 * 00393 ELSE 00394 * 00395 * Reduce to lower bidiagonal form 00396 * 00397 CALL DESCSET( DESCD, IA+MIN(M,N)-1, 1, DESCA( MB_ ), 1, 00398 $ DESCA( RSRC_ ), MYCOL, DESCA( CTXT_ ), 00399 $ DESCA( LLD_ ) ) 00400 CALL DESCSET( DESCE, 1, JA+MIN(M,N)-1, 1, DESCA( NB_ ), MYROW, 00401 $ DESCA( CSRC_ ), DESCA( CTXT_ ), 1 ) 00402 DO 20 K = 1, NB 00403 I = IA + K - 1 00404 J = JA + K - 1 00405 JWY = IW + K 00406 * 00407 * Update A(i,j:ja+n-1) 00408 * 00409 IF( K.GT.1 ) THEN 00410 CALL PSGEMV( 'Transpose', K-1, N-K+1, -ONE, Y, IY, 00411 $ JY+K-1, DESCY, A, I, JA, DESCA, DESCA( M_ ), 00412 $ ONE, A, I, J, DESCA, DESCA( M_ ) ) 00413 CALL PSGEMV( 'Transpose', K-1, N-K+1, -ONE, A, IA, J, 00414 $ DESCA, X, IX+K-1, JX, DESCX, DESCX( M_ ), 00415 $ ONE, A, I, J, DESCA, DESCA( M_ ) ) 00416 CALL PSELSET( A, I, J-1, DESCA, ALPHA ) 00417 END IF 00418 * 00419 * Generate reflection P(i) to annihilate A(i,j+1:ja+n-1) 00420 * 00421 CALL PSLARFG( N-K+1, ALPHA, I, J, A, I, J+1, DESCA, 00422 $ DESCA( M_ ), TAUP ) 00423 CALL PSELSET( D, I, 1, DESCD, ALPHA ) 00424 CALL PSELSET( A, I, J, DESCA, ONE ) 00425 * 00426 * Compute X(i+1:ia+m-1,j) 00427 * 00428 CALL PSGEMV( 'No transpose', M-K, N-K+1, ONE, A, I+1, J, 00429 $ DESCA, A, I, J, DESCA, DESCA( M_ ), ZERO, X, 00430 $ IX+K, JX+K-1, DESCX, 1 ) 00431 CALL PSGEMV( 'No transpose', K-1, N-K+1, ONE, Y, IY, JY+K-1, 00432 $ DESCY, A, I, J, DESCA, DESCA( M_ ), ZERO, 00433 $ WORK, IW, 1, DESCW, 1 ) 00434 CALL PSGEMV( 'No transpose', M-K, K-1, -ONE, A, I+1, JA, 00435 $ DESCA, WORK, IW, 1, DESCW, 1, ONE, X, IX+K, 00436 $ JX+K-1, DESCX, 1 ) 00437 CALL PSGEMV( 'No transpose', K-1, N-K+1, ONE, A, IA, J, 00438 $ DESCA, A, I, J, DESCA, DESCA( M_ ), ZERO, 00439 $ WORK, IW, 1, DESCW, 1 ) 00440 CALL PSGEMV( 'No transpose', M-K, K-1, -ONE, X, IX+K, JX, 00441 $ DESCX, WORK, IW, 1, DESCW, 1, ONE, X, IX+K, 00442 $ JX+K-1, DESCX, 1 ) 00443 * 00444 CALL PSELGET( 'Columnwise', ' ', TAU, TAUP, I, 1, DESCTP ) 00445 CALL PSSCAL( M-K, TAU, X, IX+K, JX+K-1, DESCX, 1 ) 00446 * 00447 * Update A(i+1:ia+m-1,j) 00448 * 00449 CALL PSGEMV( 'No transpose', M-K, K-1, -ONE, A, I+1, JA, 00450 $ DESCA, Y, IY, JY+K-1, DESCY, 1, ONE, A, I+1, J, 00451 $ DESCA, 1 ) 00452 CALL PSGEMV( 'No transpose', M-K, K, -ONE, X, IX+K, JX, 00453 $ DESCX, A, IA, J, DESCA, 1, ONE, A, I+1, J, 00454 $ DESCA, 1 ) 00455 CALL PSELSET( A, I, J, DESCA, ALPHA ) 00456 * 00457 * Generate reflection Q(i) to annihilate A(i+2:ia+m-1,j) 00458 * 00459 CALL PSLARFG( M-K, ALPHA, I+1, J, A, MIN( I+2, M+IA-1 ), 00460 $ J, DESCA, 1, TAUQ ) 00461 CALL PSELSET( E, 1, J, DESCE, ALPHA ) 00462 CALL PSELSET( A, I+1, J, DESCA, ONE ) 00463 * 00464 * Compute Y(ia+i:ia+n-1,j) 00465 * 00466 CALL PSGEMV( 'Transpose', M-K, N-K, ONE, A, I+1, J+1, DESCA, 00467 $ A, I+1, J, DESCA, 1, ZERO, WORK( IPY ), 1, 00468 $ JWY, DESCWY, DESCWY( M_ ) ) 00469 CALL PSGEMV( 'Transpose', M-K, K-1, ONE, A, I+1, JA, DESCA, 00470 $ A, I+1, J, DESCA, 1, ZERO, WORK, IW, 1, DESCW, 00471 $ 1 ) 00472 CALL PSGEMV( 'Transpose', K-1, N-K, -ONE, Y, IY, JY+K, 00473 $ DESCY, WORK, IW, 1, DESCW, 1, ONE, WORK( IPY ), 00474 $ 1, JWY, DESCWY, DESCWY( M_ ) ) 00475 CALL PSGEMV( 'Transpose', M-K, K, ONE, X, IX+K, JX, DESCX, 00476 $ A, I+1, J, DESCA, 1, ZERO, WORK, IW, 1, DESCW, 00477 $ 1 ) 00478 CALL PSGEMV( 'Transpose', K, N-K, -ONE, A, IA, J+1, DESCA, 00479 $ WORK, IW, 1, DESCW, 1, ONE, WORK( IPY ), 1, 00480 $ JWY, DESCWY, DESCWY( M_ ) ) 00481 * 00482 CALL PSELGET( 'Rowwise', ' ', TAU, TAUQ, 1, J, DESCTQ ) 00483 CALL PSSCAL( N-K, TAU, WORK( IPY ), 1, JWY, DESCWY, 00484 $ DESCWY( M_ ) ) 00485 CALL PSCOPY( N-K, WORK( IPY ), 1, JWY, DESCWY, DESCWY( M_ ), 00486 $ Y, IY+K-1, JY+K, DESCY, DESCY( M_ ) ) 00487 20 CONTINUE 00488 END IF 00489 * 00490 RETURN 00491 * 00492 * End of PSLABRD 00493 * 00494 END
