ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pctrtrs.f
Go to the documentation of this file.
1  SUBROUTINE pctrtrs( UPLO, TRANS, DIAG, N, NRHS, A, IA, JA, DESCA,
2  $ B, IB, JB, DESCB, INFO )
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 DIAG, TRANS, UPLO
11  INTEGER IA, IB, INFO, JA, JB, N, NRHS
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCA( * ), DESCB( * )
15  COMPLEX A( * ), B( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PCTRTRS solves a triangular system of the form
22 *
23 * sub( A ) * X = sub( B ) or sub( A )**T * X = sub( B ) or
24 *
25 * sub( A )**H * X = sub( B ),
26 *
27 * where sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1) and is a triangular
28 * distributed matrix of order N, and B(IB:IB+N-1,JB:JB+NRHS-1) is an
29 * N-by-NRHS distributed matrix denoted by sub( B ). A check is made
30 * to verify that sub( A ) is nonsingular.
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 * UPLO (global input) CHARACTER
90 * = 'U': sub( A ) is upper triangular;
91 * = 'L': sub( A ) is lower triangular.
92 *
93 * TRANS (global input) CHARACTER
94 * Specifies the form of the system of equations:
95 * = 'N': Solve sub( A ) * X = sub( B ) (No transpose)
96 * = 'T': Solve sub( A )**T * X = sub( B ) (Transpose)
97 * = 'C': Solve sub( A )**H * X = sub( B ) (Conjugate transpose)
98 *
99 * DIAG (global input) CHARACTER
100 * = 'N': sub( A ) is non-unit triangular;
101 * = 'U': sub( A ) is unit triangular.
102 *
103 * N (global input) INTEGER
104 * The number of rows and columns to be operated on i.e the
105 * order of the distributed submatrix sub( A ). N >= 0.
106 *
107 * NRHS (global input) INTEGER
108 * The number of right hand sides, i.e., the number of columns
109 * of the distributed matrix sub( B ). NRHS >= 0.
110 *
111 * A (local input) COMPLEX pointer into the local memory
112 * to an array of dimension (LLD_A,LOCc(JA+N-1) ). This array
113 * contains the local pieces of the distributed triangular
114 * matrix sub( A ). If UPLO = 'U', the leading N-by-N upper
115 * triangular part of sub( A ) contains the upper triangular
116 * matrix, and the strictly lower triangular part of sub( A )
117 * is not referenced. If UPLO = 'L', the leading N-by-N lower
118 * triangular part of sub( A ) contains the lower triangular
119 * matrix, and the strictly upper triangular part of sub( A )
120 * is not referenced. If DIAG = 'U', the diagonal elements of
121 * sub( A ) are also not referenced and are assumed to be 1.
122 *
123 * IA (global input) INTEGER
124 * The row index in the global array A indicating the first
125 * row of sub( A ).
126 *
127 * JA (global input) INTEGER
128 * The column index in the global array A indicating the
129 * first column of sub( A ).
130 *
131 * DESCA (global and local input) INTEGER array of dimension DLEN_.
132 * The array descriptor for the distributed matrix A.
133 *
134 * B (local input/local output) COMPLEX pointer into the
135 * local memory to an array of dimension
136 * (LLD_B,LOCc(JB+NRHS-1)). On entry, this array contains the
137 * local pieces of the right hand side distributed matrix
138 * sub( B ). On exit, if INFO = 0, sub( B ) is overwritten by
139 * the solution matrix X.
140 *
141 * IB (global input) INTEGER
142 * The row index in the global array B indicating the first
143 * row of sub( B ).
144 *
145 * JB (global input) INTEGER
146 * The column index in the global array B indicating the
147 * first column of sub( B ).
148 *
149 * DESCB (global and local input) INTEGER array of dimension DLEN_.
150 * The array descriptor for the distributed matrix B.
151 *
152 * INFO (output) INTEGER
153 * = 0: successful exit
154 * < 0: If the i-th argument is an array and the j-entry had
155 * an illegal value, then INFO = -(i*100+j), if the i-th
156 * argument is a scalar and had an illegal value, then
157 * INFO = -i.
158 * > 0: If INFO = i, the i-th diagonal element of sub( A ) is
159 * zero, indicating that the submatrix is singular and the
160 * solutions X have not been computed.
161 *
162 * =====================================================================
163 *
164 * .. Parameters ..
165  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
166  $ lld_, mb_, m_, nb_, n_, rsrc_
167  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
168  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
169  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
170  COMPLEX ZERO, ONE
171  parameter( zero = 0.0e+0, one = 1.0e+0 )
172 * ..
173 * .. Local Scalars ..
174  LOGICAL NOTRAN, NOUNIT, UPPER
175  INTEGER I, IAROW, IBROW, ICOFFA, ICTXT, ICURCOL,
176  $ icurrow, iroffa, iroffb, idum, ii, ioffa, j,
177  $ jblk, jj, jn, lda, ll, mycol, myrow, npcol,
178  $ nprow
179 * ..
180 * .. Local Arrays ..
181  INTEGER IDUM1( 3 ), IDUM2( 3 )
182 * ..
183 * .. External Subroutines ..
184  EXTERNAL blacs_gridinfo, chk1mat, igamx2d, infog2l,
185  $ pchk2mat, pctrsm, pxerbla
186 * ..
187 * .. External Functions ..
188  LOGICAL LSAME
189  INTEGER ICEIL, INDXG2P
190  EXTERNAL iceil, indxg2p, lsame
191 * ..
192 * .. Intrinsic Functions ..
193  INTRINSIC ichar, min, mod
194 * ..
195 * .. Executable Statements ..
196 *
197 * Get grid parameters
198 *
199  ictxt = desca( ctxt_ )
200  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
201 *
202 * Test input parameters
203 *
204  info = 0
205  IF( nprow.EQ.-1 ) THEN
206  info = -907
207  ELSE
208  upper = lsame( uplo, 'U' )
209  nounit = lsame( diag, 'N' )
210  notran = lsame( trans, 'N' )
211 *
212  CALL chk1mat( n, 4, n, 4, ia, ja, desca, 9, info )
213  CALL chk1mat( n, 4, nrhs, 5, ib, jb, descb, 13, info )
214  IF( info.EQ.0 ) THEN
215  iroffa = mod( ia-1, desca( mb_ ) )
216  icoffa = mod( ja-1, desca( nb_ ) )
217  iroffb = mod( ib-1, descb( mb_ ) )
218  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
219  $ nprow )
220  ibrow = indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
221  $ nprow )
222  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
223  info = -1
224  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND.
225  $ .NOT.lsame( trans, 'C' ) ) THEN
226  info = -2
227  ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
228  info = -3
229  ELSE IF( iroffa.NE.icoffa .OR. iroffa.NE.0 ) THEN
230  info = -8
231  ELSE IF( iroffa.NE.iroffb .OR. iarow.NE.ibrow ) THEN
232  info = -11
233  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
234  info = -904
235  ELSE IF( descb( mb_ ).NE.desca( nb_ ) ) THEN
236  info = -1304
237  END IF
238  END IF
239 *
240  IF( upper ) THEN
241  idum1( 1 ) = ichar( 'U' )
242  ELSE
243  idum1( 1 ) = ichar( 'L' )
244  END IF
245  idum2( 1 ) = 1
246  IF( notran ) THEN
247  idum1( 2 ) = ichar( 'N' )
248  ELSE IF( lsame( trans, 'T' ) ) THEN
249  idum1( 2 ) = ichar( 'T' )
250  ELSE IF( lsame( trans, 'C' ) ) THEN
251  idum1( 2 ) = ichar( 'C' )
252  END IF
253  idum2( 2 ) = 2
254  IF( nounit ) THEN
255  idum1( 3 ) = ichar( 'N' )
256  ELSE
257  idum1( 3 ) = ichar( 'D' )
258  END IF
259  idum2( 3 ) = 3
260  CALL pchk2mat( n, 4, n, 4, ia, ja, desca, 9, n, 4, nrhs, 5,
261  $ ib, jb, descb, 13, 3, idum1, idum2, info )
262  END IF
263 *
264  IF( info.NE.0 ) THEN
265  CALL pxerbla( ictxt, 'PCTRTRS', -info )
266  RETURN
267  END IF
268 *
269 * Quick return if possible
270 *
271  IF( n.EQ.0 .OR. nrhs.EQ.0 )
272  $ RETURN
273 *
274 * Check for singularity if non-unit.
275 *
276  IF( nounit ) THEN
277  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol,
278  $ ii, jj, icurrow, icurcol )
279  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
280  lda = desca( lld_ )
281  ioffa = ii + ( jj - 1 ) * lda
282 *
283 * Handle first block separately
284 *
285  jblk = jn-ja+1
286  IF( myrow.EQ.icurrow .AND. mycol.EQ.icurcol ) THEN
287  ll = ioffa
288  DO 10 i = 0, jblk-1
289  IF( a( ll ).EQ.zero .AND. info.EQ.0 )
290  $ info = i + 1
291  ll = ioffa + lda + 1
292  10 CONTINUE
293  END IF
294  IF( myrow.EQ.icurrow )
295  $ ioffa = ioffa + jblk
296  IF( mycol.EQ.icurcol )
297  $ ioffa = ioffa + jblk*lda
298  icurrow = mod( icurrow+1, nprow )
299  icurcol = mod( icurcol+1, npcol )
300 *
301 * Loop over remaining blocks of columns
302 *
303  DO 30 j = jn+1, ja+n-1, desca( nb_ )
304  jblk = min( ja+n-j, desca( nb_ ) )
305  IF( myrow.EQ.icurrow .AND. mycol.EQ.icurcol ) THEN
306  ll = ioffa
307  DO 20 i = 0, jblk-1
308  IF( a( ll ).EQ.zero .AND. info.EQ.0 )
309  $ info = j + i - ja + 1
310  ll = ioffa + lda + 1
311  20 CONTINUE
312  END IF
313  IF( myrow.EQ.icurrow )
314  $ ioffa = ioffa + jblk
315  IF( mycol.EQ.icurcol )
316  $ ioffa = ioffa + jblk*lda
317  icurrow = mod( icurrow+1, nprow )
318  icurcol = mod( icurcol+1, npcol )
319  30 CONTINUE
320  CALL igamx2d( ictxt, 'All', ' ', 1, 1, info, 1, idum, idum,
321  $ -1, -1, mycol )
322  IF( info.NE.0 )
323  $ RETURN
324  END IF
325 *
326 * Solve A * x = b, A**T * x = b, or A**H * x = b.
327 *
328  CALL pctrsm( 'Left', uplo, trans, diag, n, nrhs, one, a, ia, ja,
329  $ desca, b, ib, jb, descb )
330 *
331  RETURN
332 *
333 * End of PCTRTRS
334 *
335  END
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pchk2mat
subroutine pchk2mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, MB, MBPOS0, NB, NBPOS0, IB, JB, DESCB, DESCBPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:175
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pctrtrs
subroutine pctrtrs(UPLO, TRANS, DIAG, N, NRHS, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO)
Definition: pctrtrs.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
min
#define min(A, B)
Definition: pcgemr.c:181