OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pctrtri.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine pctrtri (uplo, diag, n, a, ia, ja, desca, info)

Function/Subroutine Documentation

◆ pctrtri()

subroutine pctrtri ( character uplo,
character diag,
integer n,
complex, dimension( * ) a,
integer ia,
integer ja,
integer, dimension( * ) desca,
integer info )

Definition at line 1 of file pctrtri.f.

2*
3* -- ScaLAPACK 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 DIAG, UPLO
10 INTEGER IA, INFO, JA, N
11* ..
12* .. Array Arguments ..
13 INTEGER DESCA( * )
14 COMPLEX A( * )
15* ..
16*
17* Purpose
18* =======
19*
20* PCTRTRI computes the inverse of a upper or lower triangular
21* distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1).
22*
23* Notes
24* =====
25*
26* Each global data object is described by an associated description
27* vector. This vector stores the information required to establish
28* the mapping between an object element and its corresponding process
29* and memory location.
30*
31* Let A be a generic term for any 2D block cyclicly distributed array.
32* Such a global array has an associated description vector DESCA.
33* In the following comments, the character _ should be read as
34* "of the global array".
35*
36* NOTATION STORED IN EXPLANATION
37* --------------- -------------- --------------------------------------
38* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
39* DTYPE_A = 1.
40* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
41* the BLACS process grid A is distribu-
42* ted over. The context itself is glo-
43* bal, but the handle (the integer
44* value) may vary.
45* M_A (global) DESCA( M_ ) The number of rows in the global
46* array A.
47* N_A (global) DESCA( N_ ) The number of columns in the global
48* array A.
49* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
50* the rows of the array.
51* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
52* the columns of the array.
53* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
54* row of the array A is distributed.
55* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
56* first column of the array A is
57* distributed.
58* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
59* array. LLD_A >= MAX(1,LOCr(M_A)).
60*
61* Let K be the number of rows or columns of a distributed matrix,
62* and assume that its process grid has dimension p x q.
63* LOCr( K ) denotes the number of elements of K that a process
64* would receive if K were distributed over the p processes of its
65* process column.
66* Similarly, LOCc( K ) denotes the number of elements of K that a
67* process would receive if K were distributed over the q processes of
68* its process row.
69* The values of LOCr() and LOCc() may be determined via a call to the
70* ScaLAPACK tool function, NUMROC:
71* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
72* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
73* An upper bound for these quantities may be computed by:
74* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
75* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
76*
77* Arguments
78* =========
79*
80* UPLO (global input) CHARACTER
81* Specifies whether the distributed matrix sub( A ) is upper
82* or lower triangular:
83* = 'U': Upper triangular,
84* = 'L': Lower triangular.
85*
86* DIAG (global input) CHARACTER
87* Specifies whether or not the distributed matrix sub( A )
88* is unit triangular:
89* = 'N': Non-unit triangular,
90* = 'U': Unit triangular.
91*
92* N (global input) INTEGER
93* The number of rows and columns to be operated on, i.e. the
94* order of the distributed submatrix sub( A ). N >= 0.
95*
96* A (local input/local output) COMPLEX pointer into the
97* local memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
98* On entry, this array contains the local pieces of the
99* triangular matrix sub( A ). If UPLO = 'U', the leading
100* N-by-N upper triangular part of the matrix sub( A ) contains
101* the upper triangular matrix to be inverted, and the strictly
102* lower triangular part of sub( A ) is not referenced.
103* If UPLO = 'L', the leading N-by-N lower triangular part of
104* the matrix sub( A ) contains the lower triangular matrix,
105* and the strictly upper triangular part of sub( A ) is not
106* referenced.
107* On exit, the (triangular) inverse of the original matrix.
108*
109* IA (global input) INTEGER
110* The row index in the global array A indicating the first
111* row of sub( A ).
112*
113* JA (global input) INTEGER
114* The column index in the global array A indicating the
115* first column of sub( A ).
116*
117* DESCA (global and local input) INTEGER array of dimension DLEN_.
118* The array descriptor for the distributed matrix A.
119*
120* INFO (global output) INTEGER
121* = 0: successful exit
122* < 0: If the i-th argument is an array and the j-entry had
123* an illegal value, then INFO = -(i*100+j), if the i-th
124* argument is a scalar and had an illegal value, then
125* INFO = -i.
126* > 0: If INFO = K, A(IA+K-1,JA+K-1) is exactly zero. The
127* triangular matrix sub( A ) is singular and its
128* inverse can not be computed.
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 COMPLEX ZERO, ONE
139 parameter( zero = 0.0e+0, one = 1.0e+0 )
140* ..
141* .. Local Scalars ..
142 LOGICAL NOUNIT, UPPER
143 INTEGER I, ICOFF, ICTXT, IROFF, ICURCOL, ICURROW,
144 $ IDUMMY, II, IOFFA, J, JB, JJ, JN, LDA, MYCOL,
145 $ MYROW, NN, NPCOL, NPROW
146* ..
147* .. Local Arrays ..
148 INTEGER IDUM1( 2 ), IDUM2( 2 )
149* ..
150* .. External Subroutines ..
151 EXTERNAL blacs_gridinfo, chk1mat, igamx2d, infog2l,
152 $ pchk1mat, pctrti2, pctrmm, pctrsm,
153 $ pxerbla
154* ..
155* .. External Functions ..
156 LOGICAL LSAME
157 INTEGER ICEIL
158 EXTERNAL iceil, lsame
159* ..
160* .. Intrinsic Functions ..
161 INTRINSIC ichar, min, mod
162* ..
163* .. Executable Statements ..
164*
165* Get grid parameters
166*
167 ictxt = desca( ctxt_ )
168 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
169*
170* Test input parameters
171*
172 info = 0
173 IF( nprow.EQ.-1 ) THEN
174 info = -(700+ctxt_)
175 ELSE
176 upper = lsame( uplo, 'U' )
177 nounit = lsame( diag, 'n' )
178*
179 CALL CHK1MAT( N, 3, N, 3, IA, JA, DESCA, 7, INFO )
180.EQ. IF( INFO0 ) THEN
181 IROFF = MOD( IA-1, DESCA( MB_ ) )
182 ICOFF = MOD( JA-1, DESCA( NB_ ) )
183.NOT..AND..NOT. IF( UPPER LSAME( UPLO, 'l' ) ) THEN
184 INFO = -1
185.NOT..AND..NOT. ELSE IF( NOUNIT LSAME( DIAG, 'u' ) ) THEN
186 INFO = -2
187.NE..OR..NE. ELSE IF( IROFFICOFF IROFF0 ) THEN
188 INFO = -6
189.NE. ELSE IF( DESCA( MB_ )DESCA( NB_ ) ) THEN
190 INFO = -(700+NB_)
191 END IF
192 END IF
193*
194 IF( UPPER ) THEN
195 IDUM1( 1 ) = ICHAR( 'u' )
196 ELSE
197 IDUM1( 1 ) = ICHAR( 'l' )
198 END IF
199 IDUM2( 1 ) = 1
200 IF( NOUNIT ) THEN
201 IDUM1( 2 ) = ICHAR( 'n' )
202 ELSE
203 IDUM1( 2 ) = ICHAR( 'u' )
204 END IF
205 IDUM2( 2 ) = 2
206*
207 CALL PCHK1MAT( N, 3, N, 3, IA, JA, DESCA, 7, 2, IDUM1, IDUM2,
208 $ INFO )
209 END IF
210*
211.NE. IF( INFO0 ) THEN
212 CALL PXERBLA( ICTXT, 'pctrtri', -INFO )
213 RETURN
214 END IF
215*
216* Quick return if possible
217*
218.EQ. IF( N0 )
219 $ RETURN
220*
221* Check for singularity if non-unit.
222*
223 JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 )
224 IF( NOUNIT ) THEN
225 CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL,
226 $ II, JJ, ICURROW, ICURCOL )
227*
228* Handle first block separately
229*
230 JB = JN-JA+1
231 LDA = DESCA( LLD_ )
232.EQ..AND..EQ. IF( MYROWICURROW MYCOLICURCOL ) THEN
233 IOFFA = II+(JJ-1)*LDA
234 DO 10 I = 0, JB-1
235.EQ..AND..EQ. IF( A( IOFFA )ZERO INFO0 )
236 $ INFO = I + 1
237 IOFFA = IOFFA + LDA + 1
238 10 CONTINUE
239 END IF
240.EQ. IF( MYROWICURROW )
241 $ II = II + JB
242.EQ. IF( MYCOLICURCOL )
243 $ JJ = JJ + JB
244 ICURROW = MOD( ICURROW+1, NPROW )
245 ICURCOL = MOD( ICURCOL+1, NPCOL )
246*
247* Loop over remaining blocks of columns
248*
249 DO 30 J = JN+1, JA+N-1, DESCA( NB_ )
250 JB = MIN( JA+N-J, DESCA( NB_ ) )
251.EQ..AND..EQ. IF( MYROWICURROW MYCOLICURCOL ) THEN
252 IOFFA = II+(JJ-1)*LDA
253 DO 20 I = 0, JB-1
254.EQ..AND..EQ. IF( A( IOFFA )ZERO INFO0 )
255 $ INFO = J + I - JA + 1
256 IOFFA = IOFFA + LDA + 1
257 20 CONTINUE
258 END IF
259.EQ. IF( MYROWICURROW )
260 $ II = II + JB
261.EQ. IF( MYCOLICURCOL )
262 $ JJ = JJ + JB
263 ICURROW = MOD( ICURROW+1, NPROW )
264 ICURCOL = MOD( ICURCOL+1, NPCOL )
265 30 CONTINUE
266 CALL IGAMX2D( ICTXT, 'all', ' ', 1, 1, INFO, 1, IDUMMY,
267 $ IDUMMY, -1, -1, MYCOL )
268.NE. IF( INFO0 )
269 $ RETURN
270 END IF
271*
272* Use blocked code
273*
274 IF( UPPER ) THEN
275*
276* Compute inverse of upper triangular matrix
277*
278 JB = JN-JA+1
279*
280* Handle first block of column separately
281*
282 CALL PCTRTI2( UPLO, DIAG, JB, A, IA, JA, DESCA, INFO )
283*
284* Loop over remaining block of columns
285*
286 DO 40 J = JN+1, JA+N-1, DESCA( NB_ )
287 JB = MIN( DESCA( NB_ ), JA+N-J )
288 I = IA + J - JA
289*
290* Compute rows 1:j-1 of current block column
291*
292 CALL PCTRMM( 'left', UPLO, 'no transpose', DIAG, J-JA, JB,
293 $ ONE, A, IA, JA, DESCA, A, IA, J, DESCA )
294 CALL PCTRSM( 'right', UPLO, 'no transpose', DIAG, J-JA,
295 $ JB, -ONE, A, I, J, DESCA, A, IA, J, DESCA )
296*
297* Compute inverse of current diagonal block
298*
299 CALL PCTRTI2( UPLO, DIAG, JB, A, I, J, DESCA, INFO )
300*
301 40 CONTINUE
302*
303 ELSE
304*
305* Compute inverse of lower triangular matrix
306*
307 NN = ( ( JA+N-2 ) / DESCA( NB_ ) )*DESCA( NB_ ) + 1
308 DO 50 J = NN, JN+1, -DESCA( NB_ )
309 JB = MIN( DESCA( NB_ ), JA+N-J )
310 I = IA + J - JA
311.LE. IF( J+JBJA+N-1 ) THEN
312*
313* Compute rows j+jb:ja+n-1 of current block column
314*
315 CALL PCTRMM( 'left', UPLO, 'no transpose', DIAG,
316 $ JA+N-J-JB, JB, ONE, A, I+JB, J+JB, DESCA,
317 $ A, I+JB, J, DESCA )
318 CALL PCTRSM( 'right', UPLO, 'no transpose', DIAG,
319 $ JA+N-J-JB, JB, -ONE, A, I, J, DESCA,
320 $ A, I+JB, J, DESCA )
321 END IF
322*
323* Compute inverse of current diagonal block
324*
325 CALL PCTRTI2( UPLO, DIAG, JB, A, I, J, DESCA, INFO )
326*
327 50 CONTINUE
328*
329* Handle the last block of columns separately
330*
331 JB = JN-JA+1
332.LE. IF( JA+JBJA+N-1 ) THEN
333*
334* Compute rows ja+jb:ja+n-1 of current block column
335*
336 CALL PCTRMM( 'left', UPLO, 'no transpose', DIAG, N-JB, JB,
337 $ ONE, A, IA+JB, JA+JB, DESCA, A, IA+JB, JA,
338 $ DESCA )
339 CALL PCTRSM( 'right', UPLO, 'no transpose', DIAG, N-JB, JB,
340 $ -ONE, A, IA, JA, DESCA, A, IA+JB, JA, DESCA )
341 END IF
342*
343* Compute inverse of current diagonal block
344*
345 CALL PCTRTI2( UPLO, DIAG, JB, A, IA, JA, DESCA, INFO )
346*
347 END IF
348*
349 RETURN
350*
351* End PCTRTRI
352*
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
min
#define min(a, b)
Definition macros.h:20
chk1mat
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
Definition mpi.f:1577
pxerbla
subroutine pxerbla(contxt, srname, info)
Definition mpi.f:1600
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
pctrsm
subroutine pctrsm(side, uplo, transa, diag, m, n, alpha, a, ia, ja, desca, b, ib, jb, descb)
Definition mpi.f:1470
pchk1mat
subroutine pchk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, nextra, ex, expos, info)
Definition pchkxmat.f:3
pctrti2
subroutine pctrti2(uplo, diag, n, a, ia, ja, desca, info)
Definition pctrti2.f:2
pctrtri
subroutine pctrtri(uplo, diag, n, a, ia, ja, desca, info)
Definition pctrtri.f:2