OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pclarft.f
Go to the documentation of this file.
1 SUBROUTINE pclarft( DIRECT, STOREV, N, K, V, IV, JV, DESCV, TAU,
2 $ T, 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 CHARACTER DIRECT, STOREV
11 INTEGER IV, JV, K, N
12* ..
13* .. Array Arguments ..
14 INTEGER DESCV( * )
15 COMPLEX TAU( * ), T( * ), V( * ), WORK( * )
16* ..
17*
18* Purpose
19* =======
20*
21* PCLARFT forms the triangular factor T of a complex block reflector H
22* of order n, which is defined as a product of k elementary reflectors.
23*
24* If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
25*
26* If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
27*
28* If STOREV = 'C', the vector which defines the elementary reflector
29* H(i) is stored in the i-th column of the distributed matrix V, and
30*
31* H = I - V * T * V'
32*
33* If STOREV = 'R', the vector which defines the elementary reflector
34* H(i) is stored in the i-th row of the distributed matrix V, and
35*
36* H = I - V' * T * V
37*
38* Notes
39* =====
40*
41* Each global data object is described by an associated description
42* vector. This vector stores the information required to establish
43* the mapping between an object element and its corresponding process
44* and memory location.
45*
46* Let A be a generic term for any 2D block cyclicly distributed array.
47* Such a global array has an associated description vector DESCA.
48* In the following comments, the character _ should be read as
49* "of the global array".
50*
51* NOTATION STORED IN EXPLANATION
52* --------------- -------------- --------------------------------------
53* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
54* DTYPE_A = 1.
55* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
56* the BLACS process grid A is distribu-
57* ted over. The context itself is glo-
58* bal, but the handle (the integer
59* value) may vary.
60* M_A (global) DESCA( M_ ) The number of rows in the global
61* array A.
62* N_A (global) DESCA( N_ ) The number of columns in the global
63* array A.
64* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
65* the rows of the array.
66* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
67* the columns of the array.
68* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
69* row of the array A is distributed.
70* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
71* first column of the array A is
72* distributed.
73* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
74* array. LLD_A >= MAX(1,LOCr(M_A)).
75*
76* Let K be the number of rows or columns of a distributed matrix,
77* and assume that its process grid has dimension p x q.
78* LOCr( K ) denotes the number of elements of K that a process
79* would receive if K were distributed over the p processes of its
80* process column.
81* Similarly, LOCc( K ) denotes the number of elements of K that a
82* process would receive if K were distributed over the q processes of
83* its process row.
84* The values of LOCr() and LOCc() may be determined via a call to the
85* ScaLAPACK tool function, NUMROC:
86* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
87* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
88* An upper bound for these quantities may be computed by:
89* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
90* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
91*
92* Arguments
93* =========
94*
95* DIRECT (global input) CHARACTER*1
96* Specifies the order in which the elementary reflectors are
97* multiplied to form the block reflector:
98* = 'F': H = H(1) H(2) . . . H(k) (Forward)
99* = 'B': H = H(k) . . . H(2) H(1) (Backward)
100*
101* STOREV (global input) CHARACTER*1
102* Specifies how the vectors which define the elementary
103* reflectors are stored (see also Further Details):
104* = 'C': columnwise
105* = 'R': rowwise
106*
107* N (global input) INTEGER
108* The order of the block reflector H. N >= 0.
109*
110* K (global input) INTEGER
111* The order of the triangular factor T (= the number of
112* elementary reflectors). 1 <= K <= MB_V (= NB_V).
113*
114* V (input/output) COMPLEX pointer into the local memory
115* to an array of local dimension (LOCr(IV+N-1),LOCc(JV+K-1))
116* if STOREV = 'C', and (LOCr(IV+K-1),LOCc(JV+N-1)) if
117* STOREV = 'R'. The distributed matrix V contains the
118* Householder vectors. See further details.
119*
120* IV (global input) INTEGER
121* The row index in the global array V indicating the first
122* row of sub( V ).
123*
124* JV (global input) INTEGER
125* The column index in the global array V indicating the
126* first column of sub( V ).
127*
128* DESCV (global and local input) INTEGER array of dimension DLEN_.
129* The array descriptor for the distributed matrix V.
130*
131* TAU (local input) COMPLEX, array, dimension LOCr(IV+K-1)
132* if INCV = M_V, and LOCc(JV+K-1) otherwise. This array
133* contains the Householder scalars related to the Householder
134* vectors. TAU is tied to the distributed matrix V.
135*
136* T (local output) COMPLEX array, dimension (NB_V,NB_V)
137* if STOREV = 'Col', and (MB_V,MB_V) otherwise. It contains
138* the k-by-k triangular factor of the block reflector asso-
139* ciated with V. If DIRECT = 'F', T is upper triangular;
140* if DIRECT = 'B', T is lower triangular.
141*
142* WORK (local workspace) COMPLEX array,
143* dimension (K*(K-1)/2)
144*
145* Further Details
146* ===============
147*
148* The shape of the matrix V and the storage of the vectors which define
149* the H(i) is best illustrated by the following example with n = 5 and
150* k = 3. The elements equal to 1 are not stored; the corresponding
151* array elements are modified but restored on exit. The rest of the
152* array is not used.
153*
154* DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
155*
156* V( IV:IV+N-1, ( 1 ) V( IV:IV+K-1, ( 1 v1 v1 v1 v1 )
157* JV:JV+K-1 ) = ( v1 1 ) JV:JV+N-1 ) = ( 1 v2 v2 v2 )
158* ( v1 v2 1 ) ( 1 v3 v3 )
159* ( v1 v2 v3 )
160* ( v1 v2 v3 )
161*
162* DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
163*
164* V( IV:IV+N-1, ( v1 v2 v3 ) V( IV:IV+K-1, ( v1 v1 1 )
165* JV:JV+K-1 ) = ( v1 v2 v3 ) JV:JV+N-1 ) = ( v2 v2 v2 1 )
166* ( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
167* ( 1 v3 )
168* ( 1 )
169*
170* =====================================================================
171*
172* .. Parameters ..
173 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
174 $ lld_, mb_, m_, nb_, n_, rsrc_
175 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
176 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
177 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
178 COMPLEX ONE, ZERO
179 parameter( one = ( 1.0e+0, 0.0e+0 ),
180 $ zero = ( 0.0e+0, 0.0e+0 ) )
181* ..
182* .. Local Scalars ..
183 LOGICAL FORWARD
184 INTEGER ICOFF, ICTXT, II, IIV, IROFF, IVCOL, IVROW,
185 $ itmp0, itmp1, iw, jj, jjv, ldv, micol, mirow,
186 $ mycol, myrow, np, npcol, nprow, nq
187 COMPLEX VII
188* ..
189* .. External Subroutines ..
190 EXTERNAL blacs_gridinfo, ccopy, cgemv, cgsum2d,
191 $ clacgv, claset, ctrmv, infog2l
192* ..
193* .. External Functions ..
194 LOGICAL LSAME
195 INTEGER INDXG2P, NUMROC
196 EXTERNAL indxg2p, lsame, numroc
197* ..
198* .. Intrinsic Functions ..
199 INTRINSIC mod
200* ..
201* .. Executable Statements ..
202*
203* Quick return if possible
204*
205 IF( n.LE.0 .OR. k.LE.0 )
206 $ RETURN
207*
208 ictxt = descv( ctxt_ )
209 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
210*
211 forward = lsame( direct, 'F' )
212 CALL infog2l( iv, jv, descv, nprow, npcol, myrow, mycol,
213 $ iiv, jjv, ivrow, ivcol )
214*
215 IF( lsame( storev, 'C' ) .AND. mycol.EQ.ivcol ) THEN
216*
217 iw = 1
218 ldv = descv( lld_ )
219 iroff = mod( iv-1, descv( mb_ ) )
220*
221 IF( forward ) THEN
222*
223* DIRECT = 'Forward', STOREV = 'Columnwise'
224*
225 np = numroc( n+iroff, descv( mb_ ), myrow, ivrow, nprow )
226 IF( myrow.EQ.ivrow ) THEN
227 np = np - iroff
228 ii = iiv + 1
229 ELSE
230 ii = iiv
231 END IF
232 IF( iroff+1.EQ.descv( mb_ ) ) THEN
233 mirow = mod( ivrow+1, nprow )
234 ELSE
235 mirow = ivrow
236 END IF
237 itmp0 = 0
238*
239 DO 10 jj = jjv+1, jjv+k-1
240*
241 IF( myrow.EQ.mirow ) THEN
242 vii = v( ii+(jj-1)*ldv )
243 v( ii+(jj-1)*ldv ) = one
244 END IF
245*
246* T(1:i-1,i) = -tau( jv+i-1 ) *
247* V(iv+i-1:iv+n-1,jv:jv+i-2)' * V(iv+i-1:iv+n-1,jv+i-1)
248*
249 itmp0 = itmp0 + 1
250 IF( np-ii+iiv.GT.0 ) THEN
251 CALL cgemv( 'Conjugate transpose', np-ii+iiv, itmp0,
252 $ -tau( jj ), v( ii+(jjv-1)*ldv ), ldv,
253 $ v( ii+(jj-1)*ldv ), 1, zero,
254 $ work( iw ), 1 )
255 ELSE
256 CALL claset( 'All', itmp0, 1, zero, zero, work( iw ),
257 $ itmp0 )
258 END IF
259*
260 iw = iw + itmp0
261 IF( myrow.EQ.mirow ) THEN
262 v( ii+(jj-1)*ldv ) = vii
263 ii = ii + 1
264 END IF
265*
266 IF( mod( iv+itmp0, descv( mb_ ) ).EQ.0 )
267 $ mirow = mod( mirow+1, nprow )
268*
269 10 CONTINUE
270*
271 CALL cgsum2d( ictxt, 'Columnwise', ' ', iw-1, 1, work, iw-1,
272 $ ivrow, mycol )
273*
274 IF( myrow.EQ.ivrow ) THEN
275*
276 iw = 1
277 itmp0 = 0
278 itmp1 = 1
279*
280 t( itmp1 ) = tau( jjv )
281*
282 DO 20 jj = jjv+1, jjv+k-1
283*
284* T(1:j-1,j) = T(1:j-1,1:j-1) * T(1:j-1,j)
285*
286 itmp0 = itmp0 + 1
287 itmp1 = itmp1 + descv( nb_ )
288 CALL ccopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
289 iw = iw + itmp0
290*
291 CALL ctrmv( 'upper', 'no transpose', 'non-unit',
292 $ ITMP0, T, DESCV( NB_ ), T( ITMP1 ), 1 )
293 T(ITMP1+ITMP0) = TAU( JJ )
294*
295 20 CONTINUE
296*
297 END IF
298*
299 ELSE
300*
301* DIRECT = 'Backward', STOREV = 'Columnwise'
302*
303 NP = NUMROC( N+IROFF-1, DESCV( MB_ ), MYROW, IVROW, NPROW )
304.EQ. IF( MYROWIVROW )
305 $ NP = NP - IROFF
306 MIROW = INDXG2P( IV+N-2, DESCV( MB_ ), MYROW,
307 $ DESCV( RSRC_ ), NPROW )
308 II = IIV + NP - 1
309 ITMP0 = 0
310*
311 DO 30 JJ = JJV+K-2, JJV, -1
312*
313.EQ. IF( MYROWMIROW ) THEN
314 VII = V( II+(JJ-1)*LDV )
315 V( II+(JJ-1)*LDV ) = ONE
316 END IF
317*
318* T(1:i-1,i) = -tau( jv+i-1 ) *
319* V(iv:iv+n-k+i-1,jv+i:jv+k-1)' * V(iv:iv+n-k+i-1,jv+i-1)
320*
321 ITMP0 = ITMP0 + 1
322.GT. IF( II-IIV+10 ) THEN
323 CALL CGEMV( 'conjugate transpose', II-IIV+1, ITMP0,
324 $ -TAU( JJ ), V( IIV+JJ*LDV ), LDV,
325 $ V( IIV+(JJ-1)*LDV ), 1, ZERO,
326 $ WORK( IW ), 1 )
327 ELSE
328 CALL CLASET( 'all', ITMP0, 1, ZERO, ZERO, WORK( IW ),
329 $ ITMP0 )
330 END IF
331*
332 IW = IW + ITMP0
333.EQ. IF( MYROWMIROW ) THEN
334 V( II+(JJ-1)*LDV ) = VII
335 II = II - 1
336 END IF
337*
338.EQ. IF( MOD( IV+N-ITMP0-2, DESCV(MB_) )0 )
339 $ MIROW = MOD( MIROW+NPROW-1, NPROW )
340*
341 30 CONTINUE
342*
343 CALL CGSUM2D( ICTXT, 'columnwise', ' ', IW-1, 1, WORK, IW-1,
344 $ IVROW, MYCOL )
345*
346.EQ. IF( MYROWIVROW ) THEN
347*
348 IW = 1
349 ITMP0 = 0
350 ITMP1 = K + 1 + (K-1) * DESCV( NB_ )
351*
352 T( ITMP1-1 ) = TAU( JJV+K-1 )
353*
354 DO 40 JJ = JJV+K-2, JJV, -1
355*
356* T(j+1:k,j) = T(j+1:k,j+1:k) * T(j+1:k,j)
357*
358 ITMP0 = ITMP0 + 1
359 ITMP1 = ITMP1 - DESCV( NB_ ) - 1
360 CALL CCOPY( ITMP0, WORK( IW ), 1, T( ITMP1 ), 1 )
361 IW = IW + ITMP0
362*
363 CALL CTRMV( 'lower', 'no transpose', 'non-unit',
364 $ ITMP0, T( ITMP1+DESCV( NB_ ) ),
365 $ DESCV( NB_ ), T( ITMP1 ), 1 )
366 T( ITMP1-1 ) = TAU( JJ )
367*
368 40 CONTINUE
369*
370 END IF
371*
372 END IF
373*
374 ELSE IF( LSAME( STOREV, 'r.AND..EQ.' ) MYROWIVROW ) THEN
375*
376 IW = 1
377 LDV = DESCV( LLD_ )
378 ICOFF = MOD( JV-1, DESCV( NB_ ) )
379*
380 IF( FORWARD ) THEN
381*
382* DIRECT = 'Forward', STOREV = 'Rowwise'
383*
384 NQ = NUMROC( N+ICOFF, DESCV( NB_ ), MYCOL, IVCOL, NPCOL )
385.EQ. IF( MYCOLIVCOL ) THEN
386 NQ = NQ - ICOFF
387 JJ = JJV + 1
388 ELSE
389 JJ = JJV
390 END IF
391.EQ. IF( ICOFF+1DESCV( NB_ ) ) THEN
392 MICOL = MOD( IVCOL+1, NPCOL )
393 ELSE
394 MICOL = IVCOL
395 END IF
396 ITMP0 = 0
397*
398 DO 50 II = IIV+1, IIV+K-1
399*
400.EQ. IF( MYCOLMICOL ) THEN
401 VII = V( II+(JJ-1)*LDV )
402 V( II+(JJ-1)*LDV ) = ONE
403 END IF
404*
405* T(1:i-1,i) = -tau( iv+i-1 ) *
406* V(iv+i-1,jv+i-1:jv+n-1) * V(iv:iv+i-2,jv+i-1:jv+n-1)'
407*
408 ITMP0 = ITMP0 + 1
409.GT. IF( NQ-JJ+JJV0 ) THEN
410 CALL CLACGV( NQ-JJ+JJV, V( II+(JJ-1)*LDV ), LDV )
411 CALL CGEMV( 'no transpose', ITMP0, NQ-JJ+JJV,
412 $ -TAU(II), V( IIV+(JJ-1)*LDV ), LDV,
413 $ V( II+(JJ-1)*LDV ), LDV, ZERO,
414 $ WORK( IW ), 1 )
415 CALL CLACGV( NQ-JJ+JJV, V( II+(JJ-1)*LDV ), LDV )
416 ELSE
417 CALL CLASET( 'all', ITMP0, 1, ZERO, ZERO,
418 $ WORK( IW ), ITMP0 )
419 END IF
420*
421 IW = IW + ITMP0
422.EQ. IF( MYCOLMICOL ) THEN
423 V( II+(JJ-1)*LDV ) = VII
424 JJ = JJ + 1
425 END IF
426*
427.EQ. IF( MOD( JV+ITMP0, DESCV( NB_ ) )0 )
428 $ MICOL = MOD( MICOL+1, NPCOL )
429*
430 50 CONTINUE
431*
432 CALL CGSUM2D( ICTXT, 'rowwise', ' ', IW-1, 1, WORK, IW-1,
433 $ MYROW, IVCOL )
434*
435.EQ. IF( MYCOLIVCOL ) THEN
436*
437 IW = 1
438 ITMP0 = 0
439 ITMP1 = 1
440*
441 T( ITMP1 ) = TAU( IIV )
442*
443 DO 60 II = IIV+1, IIV+K-1
444*
445* T(1:i-1,i) = T(1:i-1,1:i-1) * T(1:i-1,i)
446*
447 ITMP0 = ITMP0 + 1
448 ITMP1 = ITMP1 + DESCV( MB_ )
449 CALL CCOPY( ITMP0, WORK( IW ), 1, T( ITMP1 ), 1 )
450 IW = IW + ITMP0
451*
452 CALL CTRMV( 'upper', 'no transpose', 'non-unit',
453 $ ITMP0, T, DESCV( MB_ ), T( ITMP1 ), 1 )
454 T( ITMP1+ITMP0 ) = TAU( II )
455*
456 60 CONTINUE
457*
458 END IF
459*
460 ELSE
461*
462* DIRECT = 'Backward', STOREV = 'Rowwise'
463*
464 NQ = NUMROC( N+ICOFF-1, DESCV( NB_ ), MYCOL, IVCOL, NPCOL )
465.EQ. IF( MYCOLIVCOL )
466 $ NQ = NQ - ICOFF
467 MICOL = INDXG2P( JV+N-2, DESCV( NB_ ), MYCOL,
468 $ DESCV( CSRC_ ), NPCOL )
469 JJ = JJV + NQ - 1
470 ITMP0 = 0
471*
472 DO 70 II = IIV+K-2, IIV, -1
473*
474.EQ. IF( MYCOLMICOL ) THEN
475 VII = V( II+(JJ-1)*LDV )
476 V( II+(JJ-1)*LDV ) = ONE
477 END IF
478*
479* T(i+1:k,i) = -tau( iv+i-1 ) *
480* V(iv+i:iv+k-1,jv:jv+n-k+i-1)' * V(iv+i-1,jv:jv+n-k+i-1)'
481*
482 ITMP0 = ITMP0 + 1
483.GT. IF( JJ-JJV+10 ) THEN
484 CALL CLACGV( JJ-JJV+1, V( II+(JJV-1)*LDV ), LDV )
485 CALL CGEMV( 'no transpose', itmp0, jj-jjv+1,
486 $ -tau( ii ), v( ii+1+(jjv-1)*ldv ), ldv,
487 $ v( ii+(jjv-1)*ldv ), ldv, zero,
488 $ work( iw ), 1 )
489 CALL clacgv( jj-jjv+1, v( ii+(jjv-1)*ldv ), ldv )
490 ELSE
491 CALL claset( 'All', itmp0, 1, zero, zero,
492 $ work( iw ), itmp0 )
493 END IF
494*
495 iw = iw + itmp0
496 IF( mycol.EQ.micol ) THEN
497 v( ii+(jj-1)*ldv ) = vii
498 jj = jj - 1
499 END IF
500*
501 IF( mod( jv+n-itmp0-2, descv( nb_ ) ).EQ.0 )
502 $ micol = mod( micol+npcol-1, npcol )
503*
504 70 CONTINUE
505*
506 CALL cgsum2d( ictxt, 'Rowwise', ' ', iw-1, 1, work, iw-1,
507 $ myrow, ivcol )
508*
509 IF( mycol.EQ.ivcol ) THEN
510*
511 iw = 1
512 itmp0 = 0
513 itmp1 = k + 1 + (k-1) * descv( mb_ )
514*
515 t( itmp1-1 ) = tau( iiv+k-1 )
516*
517 DO 80 ii = iiv+k-2, iiv, -1
518*
519* T(i+1:k,i) = T(i+1:k,i+1:k) * T(i+1:k,i)
520*
521 itmp0 = itmp0 + 1
522 itmp1 = itmp1 - descv( mb_ ) - 1
523 CALL ccopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
524 iw = iw + itmp0
525*
526 CALL ctrmv( 'Lower', 'No transpose', 'Non-unit',
527 $ itmp0, t( itmp1+descv( mb_ ) ),
528 $ descv( mb_ ), t( itmp1 ), 1 )
529 t( itmp1-1 ) = tau( ii )
530*
531 80 CONTINUE
532*
533 END IF
534*
535 END IF
536*
537 END IF
538*
539 RETURN
540*
541* End of PCLARFT
542*
543 END
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
ctrmv
subroutine ctrmv(uplo, trans, diag, n, a, lda, x, incx)
CTRMV
Definition ctrmv.f:147
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158
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
pclarft
subroutine pclarft(direct, storev, n, k, v, iv, jv, descv, tau, t, work)
Definition pclarft.f:3