SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pdlase2.f
Go to the documentation of this file.
1 SUBROUTINE pdlase2( UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA )
2*
3* -- ScaLAPACK auxiliary routine (version 1.7) --
4* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5* and University of California, Berkeley.
6* May 1, 1997
7*
8* .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER IA, JA, M, N
11 DOUBLE PRECISION ALPHA, BETA
12* ..
13* .. Array Arguments ..
14 INTEGER DESCA( * )
15 DOUBLE PRECISION A( * )
16* ..
17*
18* Purpose
19* =======
20*
21* PDLASE2 initializes an M-by-N distributed matrix sub( A ) denoting
22* A(IA:IA+M-1,JA:JA+N-1) to BETA on the diagonal and ALPHA on the
23* offdiagonals. PDLASE2 requires that only dimension of the matrix
24* operand is distributed.
25*
26* Notes
27* =====
28*
29* Each global data object is described by an associated description
30* vector. This vector stores the information required to establish
31* the mapping between an object element and its corresponding process
32* and memory location.
33*
34* Let A be a generic term for any 2D block cyclicly distributed array.
35* Such a global array has an associated description vector DESCA.
36* In the following comments, the character _ should be read as
37* "of the global array".
38*
39* NOTATION STORED IN EXPLANATION
40* --------------- -------------- --------------------------------------
41* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
42* DTYPE_A = 1.
43* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
44* the BLACS process grid A is distribu-
45* ted over. The context itself is glo-
46* bal, but the handle (the integer
47* value) may vary.
48* M_A (global) DESCA( M_ ) The number of rows in the global
49* array A.
50* N_A (global) DESCA( N_ ) The number of columns in the global
51* array A.
52* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
53* the rows of the array.
54* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
55* the columns of the array.
56* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
57* row of the array A is distributed.
58* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
59* first column of the array A is
60* distributed.
61* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
62* array. LLD_A >= MAX(1,LOCr(M_A)).
63*
64* Let K be the number of rows or columns of a distributed matrix,
65* and assume that its process grid has dimension p x q.
66* LOCr( K ) denotes the number of elements of K that a process
67* would receive if K were distributed over the p processes of its
68* process column.
69* Similarly, LOCc( K ) denotes the number of elements of K that a
70* process would receive if K were distributed over the q processes of
71* its process row.
72* The values of LOCr() and LOCc() may be determined via a call to the
73* ScaLAPACK tool function, NUMROC:
74* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
75* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
76* An upper bound for these quantities may be computed by:
77* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
78* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
79*
80* Arguments
81* =========
82*
83* UPLO (global input) CHARACTER
84* Specifies the part of the distributed matrix sub( A ) to be
85* set:
86* = 'U': Upper triangular part is set; the strictly lower
87* triangular part of sub( A ) is not changed;
88* = 'L': Lower triangular part is set; the strictly upper
89* triangular part of sub( A ) is not changed;
90* Otherwise: All of the matrix sub( A ) is set.
91*
92* M (global input) INTEGER
93* The number of rows to be operated on i.e the number of rows
94* of the distributed submatrix sub( A ). M >= 0.
95*
96* N (global input) INTEGER
97* The number of columns to be operated on i.e the number of
98* columns of the distributed submatrix sub( A ). N >= 0.
99*
100* ALPHA (global input) DOUBLE PRECISION
101* The constant to which the offdiagonal elements are to be
102* set.
103*
104* BETA (global input) DOUBLE PRECISION
105* The constant to which the diagonal elements are to be set.
106*
107* A (local output) DOUBLE PRECISION pointer into the local memory
108* to an array of dimension (LLD_A,LOCc(JA+N-1)). This array
109* contains the local pieces of the distributed matrix sub( A )
110* to be set. On exit, the leading M-by-N submatrix sub( A )
111* is set as follows:
112*
113* if UPLO = 'U', A(IA+i-1,JA+j-1) = ALPHA, 1<=i<=j-1, 1<=j<=N,
114* if UPLO = 'L', A(IA+i-1,JA+j-1) = ALPHA, j+1<=i<=M, 1<=j<=N,
115* otherwise, A(IA+i-1,JA+j-1) = ALPHA, 1<=i<=M, 1<=j<=N,
116* IA+i.NE.JA+j,
117* and, for all UPLO, A(IA+i-1,JA+i-1) = BETA, 1<=i<=min(M,N).
118*
119* IA (global input) INTEGER
120* The row index in the global array A indicating the first
121* row of sub( A ).
122*
123* JA (global input) INTEGER
124* The column index in the global array A indicating the
125* first column of sub( A ).
126*
127* DESCA (global and local input) INTEGER array of dimension DLEN_.
128* The array descriptor for the distributed matrix A.
129*
130* =====================================================================
131*
132* .. Parameters ..
133 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
134 $ LLD_, MB_, M_, NB_, N_, RSRC_
135 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
136 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
137 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
138* ..
139* .. Local Scalars ..
140 INTEGER HEIGHT, IACOL, IAROW, IBASE, ICOFFA, II, IIA,
141 $ IIBEG, IIEND, IINXT, ILEFT, IRIGHT, IROFFA,
142 $ ITOP, JJ, JJA, JJBEG, JJEND, JJNXT, LDA, MBA,
143 $ MP, MPA, MYCOL, MYDIST, MYROW, NBA, NPCOL,
144 $ NPROW, NQ, NQA, WIDE
145* ..
146* .. External Subroutines ..
147 EXTERNAL blacs_gridinfo, dlaset, infog2l
148* ..
149* .. External Functions ..
150 LOGICAL LSAME
151 INTEGER ICEIL, NUMROC
152 EXTERNAL iceil, lsame, numroc
153* ..
154* .. Intrinsic Functions ..
155 INTRINSIC max, min, mod
156* ..
157* .. Executable Statements ..
158*
159 IF( m.EQ.0 .OR. n.EQ.0 )
160 $ RETURN
161*
162* Get grid parameters
163*
164 CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
165*
166 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
167 $ iarow, iacol )
168 mba = desca( mb_ )
169 nba = desca( nb_ )
170 lda = desca( lld_ )
171 iroffa = mod( ia-1, mba )
172 icoffa = mod( ja-1, nba )
173*
174 IF( n.LE.( nba-icoffa ) ) THEN
175*
176* It is assumed that the local columns JJA:JJA+N-1 of the matrix
177* A are in the same process column (IACOL).
178*
179* N
180* JJA JJA+N-1
181* / --------------------- \
182* IROFFA| | | |
183* \ |...................| | ( IAROW )
184* IIA |x | | MB_A
185* | x | |
186* |--x----------------| /
187* | x |
188* | x | ITOP
189* | x | |
190* | x | /-------\
191* |-------x-----------| |-------x-----------|
192* | x | | x |
193* | x | | x |
194* | x | | x |
195* | x | | x |
196* |------------x------| |------------x------|
197* | x | \____________/
198* | x | |
199* | x | IBASE
200* | x |
201* |-----------------x-| Local picture
202* | x|
203* | |
204* | |
205* | |
206* |-------------------|
207* | |
208* . .
209* . .
210* . (IACOL) .
211*
212 IF( mycol.EQ.iacol ) THEN
213*
214 mpa = numroc( m+iroffa, mba, myrow, iarow, nprow )
215 IF( mpa.LE.0 )
216 $ RETURN
217 IF( myrow.EQ.iarow )
218 $ mpa = mpa - iroffa
219 mydist = mod( myrow-iarow+nprow, nprow )
220 itop = mydist * mba - iroffa
221*
222 IF( lsame( uplo, 'U' ) ) THEN
223*
224 itop = max( 0, itop )
225 iibeg = iia
226 iiend = iia + mpa - 1
227 iinxt = min( iceil( iibeg, mba ) * mba, iiend )
228*
229 10 CONTINUE
230 IF( ( n-itop ).GT.0 ) THEN
231 CALL dlaset( uplo, iinxt-iibeg+1, n-itop, alpha, beta,
232 $ a( iibeg+(jja+itop-1)*lda ), lda )
233 mydist = mydist + nprow
234 itop = mydist * mba - iroffa
235 iibeg = iinxt +1
236 iinxt = min( iinxt+mba, iiend )
237 GO TO 10
238 END IF
239*
240 ELSE IF( lsame( uplo, 'L' ) ) THEN
241*
242 ii = iia
243 jj = jja
244 mp = mpa
245 ibase = min( itop+mba, n )
246 itop = min( max( 0, itop ), n )
247*
248 20 CONTINUE
249 IF( jj.LE.( jja+n-1 ) ) THEN
250 height = ibase - itop
251 CALL dlaset( 'All', mp, itop-jj+jja, alpha, alpha,
252 $ a( ii+(jj-1)*lda ), lda )
253 CALL dlaset( uplo, mp, height, alpha, beta,
254 $ a( ii+(jja+itop-1)*lda ), lda )
255 mp = max( 0, mp - height )
256 ii = ii + height
257 jj = jja + ibase
258 mydist = mydist + nprow
259 itop = mydist * mba - iroffa
260 ibase = min( itop + mba, n )
261 itop = min( itop, n )
262 GO TO 20
263 END IF
264*
265 ELSE
266*
267 ii = iia
268 jj = jja
269 mp = mpa
270 ibase = min( itop+mba, n )
271 itop = min( max( 0, itop ), n )
272*
273 30 CONTINUE
274 IF( jj.LE.( jja+n-1 ) ) THEN
275 height = ibase - itop
276 CALL dlaset( 'All', mpa, itop-jj+jja, alpha, alpha,
277 $ a( iia+(jj-1)*lda ), lda )
278 CALL dlaset( 'All', mpa-mp, height, alpha, alpha,
279 $ a( iia+(jja+itop-1)*lda ), lda )
280 CALL dlaset( 'All', mp, height, alpha, beta,
281 $ a( ii+(jja+itop-1)*lda ), lda )
282 mp = max( 0, mp - height )
283 ii = ii + height
284 jj = jja + ibase
285 mydist = mydist + nprow
286 itop = mydist * mba - iroffa
287 ibase = min( itop + mba, n )
288 itop = min( itop, n )
289 GO TO 30
290 END IF
291*
292 END IF
293*
294 END IF
295*
296 ELSE IF( m.LE.( mba-iroffa ) ) THEN
297*
298* It is assumed that the local rows IIA:IIA+M-1 of the matrix A
299* are in the same process row (IAROW).
300*
301* ICOFFA
302* / \JJA
303* IIA ------------------ .... --------
304* | .x | | | / | | \
305* | . x | | | ILEFT| | | |
306* | . x | | | | | |
307* | . x | | \ x | |
308* | . |x | | |x | | IRIGHT
309* | . | x | | | x | |
310* (IAROW) | . | x | | | x | |
311* | . | x| | | x| |
312* | . | x | | x /
313* | . | |x | | |
314* | . | | x | | |
315* | . | | x | | |
316* | . | | x| | |
317* IIA+M-1 ------------------ .... -------
318* NB_A
319* (IACOL) Local picture
320*
321 IF( myrow.EQ.iarow ) THEN
322*
323 nqa = numroc( n+icoffa, nba, mycol, iacol, npcol )
324 IF( nqa.LE.0 )
325 $ RETURN
326 IF( mycol.EQ.iacol )
327 $ nqa = nqa - icoffa
328 mydist = mod( mycol-iacol+npcol, npcol )
329 ileft = mydist * nba - icoffa
330*
331 IF( lsame( uplo, 'L' ) ) THEN
332*
333 ileft = max( 0, ileft )
334 jjbeg = jja
335 jjend = jja + nqa - 1
336 jjnxt = min( iceil( jjbeg, nba ) * nba, jjend )
337*
338 40 CONTINUE
339 IF( ( m-ileft ).GT.0 ) THEN
340 CALL dlaset( uplo, m-ileft, jjnxt-jjbeg+1, alpha,
341 $ beta, a( iia+ileft+(jjbeg-1)*lda ), lda )
342 mydist = mydist + npcol
343 ileft = mydist * nba - icoffa
344 jjbeg = jjnxt +1
345 jjnxt = min( jjnxt+nba, jjend )
346 GO TO 40
347 END IF
348*
349 ELSE IF( lsame( uplo, 'U' ) ) THEN
350*
351 ii = iia
352 jj = jja
353 nq = nqa
354 iright = min( ileft+nba, m )
355 ileft = min( max( 0, ileft ), m )
356*
357 50 CONTINUE
358 IF( ii.LE.( iia+m-1 ) ) THEN
359 wide = iright - ileft
360 CALL dlaset( 'All', ileft-ii+iia, nq, alpha, alpha,
361 $ a( ii+(jj-1)*lda ), lda )
362 CALL dlaset( uplo, wide, nq, alpha, beta,
363 $ a( iia+ileft+(jj-1)*lda ), lda )
364 nq = max( 0, nq - wide )
365 ii = iia + iright
366 jj = jj + wide
367 mydist = mydist + npcol
368 ileft = mydist * nba - icoffa
369 iright = min( ileft + nba, m )
370 ileft = min( ileft, m )
371 GO TO 50
372 END IF
373*
374 ELSE
375*
376 ii = iia
377 jj = jja
378 nq = nqa
379 iright = min( ileft+nba, m )
380 ileft = min( max( 0, ileft ), m )
381*
382 60 CONTINUE
383 IF( ii.LE.( iia+m-1 ) ) THEN
384 wide = iright - ileft
385 CALL dlaset( 'All', ileft-ii+iia, nqa, alpha, alpha,
386 $ a( ii+(jja-1)*lda ), lda )
387 CALL dlaset( 'All', wide, nqa-nq, alpha, alpha,
388 $ a( iia+ileft+(jja-1)*lda ), lda )
389 CALL dlaset( 'All', wide, nq, alpha, beta,
390 $ a( iia+ileft+(jj-1)*lda ), lda )
391 nq = max( 0, nq - wide )
392 ii = iia + iright
393 jj = jj + wide
394 mydist = mydist + npcol
395 ileft = mydist * nba - icoffa
396 iright = min( ileft + nba, m )
397 ileft = min( ileft, m )
398 GO TO 60
399 END IF
400*
401 END IF
402*
403 END IF
404*
405 END IF
406*
407 RETURN
408*
409* End of PDLASE2
410*
411 END
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
#define max(A, B)
Definition pcgemr.c:180
#define min(A, B)
Definition pcgemr.c:181
subroutine pdlase2(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pdlase2.f:2