OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pdlabrd.f
Go to the documentation of this file.
1 SUBROUTINE pdlabrd( 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 DOUBLE PRECISION A( * ), D( * ), E( * ), TAUP( * ),
15 $ tauq( * ), x( * ), y( * ), work( * )
16* ..
17*
18* Purpose
19* =======
20*
21* PDLABRD 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 PDGEBRD.
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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DOUBLE PRECISION ONE, ZERO
254 parameter( one = 1.0d+0, zero = 0.0d+0 )
255* ..
256* .. Local Scalars ..
257 INTEGER I, IACOL, IAROW, ICTXT, II, IPY, IW, J, JJ,
258 $ jwy, k, mycol, myrow, npcol, nprow
259 DOUBLE PRECISION 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, pdcopy,
266 $ pdelget, pdelset, pdgemv, pdlarfg,
267 $ pdscal
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 pdgemv( '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 pdgemv( '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 pdelset( 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 pdlarfg( m-k+1, alpha, i, j, a, i+1, j, desca, 1,
326 $ tauq )
327 CALL pdelset( d, 1, j, descd, alpha )
328 CALL pdelset( a, i, j, desca, one )
329*
330* Compute Y(IA+I:IA+N-1,J)
331*
332 CALL pdgemv( '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 pdgemv( '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 pdgemv( '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 pdgemv( '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 pdgemv( '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 pdelget( 'Rowwise', ' ', TAU, TAUQ, 1, J, DESCTQ )
349 CALL PDSCAL( N-K, TAU, WORK( IPY ), 1, JWY, DESCWY,
350 $ DESCWY( M_ ) )
351 CALL PDCOPY( 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 PDGEMV( '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 PDGEMV( '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 PDELSET( A, I, J, DESCA, ALPHA )
363*
364* Generate reflection P(i) to annihilate A(i,j+2:ja+n-1)
365*
366 CALL PDLARFG( N-K, ALPHA, I, J+1, A, I,
367 $ MIN( J+2, N+JA-1 ), DESCA, DESCA( M_ ), TAUP )
368 CALL PDELSET( E, I, 1, DESCE, ALPHA )
369 CALL PDELSET( A, I, J+1, DESCA, ONE )
370*
371* Compute X(I+1:IA+M-1,J)
372*
373 CALL PDGEMV( '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 PDGEMV( '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 PDGEMV( '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 PDGEMV( '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 PDGEMV( '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 PDELGET( 'columnwise', ' ', TAU, TAUP, I, 1, DESCTP )
390 CALL PDSCAL( 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.GT. IF( K1 ) THEN
410 CALL PDGEMV( '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 PDGEMV( '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 PDELSET( 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 PDLARFG( N-K+1, ALPHA, I, J, A, I, J+1, DESCA,
422 $ DESCA( M_ ), TAUP )
423 CALL PDELSET( D, I, 1, DESCD, ALPHA )
424 CALL PDELSET( A, I, J, DESCA, ONE )
425*
426* Compute X(i+1:ia+m-1,j)
427*
428 CALL PDGEMV( '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 PDGEMV( '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 PDGEMV( '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 PDGEMV( '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 PDGEMV( '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 PDELGET( 'columnwise', ' ', TAU, TAUP, I, 1, DESCTP )
445 CALL PDSCAL( M-K, TAU, X, IX+K, JX+K-1, DESCX, 1 )
446*
447* Update A(i+1:ia+m-1,j)
448*
449 CALL PDGEMV( '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 PDGEMV( '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 PDELSET( A, I, J, DESCA, ALPHA )
456*
457* Generate reflection Q(i) to annihilate A(i+2:ia+m-1,j)
458*
459 CALL PDLARFG( M-K, ALPHA, I+1, J, A, MIN( I+2, M+IA-1 ),
460 $ J, DESCA, 1, TAUQ )
461 CALL PDELSET( E, 1, J, DESCE, ALPHA )
462 CALL PDELSET( A, I+1, J, DESCA, ONE )
463*
464* Compute Y(ia+i:ia+n-1,j)
465*
466 CALL PDGEMV( '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 PDGEMV( '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 PDGEMV( '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 PDGEMV( '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 PDGEMV( '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 PDELGET( 'rowwise', ' ', TAU, TAUQ, 1, J, DESCTQ )
483 CALL PDSCAL( N-K, TAU, WORK( IPY ), 1, JWY, DESCWY,
484 $ DESCWY( M_ ) )
485 CALL PDCOPY( 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 PDLABRD
493*
494 END
min
#define min(a, b)
Definition macros.h:20
pdscal
subroutine pdscal(n, alpha, x, ix, jx, descx, incx)
Definition mpi.f:978
descset
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
Definition mpi.f:1610
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition mpi.f:937
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
pdelget
subroutine pdelget(scope, top, alpha, a, ia, ja, desca)
Definition pdelget.f:2
pdelset
subroutine pdelset(a, ia, ja, desca, alpha)
Definition pdelset.f:2
pdlabrd
subroutine pdlabrd(m, n, nb, a, ia, ja, desca, d, e, tauq, taup, x, ix, jx, descx, y, iy, jy, descy, work)
Definition pdlabrd.f:3
pdlarfg
subroutine pdlarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)
Definition pdlarfg.f:3