|
|
| |
| # lines: |
795 | | # code: |
795 | | # comment: | 0 | |
# blank: | 0 |
| # Variables: | 60 |
| # Callers: | 1 |
| # Callings: | 1 |
| # Words: | 231 |
| # Keywords: | 137 |
|
|
|
|
|
..
.. Array Arguments ..
..
Purpose
=======
PCDTTRS solves a system of linear equations
A(1:N, JA:JA+N-1) * X = B(IB:IB+N-1, 1:NRHS)
or
A(1:N, JA:JA+N-1)' * X = B(IB:IB+N-1, 1:NRHS)
where A(1:N, JA:JA+N-1) is the matrix used to produce the factors
stored in A(1:N,JA:JA+N-1) and AF by PCDTTRF.
A(1:N, JA:JA+N-1) is an N-by-N complex
tridiagonal diagonally dominant-like distributed
matrix.
Routine PCDTTRF MUST be called first.
=====================================================================
Arguments
=========
TRANS (global input) CHARACTER
= 'N': Solve with A(1:N, JA:JA+N-1);
= 'C': Solve with conjugate_transpose( A(1:N, JA:JA+N-1) );
N (global input) INTEGER
The number of rows and columns to be operated on, i.e. the
order of the distributed submatrix A(1:N, JA:JA+N-1). N >= 0.
NRHS (global input) INTEGER
The number of right hand sides, i.e., the number of columns
of the distributed submatrix B(IB:IB+N-1, 1:NRHS).
NRHS >= 0.
DL (local input/local output) COMPLEX pointer to local
part of global vector storing the lower diagonal of the
matrix. Globally, DL(1) is not referenced, and DL must be
aligned with D.
Must be of size >= DESCA( NB_ ).
On exit, this array contains information containing the
factors of the matrix.
D (local input/local output) COMPLEX pointer to local
part of global vector storing the main diagonal of the
matrix.
On exit, this array contains information containing the
factors of the matrix.
Must be of size >= DESCA( NB_ ).
DU (local input/local output) COMPLEX pointer to local
part of global vector storing the upper diagonal of the
matrix. Globally, DU(n) is not referenced, and DU must be
aligned with D.
On exit, this array contains information containing the
factors of the matrix.
Must be of size >= DESCA( NB_ ).
JA (global input) INTEGER
The index in the global array A that points to the start of
the matrix to be operated on (which may be either all of A
or a submatrix of A).
DESCA (global and local input) INTEGER array of dimension DLEN.
if 1D type (DTYPE_A=501 or 502), DLEN >= 7;
if 2D type (DTYPE_A=1), DLEN >= 9.
The array descriptor for the distributed matrix A.
Contains information of mapping of A to memory. Please
see NOTES below for full description and options.
B (local input/local output) COMPLEX pointer into
local memory to an array of local lead dimension lld_b>=NB.
On entry, this array contains the
the local pieces of the right hand sides
B(IB:IB+N-1, 1:NRHS).
On exit, this contains the local piece of the solutions
distributed matrix X.
IB (global input) INTEGER
The row index in the global array B that points to the first
row of the matrix to be operated on (which may be either
all of B or a submatrix of B).
DESCB (global and local input) INTEGER array of dimension DLEN.
if 1D type (DTYPE_B=502), DLEN >=7;
if 2D type (DTYPE_B=1), DLEN >= 9.
The array descriptor for the distributed matrix B.
Contains information of mapping of B to memory. Please
see NOTES below for full description and options.
AF (local output) COMPLEX array, dimension LAF.
Auxiliary Fillin Space.
Fillin is created during the factorization routine
PCDTTRF and this is stored in AF. If a linear system
is to be solved using PCDTTRS after the factorization
routine, AF *must not be altered* after the factorization.
LAF (local input) INTEGER
Size of user-input Auxiliary Fillin space AF. Must be >=
2*(NB+2)
If LAF is not large enough, an error code will be returned
and the minimum acceptable size will be returned in AF( 1 )
WORK (local workspace/local output)
COMPLEX temporary workspace. This space may
be overwritten in between calls to routines. WORK must be
the size given in LWORK.
On exit, WORK( 1 ) contains the minimal LWORK.
LWORK (local input or global input) INTEGER
Size of user-input workspace WORK.
If LWORK is too small, the minimal acceptable size will be
returned in WORK(1) and an error code is returned. LWORK>=
10*NPCOL+4*NRHS
INFO (local output) INTEGER
= 0: successful exit
< 0: If the i-th argument is an array and the j-entry had
an illegal value, then INFO = -(i*100+j), if the i-th
argument is a scalar and had an illegal value, then
INFO = -i.
=====================================================================
Restrictions
============
The following are restrictions on the input parameters. Some of these
are temporary and will be removed in future releases, while others
may reflect fundamental technical limitations.
Non-cyclic restriction: VERY IMPORTANT!
P*NB>= mod(JA-1,NB)+N.
The mapping for matrices must be blocked, reflecting the nature
of the divide and conquer algorithm as a task-parallel algorithm.
This formula in words is: no processor may have more than one
chunk of the matrix.
Blocksize cannot be too small:
If the matrix spans more than one processor, the following
restriction on NB, the size of each block on each processor,
must hold:
NB >= 2
The bulk of parallel computation is done on the matrix of size
O(NB) on each processor. If this is too small, divide and conquer
is a poor choice of algorithm.
Submatrix reference:
JA = IB
Alignment restriction that prevents unnecessary communication.
=====================================================================
Notes
=====
If the factorization routine and the solve routine are to be called
separately (to solve various sets of righthand sides using the same
coefficient matrix), the auxiliary space AF *must not be altered*
between calls to the factorization routine and the solve routine.
The best algorithm for solving banded and tridiagonal linear systems
depends on a variety of parameters, especially the bandwidth.
Currently, only algorithms designed for the case N/P >> bw are
implemented. These go by many names, including Divide and Conquer,
Partitioning, domain decomposition-type, etc.
For tridiagonal matrices, it is obvious: N/P >> bw(=1), and so D&C
algorithms are the appropriate choice.
Algorithm description: Divide and Conquer
The Divide and Conqer algorithm assumes the matrix is narrowly
banded compared with the number of equations. In this situation,
it is best to distribute the input matrix A one-dimensionally,
with columns atomic and rows divided amongst the processes.
The basic algorithm divides the tridiagonal matrix up into
P pieces with one stored on each processor,
and then proceeds in 2 phases for the factorization or 3 for the
solution of a linear system.
1) Local Phase:
The individual pieces are factored independently and in
parallel. These factors are applied to the matrix creating
fillin, which is stored in a non-inspectable way in auxiliary
space AF. Mathematically, this is equivalent to reordering
the matrix A as P A P^T and then factoring the principal
leading submatrix of size equal to the sum of the sizes of
the matrices factored on each processor. The factors of
these submatrices overwrite the corresponding parts of A
in memory.
2) Reduced System Phase:
A small ((P-1)) system is formed representing
interaction of the larger blocks, and is stored (as are its
factors) in the space AF. A parallel Block Cyclic Reduction
algorithm is used. For a linear system, a parallel front solve
followed by an analagous backsolve, both using the structure
of the factored matrix, are performed.
3) Backsubsitution Phase:
For a linear system, a local backsubstitution is performed on
each processor in parallel.
Descriptors
===========
Descriptors now have *types* and differ from ScaLAPACK 1.0.
Note: tridiagonal codes can use either the old two dimensional
or new one-dimensional descriptors, though the processor grid in
both cases *must be one-dimensional*. We describe both types below.
Each global data object is described by an associated description
vector. This vector stores the information required to establish
the mapping between an object element and its corresponding process
and memory location.
Let A be a generic term for any 2D block cyclicly distributed array.
Such a global array has an associated description vector DESCA.
In the following comments, the character _ should be read as
"of the global array".
NOTATION STORED IN EXPLANATION
--------------- -------------- --------------------------------------
DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
DTYPE_A = 1.
CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
the BLACS process grid A is distribu-
ted over. The context itself is glo-
bal, but the handle (the integer
value) may vary.
M_A (global) DESCA( M_ ) The number of rows in the global
array A.
N_A (global) DESCA( N_ ) The number of columns in the global
array A.
MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
the rows of the array.
NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
the columns of the array.
RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
row of the array A is distributed.
CSRC_A (global) DESCA( CSRC_ ) The process column over which the
first column of the array A is
distributed.
LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
array. LLD_A >= MAX(1,LOCr(M_A)).
Let K be the number of rows or columns of a distributed matrix,
and assume that its process grid has dimension p x q.
LOCr( K ) denotes the number of elements of K that a process
would receive if K were distributed over the p processes of its
process column.
Similarly, LOCc( K ) denotes the number of elements of K that a
process would receive if K were distributed over the q processes of
its process row.
The values of LOCr() and LOCc() may be determined via a call to the
ScaLAPACK tool function, NUMROC:
LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
An upper bound for these quantities may be computed by:
LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
One-dimensional descriptors:
One-dimensional descriptors are a new addition to ScaLAPACK since
version 1.0. They simplify and shorten the descriptor for 1D
arrays.
Since ScaLAPACK supports two-dimensional arrays as the fundamental
object, we allow 1D arrays to be distributed either over the
first dimension of the array (as if the grid were P-by-1) or the
2nd dimension (as if the grid were 1-by-P). This choice is
indicated by the descriptor type (501 or 502)
as described below.
However, for tridiagonal matrices, since the objects being
distributed are the individual vectors storing the diagonals, we
have adopted the convention that both the P-by-1 descriptor and
the 1-by-P descriptor are allowed and are equivalent for
tridiagonal matrices. Thus, for tridiagonal matrices,
DTYPE_A = 501 or 502 can be used interchangeably
without any other change.
We require that the distributed vectors storing the diagonals of a
tridiagonal matrix be aligned with each other. Because of this, a
single descriptor, DESCA, serves to describe the distribution of
of all diagonals simultaneously.
IMPORTANT NOTE: the actual BLACS grid represented by the
CTXT entry in the descriptor may be *either* P-by-1 or 1-by-P
irrespective of which one-dimensional descriptor type
(501 or 502) is input.
This routine will interpret the grid properly either way.
ScaLAPACK routines *do not support intercontext operations* so that
the grid passed to a single ScaLAPACK routine *must be the same*
for all array descriptors passed to that routine.
NOTE: In all cases where 1D descriptors are used, 2D descriptors
may also be used, since a one-dimensional array is a special case
of a two-dimensional array with one dimension of size unity.
The two-dimensional array used in this case *must* be of the
proper orientation:
If the appropriate one-dimensional descriptor is DTYPEA=501
(1 by P type), then the two dimensional descriptor must
have a CTXT value that refers to a 1 by P BLACS grid;
If the appropriate one-dimensional descriptor is DTYPEA=502
(P by 1 type), then the two dimensional descriptor must
have a CTXT value that refers to a P by 1 BLACS grid.
Summary of allowed descriptors, types, and BLACS grids:
DTYPE 501 502 1 1
BLACS grid 1xP or Px1 1xP or Px1 1xP Px1
-----------------------------------------------------
A OK OK OK NO
B NO OK NO OK
Note that a consequence of this chart is that it is not possible
for *both* DTYPE_A and DTYPE_B to be 2D_type(1), as these lead
to opposite requirements for the orientation of the BLACS grid,
and as noted before, the *same* BLACS context must be used in
all descriptors in a single ScaLAPACK subroutine call.
Let A be a generic term for any 1D block cyclicly distributed array.
Such a global array has an associated description vector DESCA.
In the following comments, the character _ should be read as
"of the global array".
NOTATION STORED IN EXPLANATION
--------------- ---------- ------------------------------------------
DTYPE_A(global) DESCA( 1 ) The descriptor type. For 1D grids,
TYPE_A = 501: 1-by-P grid.
TYPE_A = 502: P-by-1 grid.
CTXT_A (global) DESCA( 2 ) The BLACS context handle, indicating
the BLACS process grid A is distribu-
ted over. The context itself is glo-
bal, but the handle (the integer
value) may vary.
N_A (global) DESCA( 3 ) The size of the array dimension being
distributed.
NB_A (global) DESCA( 4 ) The blocking factor used to distribute
the distributed dimension of the array.
SRC_A (global) DESCA( 5 ) The process row or column over which the
first row or column of the array
is distributed.
Ignored DESCA( 6 ) Ignored for tridiagonal matrices.
Reserved DESCA( 7 ) Reserved for future use.
=====================================================================
Code Developer: Andrew J. Cleary, University of Tennessee.
Current address: Lawrence Livermore National Labs.
This version released: August, 2001.
=====================================================================
..
.. Parameters ..
|
|
|
|
001 SUBROUTINE PCDTTRS( TRANS , N , NRHS , DL , D , DU , JA , DESCA , B , IB ,
002 $DESCB , AF , LAF , WORK , LWORK , INFO )
003
004 * -- ScaLAPACK routine(version 1.7) --
005 * University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
006 * and University of California , Berkeley.
007 * August 7 , 2001
008
009 * .. Scalar Arguments ..
010 CHARACTER TRANS
011 INTEGER IB , INFO , JA , LAF , LWORK , N , NRHS
012 REAL ONE , ZERO
013 PARAMETER( ONE = 1.0E + 0 )
014 PARAMETER( ZERO = 0.0E + 0 )
015 COMPLEX CONE , CZERO
016 PARAMETER( CONE =( 1.0E + 0 , 0.0E + 0 ) )
017 PARAMETER( CZERO =( 0.0E + 0 , 0.0E + 0 ) )
018 INTEGER INT_ONE
019 PARAMETER( INT_ONE = 1 )
020 INTEGER DESCMULT , BIGNUM
021 PARAMETER(DESCMULT = 100 , BIGNUM = DESCMULT * DESCMULT)
022 INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
023 $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
024 PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
025 $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
026 $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
027 * ..
028 * .. Local Scalars ..
029 INTEGER CSRC , FIRST_PROC , ICTXT , ICTXT_NEW , ICTXT_SAVE ,
030 $IDUM2 , IDUM3 , JA_NEW , LLDA , LLDB , MYCOL , MYROW ,
031 $MY_NUM_COLS , NB , NP , NPCOL , NPROW , NP_SAVE ,
032 $ODD_SIZE , PART_OFFSET , PART_SIZE ,
033 $RETURN_CODE , STORE_M_B , STORE_N_A , TEMP ,
034 $WORK_SIZE_MIN
035 * ..
036 * .. Local Arrays ..
037 INTEGER DESCA_1XP( 7 ) , DESCB_PX1( 7 ) ,
038 $PARAM_CHECK( 15 , 3 )
039 * ..
040 * .. External Subroutines ..
041 EXTERNAL BLACS_GRIDINFO , DESC_CONVERT , GLOBCHK ,
042 $PCDTTRSV , PXERBLA , RESHAPE
043 * ..
044 * .. External Functions ..
045 LOGICAL LSAME
046 INTEGER NUMROC
047 COMPLEX CDOTC
048 EXTERNAL CDOTC , LSAME , NUMROC
049 * ..
050 * .. Intrinsic Functions ..
051 INTRINSIC ICHAR , MIN , MOD
052 * ..
053 * .. Executable Statements ..
054
055 * Test the input parameters
056
057 INFO = 0
058
059 * Convert descriptor into standard form for easy access to
060 * parameters , check that grid is of right shape.
061
062 DESCA_1XP( 1 ) = 501
063 DESCB_PX1( 1 ) = 502
064
065 TEMP = DESCA( DTYPE_ )
066 IF( TEMP .EQ. 502 ) THEN
067 * Temporarily set the descriptor type to 1xP type
067
068 DESCA( DTYPE_ ) = 501
069 ENDIF
070
071 CALL DESC_CONVERT( DESCA , DESCA_1XP , RETURN_CODE )
072
073 DESCA( DTYPE_ ) = TEMP
074
075 IF( RETURN_CODE .NE. 0) THEN
075
076 INFO = - ( 8*100 + 2 )
077 ENDIF
078
079 CALL DESC_CONVERT( DESCB , DESCB_PX1 , RETURN_CODE )
080
081 IF( RETURN_CODE .NE. 0) THEN
081
082 INFO = - ( 11*100 + 2 )
083 ENDIF
084
085 * Consistency checks for DESCA and DESCB.
086
087 * Context must be the same
088 IF( DESCA_1XP( 2 ) .NE. DESCB_PX1( 2 ) ) THEN
088
089 INFO = - ( 11*100 + 2 )
090 ENDIF
091
092 * These are alignment restrictions that may or may not be removed
093 * in future releases. - Andy Cleary , April 14 , 1996.
094
095 * Block sizes must be the same
096 IF( DESCA_1XP( 4 ) .NE. DESCB_PX1( 4 ) ) THEN
096
097 INFO = - ( 11*100 + 4 )
098 ENDIF
099
100 * Source processor must be the same
101
102 IF( DESCA_1XP( 5 ) .NE. DESCB_PX1( 5 ) ) THEN
102
103 INFO = - ( 11*100 + 5 )
104 ENDIF
105
106 * Get values out of descriptor for use in code.
107
108 ICTXT = DESCA_1XP( 2 )
109 CSRC = DESCA_1XP( 5 )
110 NB = DESCA_1XP( 4 )
111 LLDA = DESCA_1XP( 6 )
112 STORE_N_A = DESCA_1XP( 3 )
113 LLDB = DESCB_PX1( 6 )
114 STORE_M_B = DESCB_PX1( 3 )
115
116 * Get grid parameters
117
118 CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
119 NP = NPROW * NPCOL
120
121 IF( LSAME( TRANS , 'N' ) ) THEN
121
122 IDUM2 = ICHAR( 'N' )
123 ELSE IF( LSAME( TRANS , 'C' ) ) THEN
123
124 IDUM2 = ICHAR( 'C' )
125 ELSE
125
126 INFO = - 1
127 END IF
128
129 IF( LWORK .LT. - 1) THEN
129
130 INFO = - 15
131 ELSE IF( LWORK .EQ. - 1 ) THEN
131
132 IDUM3 = - 1
133 ELSE
133
134 IDUM3 = 1
135 ENDIF
136
137 IF( N .LT. 0 ) THEN
137
138 INFO = - 2
139 ENDIF
140
141 IF( N + JA - 1 .GT. STORE_N_A ) THEN
141
142 INFO = - ( 8*100 + 6 )
143 ENDIF
144
145 IF( N + IB - 1 .GT. STORE_M_B ) THEN
145
146 INFO = - ( 11*100 + 3 )
147 ENDIF
148
149 IF( LLDB .LT. NB ) THEN
149
150 INFO = - ( 11*100 + 6 )
151 ENDIF
152
153 IF( NRHS .LT. 0 ) THEN
153
154 INFO = - 3
155 ENDIF
156
157 * Current alignment restriction
158
159 IF( JA .NE. IB) THEN
159
160 INFO = - 7
161 ENDIF
162
163 * Argument checking that is specific to Divide & Conquer routine
164
165 IF( NPROW .NE. 1 ) THEN
165
166 INFO = - ( 8*100 + 2 )
167 ENDIF
168
169 IF( N .GT. NP*NB - MOD( JA - 1 , NB )) THEN
169
170 INFO = - ( 2 )
171 CALL PXERBLA( ICTXT ,
172 $ 'PCDTTRS , D&C alg. : only 1 block per proc' ,
173 $ - INFO )
174 RETURN
175 ENDIF
176
177 IF((JA + N - 1.GT.NB) .AND.( NB.LT.2*INT_ONE )) THEN
177
178 INFO = - ( 8*100 + 4 )
179 CALL PXERBLA( ICTXT ,
180 $ 'PCDTTRS , D&C alg. : NB too small' ,
181 $ - INFO )
182 RETURN
183 ENDIF
184
185 WORK_SIZE_MIN =
186 $10*NPCOL + 4*NRHS
187
188 WORK( 1 ) = WORK_SIZE_MIN
189
190 IF( LWORK .LT. WORK_SIZE_MIN ) THEN
190
191 IF( LWORK .NE. - 1 ) THEN
191
192 INFO = - 15
193 CALL PXERBLA( ICTXT ,
194 $ 'PCDTTRS : worksize error' ,
195 $ - INFO )
196 ENDIF
197 RETURN
198 ENDIF
199
200 * Pack params and positions into arrays for global consistency check
201
202 PARAM_CHECK( 15 , 1 ) = DESCB(5)
203 PARAM_CHECK( 14 , 1 ) = DESCB(4)
204 PARAM_CHECK( 13 , 1 ) = DESCB(3)
205 PARAM_CHECK( 12 , 1 ) = DESCB(2)
206 PARAM_CHECK( 11 , 1 ) = DESCB(1)
207 PARAM_CHECK( 10 , 1 ) = IB
208 PARAM_CHECK( 9 , 1 ) = DESCA(5)
209 PARAM_CHECK( 8 , 1 ) = DESCA(4)
210 PARAM_CHECK( 7 , 1 ) = DESCA(3)
211 PARAM_CHECK( 6 , 1 ) = DESCA(1)
212 PARAM_CHECK( 5 , 1 ) = JA
213 PARAM_CHECK( 4 , 1 ) = NRHS
214 PARAM_CHECK( 3 , 1 ) = N
215 PARAM_CHECK( 2 , 1 ) = IDUM3
216 PARAM_CHECK( 1 , 1 ) = IDUM2
217
218 PARAM_CHECK( 15 , 2 ) = 1105
219 PARAM_CHECK( 14 , 2 ) = 1104
220 PARAM_CHECK( 13 , 2 ) = 1103
221 PARAM_CHECK( 12 , 2 ) = 1102
222 PARAM_CHECK( 11 , 2 ) = 1101
223 PARAM_CHECK( 10 , 2 ) = 10
224 PARAM_CHECK( 9 , 2 ) = 805
225 PARAM_CHECK( 8 , 2 ) = 804
226 PARAM_CHECK( 7 , 2 ) = 803
227 PARAM_CHECK( 6 , 2 ) = 801
228 PARAM_CHECK( 5 , 2 ) = 7
229 PARAM_CHECK( 4 , 2 ) = 3
230 PARAM_CHECK( 3 , 2 ) = 2
231 PARAM_CHECK( 2 , 2 ) = 15
232 PARAM_CHECK( 1 , 2 ) = 1
233
234 * Want to find errors with MIN( ) , so if no error , set it to a big
235 * number. If there already is an error , multiply by the the
236 * descriptor multiplier.
237
238 IF( INFO.GE.0 ) THEN
238
239 INFO = BIGNUM
240 ELSE IF( INFO.LT. - DESCMULT ) THEN
240
241 INFO = - INFO
242 ELSE
242
243 INFO = - INFO * DESCMULT
244 END IF
245
246 * Check consistency across processors
247
248 CALL GLOBCHK( ICTXT , 15 , PARAM_CHECK , 15 ,
249 $PARAM_CHECK( 1 , 3 ) , INFO )
250
251 * Prepare output : set info = 0 if no error , and divide by DESCMULT
252 * if error is not in a descriptor entry.
253
254 IF( INFO.EQ.BIGNUM ) THEN
254
255 INFO = 0
256 ELSE IF( MOD( INFO , DESCMULT ) .EQ. 0 ) THEN
256
257 INFO = - INFO / DESCMULT
258 ELSE
258
259 INFO = - INFO
260 END IF
261
262 IF( INFO.LT.0 ) THEN
262
263 CALL PXERBLA( ICTXT , 'PCDTTRS' , - INFO )
264 RETURN
265 END IF
266
267 * Quick return if possible
268
269 IF( N.EQ.0 )
269
270 $ RETURN
271
272 IF( NRHS.EQ.0 )
272
273 $ RETURN
274
275 * Adjust addressing into matrix space to properly get into
276 * the beginning part of the relevant data
277
278 PART_OFFSET = NB*((JA - 1) / (NPCOL*NB) )
279
280 IF((MYCOL - CSRC) .LT.(JA - PART_OFFSET - 1) / NB ) THEN
281 PART_OFFSET = PART_OFFSET + NB
282 ENDIF
283
284 IF( MYCOL .LT. CSRC ) THEN
284
285 PART_OFFSET = PART_OFFSET - NB
286 ENDIF
287
288 * Form a new BLACS grid(the "standard form" grid) with only procs
289 * holding part of the matrix , of size 1xNP where NP is adjusted ,
290 * starting at csrc = 0 , with JA modified to reflect dropped procs.
291
292 * First processor to hold part of the matrix :
293
294 FIRST_PROC = MOD(( JA - 1 ) / NB + CSRC , NPCOL )
295
296 * Calculate new JA one while dropping off unused processors.
297
298 JA_NEW = MOD( JA - 1 , NB ) + 1
299
300 * Save and compute new value of NP
301
302 NP_SAVE = NP
303 NP =( JA_NEW + N - 2 ) / NB + 1
304
305 * Call utility routine that forms "standard-form" grid
306
307 CALL RESHAPE( ICTXT , INT_ONE , ICTXT_NEW , INT_ONE ,
308 $ FIRST_PROC , INT_ONE , NP )
309
310 * Use new context from standard grid as context.
311
312 ICTXT_SAVE = ICTXT
313 ICTXT = ICTXT_NEW
314 DESCA_1XP( 2 ) = ICTXT_NEW
315 DESCB_PX1( 2 ) = ICTXT_NEW
316
317 * Get information about new grid.
318
319 CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
320
321 * Drop out processors that do not have part of the matrix.
322
323 IF( MYROW .LT. 0 ) THEN
323
324 GOTO 1234
325 ENDIF
326
327 * ********************************
328 * Values reused throughout routine
329
330 * User - input value of partition size
331
332 PART_SIZE = NB
333
334 * Number of columns in each processor
335
336 MY_NUM_COLS = NUMROC( N , PART_SIZE , MYCOL , 0 , NPCOL )
337
338 * Offset in columns to beginning of main partition in each proc
339
340 IF( MYCOL .EQ. 0 ) THEN
340
341 PART_OFFSET = PART_OFFSET + MOD( JA_NEW - 1 , PART_SIZE )
342 MY_NUM_COLS = MY_NUM_COLS - MOD(JA_NEW - 1 , PART_SIZE )
343 ENDIF
344
345 * Size of main(or odd) partition in each processor
346
347 ODD_SIZE = MY_NUM_COLS
348 IF( MYCOL .LT. NP - 1 ) THEN
348
349 ODD_SIZE = ODD_SIZE - INT_ONE
350 ENDIF
351
352 * Begin main code
353
354 INFO = 0
355
356 * Call frontsolve routine
357
358 IF( LSAME( TRANS , 'N' ) ) THEN
359
359
360 CALL PCDTTRSV ( 'L' , 'N' , N , NRHS , DL( PART_OFFSET + 1 ) ,
361 $ D( PART_OFFSET + 1 ) , DU( PART_OFFSET + 1 ) , JA_NEW ,
362 $ DESCA_1XP , B , IB , DESCB_PX1 , AF , LAF , WORK ,
363 $ LWORK , INFO )
364
365 ELSE
366
366
367 CALL PCDTTRSV ( 'U' , 'C' , N , NRHS , DL( PART_OFFSET + 1 ) ,
368 $ D( PART_OFFSET + 1 ) , DU( PART_OFFSET + 1 ) , JA_NEW ,
369 $ DESCA_1XP , B , IB , DESCB_PX1 , AF , LAF , WORK ,
370 $ LWORK , INFO )
371
372 ENDIF
373
374 * Call backsolve routine
375
376 IF( LSAME( TRANS , 'C' ) ) THEN
377
377
378 CALL PCDTTRSV ( 'L' , 'C' , N , NRHS , DL( PART_OFFSET + 1 ) ,
379 $ D( PART_OFFSET + 1 ) , DU( PART_OFFSET + 1 ) , JA_NEW ,
380 $ DESCA_1XP , B , IB , DESCB_PX1 , AF , LAF , WORK ,
381 $ LWORK , INFO )
382
383 ELSE
384
384
385 CALL PCDTTRSV ( 'U' , 'N' , N , NRHS , DL( PART_OFFSET + 1 ) ,
386 $ D( PART_OFFSET + 1 ) , DU( PART_OFFSET + 1 ) , JA_NEW ,
387 $ DESCA_1XP , B , IB , DESCB_PX1 , AF , LAF , WORK ,
388 $ LWORK , INFO )
389
390 ENDIF
391 1000 CONTINUE
392
393 * Free BLACS space used to hold standard - form grid.
394
395 IF( ICTXT_SAVE .NE. ICTXT_NEW ) THEN
395
396 CALL BLACS_GRIDEXIT( ICTXT_NEW )
397 ENDIF
398
399 1234 CONTINUE
400
401 * Restore saved input parameters
402
403 ICTXT = ICTXT_SAVE
404 NP = NP_SAVE
405
406 * Output minimum worksize
407
408 WORK( 1 ) = WORK_SIZE_MIN
409
410 RETURN
411
412 * End of PCDTTRS
413
414 END103
41
|
|
Variables in Routine PCDTTRS()
| Summary Report |
| Data Type | Quantity | Size(byte) |
| CHARACTER | 1 | 1 |
| COMPLEX | 3 | 12 |
| INTEGER | 51 | 316 |
| LOGICAL | 1 | 1 |
| REAL | 4 | 16 |
| TOTAL | 60 | 346 |
List of Variables
CHARACTER
COMPLEX
INTEGER
| BIGNUM | BLOCK_CYCLIC_2D | CSRC | CSRC_ | CTXT_ |
| DESCA_1XP( 7 ) | DESCB_PX1( 7 ) | DESCMULT | DLEN_ | DTYPE_ |
| FIRST_PROC | IB | ICTXT | ICTXT_NEW | ICTXT_SAVE |
| IDUM2 | IDUM3 | INFO | INT_ONE | JA |
| JA_NEW | LAF | LLD_ | LLDA | LLDB |
| LWORK | M_ | MB_ | MY_NUM_COLS | MYCOL |
| MYROW | N | N_ | NB | NB_ |
| NP | NP_SAVE | NPCOL | NPROW | NRHS |
| NUMROC | ODD_SIZE | PARAM_CHECK( 15, 3 ) | PART_OFFSET | PART_SIZE |
| RETURN_CODE | RSRC_ | STORE_M_B | STORE_N_A | TEMP |
| WORK_SIZE_MIN | | | | |
LOGICAL
REAL
Variables Dependence Graph Put the mouse over a right hand side variable to display the corresponding line of the dependence | | - | | - | - | | CSRC | <--- | DESCA_1XPCSRC = DESCA_1XP( 5 ) |
| DESCA | <--- | TEMPDESCA( DTYPE_ ) = TEMP |
| DESCA_1XP | <--- | ICTXT_NEWDESCA_1XP( 2 ) = ICTXT_NEW |
| DESCB_PX1 | <--- | ICTXT_NEWDESCB_PX1( 2 ) = ICTXT_NEW |
| FIRST_PROC | <--- | JAFIRST_PROC = MOD( ( JA-1 )/NB+CSRC, NPCOL ), NBFIRST_PROC = MOD( ( JA-1 )/NB+CSRC, NPCOL ), NPCOLFIRST_PROC = MOD( ( JA-1 )/NB+CSRC, NPCOL ), CSRCFIRST_PROC = MOD( ( JA-1 )/NB+CSRC, NPCOL ) |
| ICTXT | <--- | DESCA_1XPICTXT = DESCA_1XP( 2 ), ICTXT_NEWICTXT = ICTXT_NEW, ICTXT_SAVEICTXT = ICTXT_SAVE |
| ICTXT_SAVE | <--- | ICTXTICTXT_SAVE = ICTXT |
| IDUM2 | <--- | NIDUM2 = ICHAR( 'N' ) |
| INFO | <--- | BIGNUMINFO = BIGNUM, DESCMULTINFO = -INFO * DESCMULT{2INFO = -INFO / DESCMULT}, INFOINFO = -INFO{2INFO = -INFO * DESCMULT, 3INFO = -INFO / DESCMULT, 4INFO = -INFO} |
| JA_NEW | <--- | JAJA_NEW = MOD( JA-1, NB ) + 1, NBJA_NEW = MOD( JA-1, NB ) + 1 |
| LLDA | <--- | DESCA_1XPLLDA = DESCA_1XP( 6 ) |
| LLDB | <--- | DESCB_PX1LLDB = DESCB_PX1( 6 ) |
| MY_NUM_COLS | <--- | JA_NEWMY_NUM_COLS = MY_NUM_COLS - MOD(JA_NEW-1, PART_SIZE ), MY_NUM_COLSMY_NUM_COLS = MY_NUM_COLS - MOD(JA_NEW-1, PART_SIZE ), MYCOLMY_NUM_COLS = NUMROC( N, PART_SIZE, MYCOL, 0, NPCOL ), NMY_NUM_COLS = NUMROC( N, PART_SIZE, MYCOL, 0, NPCOL ), NPCOLMY_NUM_COLS = NUMROC( N, PART_SIZE, MYCOL, 0, NPCOL ), NUMROCMY_NUM_COLS = NUMROC( N, PART_SIZE, MYCOL, 0, NPCOL ), PART_SIZEMY_NUM_COLS = NUMROC( N, PART_SIZE, MYCOL, 0, NPCOL ){2MY_NUM_COLS = MY_NUM_COLS - MOD(JA_NEW-1, PART_SIZE )} |
| NB | <--- | DESCA_1XPNB = DESCA_1XP( 4 ) |
| NP | <--- | JA_NEWNP = ( JA_NEW+N-2 )/NB + 1, NNP = ( JA_NEW+N-2 )/NB + 1, NBNP = ( JA_NEW+N-2 )/NB + 1, NP_SAVENP = NP_SAVE, NPCOLNP = NPROW * NPCOL, NPROWNP = NPROW * NPCOL |
| NP_SAVE | <--- | NPNP_SAVE = NP |
| ODD_SIZE | <--- | INT_ONEODD_SIZE = ODD_SIZE - INT_ONE, MY_NUM_COLSODD_SIZE = MY_NUM_COLS, ODD_SIZEODD_SIZE = ODD_SIZE - INT_ONE |
| PARAM_CHECK | <--- | IBPARAM_CHECK( 10, 1 ) = IB, IDUM2PARAM_CHECK( 1, 1 ) = IDUM2, IDUM3PARAM_CHECK( 2, 1 ) = IDUM3, JAPARAM_CHECK( 5, 1 ) = JA, NPARAM_CHECK( 3, 1 ) = N, NRHSPARAM_CHECK( 4, 1 ) = NRHS, DESCAPARAM_CHECK( 9, 1 ) = DESCA(5){2PARAM_CHECK( 8, 1 ) = DESCA(4), 3PARAM_CHECK( 7, 1 ) = DESCA(3), 4PARAM_CHECK( 6, 1 ) = DESCA(1)} |
| PART_OFFSET | <--- | JAPART_OFFSET = NB*( (JA-1)/(NPCOL*NB) ), JA_NEWPART_OFFSET = PART_OFFSET+MOD( JA_NEW-1, PART_SIZE ), NBPART_OFFSET = NB*( (JA-1)/(NPCOL*NB) ){2PART_OFFSET = PART_OFFSET + NB, 3PART_OFFSET = PART_OFFSET - NB}, NPCOLPART_OFFSET = NB*( (JA-1)/(NPCOL*NB) ), PART_OFFSETPART_OFFSET = PART_OFFSET + NB{2PART_OFFSET = PART_OFFSET - NB, 3PART_OFFSET = PART_OFFSET+MOD( JA_NEW-1, PART_SIZE )}, PART_SIZEPART_OFFSET = PART_OFFSET+MOD( JA_NEW-1, PART_SIZE ) |
| PART_SIZE | <--- | NBPART_SIZE = NB |
| STORE_M_B | <--- | DESCB_PX1STORE_M_B = DESCB_PX1( 3 ) |
| STORE_N_A | <--- | DESCA_1XPSTORE_N_A = DESCA_1XP( 3 ) |
| TEMP | <--- | DTYPE_TEMP = DESCA( DTYPE_ ), DESCATEMP = DESCA( DTYPE_ ) |
| WORK | <--- | WORK_SIZE_MINWORK( 1 ) = WORK_SIZE_MIN{2WORK( 1 ) = WORK_SIZE_MIN} |
| WORK_SIZE_MIN | <--- | NPCOLWORK_SIZE_MIN =, NRHSWORK_SIZE_MIN = |
|
|
Analysis elements of the routine PCDTTRS() Put the mouse over each element to display detailed matching information
Assigned variables |
| | | BIGNUM , BLOCK_CYCLIC_2D , CONE , CSRC , CSRC_ , CTXT_ , CZERO , DESCA_1XP , DESCB_PX1 , DESCMULT , DLEN_ , DTYPE_ , FIRST_PROC , ICTXT , ICTXT_SAVE , IDUM2 , IDUM3 , INFO , INT_ONE , JA_NEW , LLD_ , LLDA , LLDB , M_ , MB_ , MY_NUM_COLS , N_ , NB , NB_ , NP , NP_SAVE , ODD_SIZE , ONE , PARAM_CHECK , PART_OFFSET , PART_SIZE , RSRC_ , STORE_M_B , STORE_N_A , TEMP , WORK , WORK_SIZE_MIN , ZERO |
|
Active variables |
| | | AF , B , BIGNUM , BLOCK_CYCLIC_2D , CDOTC , CONE , CSRC , CSRC_ , CTXT_ , CZERO , D , DESCA , DESCA_1XP , DESCB , DESCB_PX1 , DESCMULT , DL , DLEN_ , DTYPE_ , DU , FIRST_PROC , IB , ICTXT , ICTXT_NEW , ICTXT_SAVE , IDUM2 , IDUM3 , INFO , INT_ONE , JA , JA_NEW , LAF , LLD_ , LLDA , LLDB , LSAME , LWORK , M_ , MB_ , MY_NUM_COLS , MYCOL , MYROW , N , N_ , NB , NB_ , NP , NP_SAVE , NPCOL , NPROW , NRHS , NUMROC , ODD_SIZE , ONE , PARAM_CHECK , PART_OFFSET , PART_SIZE , RETURN_CODE , RSRC_ , STORE_M_B , STORE_N_A , TEMP , TRANS , WORK , WORK_SIZE_MIN , ZERO |
|
Allocated variables [ statement : associated variable ] |
| | new | : a, about, Calculate, compute, Use |
|
Desallocated variables [ statement : associated variable ] |
| | free | : BLACS |
|
Accessed arrays [ array name : associated index ] |
| | D | : PART_OFFSET+1 , PART_OFFSET+1 , PART_OFFSET+1 , PART_OFFSET+1 |
| | DESCA | : 1 , 3 , 4 , 5 , DTYPE_ , DTYPE_ , DTYPE_ |
| | DESCA_1XP | : 1 , 2 , 2 , 2 , 3 , 4 , 4 , 5 , 5 , 6 , 7 |
| | DESCB | : 1 , 2 , 3 , 4 , 5 |
| | DESCB_PX1 | : 1 , 2 , 2 , 3 , 4 , 5 , 6 , 7 |
| | DL | : PART_OFFSET+1 , PART_OFFSET+1 , PART_OFFSET+1 , PART_OFFSET+1 |
| | DU | : PART_OFFSET+1 , PART_OFFSET+1 , PART_OFFSET+1 , PART_OFFSET+1 |
| | LSAME | : TRANS, 'C' , TRANS, 'C' , TRANS, 'N' , TRANS, 'N' |
| | NUMROC | : N, PART_SIZE, MYCOL, 0, NPCOL |
| | PARAM_CHECK | : 1, 1 , 1, 2 , 1, 3 , 10, 1 , 10, 2 , 11, 1 , 11, 2 , 12, 1 , 12, 2 , 13, 1 , 13, 2 , 14, 1 , 14, 2 , 15, 1 , 15, 2 , 15, 3 , 2, 1 , 2, 2 , 3, 1 , 3, 2 , 4, 1 , 4, 2 , 5, 1 , 5, 2 , 6, 1 , 6, 2 , 7, 1 , 7, 2 , 8, 1 , 8, 2 , 9, 1 , 9, 2 |
| | WORK | : 1 , 1 |
|
Conditional statements [ statement : associated predicate ] |
| | do | : ( not have part of the matrix. ) |
| | for | : ( easy access to ) , ( DESCA and DESCB. ) , ( use in code. ) , ( global consistency check ) |
| | if | : ( TEMP .EQ. 502 ) , ( RETURN_CODE .NE. 0 ) , ( RETURN_CODE .NE. 0 ) , ( (DESCA_1XP( 2 ) .NE. DESCB_PX1( 2 ) ) ) , ( (DESCA_1XP( 4 ) .NE. DESCB_PX1( 4 ) ) ) , ( (DESCA_1XP( 5 ) .NE. DESCB_PX1( 5 ) ) ) , ( (LSAME( TRANS , 'N' ) ) ) , ( (LSAME( TRANS , 'C' ) ) ) , ( LWORK .LT. - 1 ) , ( LWORK .EQ. - 1 ) , ( N .LT. 0 ) , ( N+JA-1 .GT. STORE_N_A ) , ( N+IB-1 .GT. STORE_M_B ) , ( LLDB .LT. NB ) , ( NRHS .LT. 0 ) , ( JA .NE. IB ) , ( NPROW .NE. 1 ) , ( (N .GT. NP*NB - MOD( JA - 1 , NB )) ) , ( ((JA+N - 1.GT.NB) .AND. ( NB.LT.2*INT_ONE )) ) , ( LWORK .LT. WORK_SIZE_MIN ) , ( LWORK .NE. - 1 ) , ( no error , set it to a big ) , ( there already is an error , multiply by the the ) , ( INFO.GE.0 ) , ( INFO.LT. - DESCMULT ) , ( no error , and divide by DESCMULT ) , ( error is not in a descriptor entry. ) , ( INFO.EQ.BIGNUM ) , ( (MOD( INFO , DESCMULT ) .EQ. 0 ) ) , ( INFO.LT.0 ) , ( possible ) , ( N.EQ.0 ) , ( NRHS.EQ.0 ) , ( ((MYCOL - CSRC) .LT. (JA - PART_OFFSET - 1) / NB ) ) , ( MYCOL .LT. CSRC ) , ( MYROW .LT. 0 ) , ( MYCOL .EQ. 0 ) , ( MYCOL .LT. NP - 1 ) , ( (LSAME( TRANS , 'N' ) ) ) , ( (LSAME( TRANS , 'C' ) ) ) , ( ICTXT_SAVE .NE. ICTXT_NEW ) |
| | while | : ( dropping off unused processors. ) |
|
| List of variables | BIGNUM BLOCK_CYCLIC_2D CDOTC CONE CSRC CSRC_ CTXT_
| CZERO DESCA DESCA_1XP( 7 ) DESCB_PX1( 7 ) DESCMULT DLEN_ DTYPE_ FIRST_PROC
| IB ICTXT ICTXT_NEW ICTXT_SAVE IDUM2 IDUM3 INFO INT_ONE
| JA JA_NEW LAF LLD_ LLDA LLDB LSAME LWORK
| M_ MB_ MY_NUM_COLS MYCOL MYROW N N_ NB
| NB_ NP NP_SAVE NPCOL NPROW NRHS NUMROC ODD_SIZE
| ONE PARAM_CHECK( 15, 3 ) PART_OFFSET PART_SIZE RETURN_CODE RSRC_ STORE_M_B STORE_N_A
| TEMP TRANS WORK WORK_SIZE_MIN ZERO | | close
| |
BIGNUM
BLOCK_CYCLIC_2D
CDOTC
CONE
CSRC
CSRC_
CTXT_
CZERO
DESCA
DESCA_1XP( 7 )
DESCB_PX1( 7 )
DESCMULT
DLEN_
DTYPE_
FIRST_PROC
IB
ICTXT
ICTXT_NEW
ICTXT_SAVE
IDUM2
IDUM3
INFO
INT_ONE
JA
JA_NEW
LAF
LLD_
LLDA
LLDB
LSAME
LWORK
M_
MB_
MY_NUM_COLS
MYCOL
MYROW
N
N_
NB
NB_
NP
NP_SAVE
NPCOL
NPROW
NRHS
NUMROC
ODD_SIZE
ONE
PARAM_CHECK( 15, 3 )
PART_OFFSET
PART_SIZE
RETURN_CODE
RSRC_
STORE_M_B
STORE_N_A
TEMP
TRANS
WORK
WORK_SIZE_MIN
ZERO
33
| |