ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdgetri.f
Go to the documentation of this file.
1  SUBROUTINE pdgetri( N, A, IA, JA, DESCA, IPIV, WORK, LWORK,
2  $ IWORK, LIWORK, INFO )
3 *
4 * -- ScaLAPACK routine (version 1.7.4) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * v1.7.4: May 10, 2006
8 * v1.7: May 1, 1997
9 *
10 * .. Scalar Arguments ..
11  INTEGER IA, INFO, JA, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCA( * ), IPIV( * ), IWORK( * )
15  DOUBLE PRECISION A( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PDGETRI computes the inverse of a distributed matrix using the LU
22 * factorization computed by PDGETRF. This method inverts U and then
23 * computes the inverse of sub( A ) = A(IA:IA+N-1,JA:JA+N-1) denoted
24 * InvA by solving the system InvA*L = inv(U) for InvA.
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 * N (global input) INTEGER
84 * The number of rows and columns to be operated on, i.e. the
85 * order of the distributed submatrix sub( A ). N >= 0.
86 *
87 * A (local input/local output) DOUBLE PRECISION pointer into the
88 * local memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
89 * On entry, the local pieces of the L and U obtained by the
90 * factorization sub( A ) = P*L*U computed by PDGETRF. On
91 * exit, if INFO = 0, sub( A ) contains the inverse of the
92 * original distributed matrix sub( A ).
93 *
94 * IA (global input) INTEGER
95 * The row index in the global array A indicating the first
96 * row of sub( A ).
97 *
98 * JA (global input) INTEGER
99 * The column index in the global array A indicating the
100 * first column of sub( A ).
101 *
102 * DESCA (global and local input) INTEGER array of dimension DLEN_.
103 * The array descriptor for the distributed matrix A.
104 *
105 * IPIV (local input) INTEGER array, dimension LOCr(M_A)+MB_A
106 * keeps track of the pivoting information. IPIV(i) is the
107 * global row index the local row i was swapped with. This
108 * array is tied to the distributed matrix A.
109 *
110 * WORK (local workspace/local output) DOUBLE PRECISION array,
111 * dimension (LWORK)
112 * On exit, WORK(1) returns the minimal and optimal LWORK.
113 *
114 * LWORK (local or global input) INTEGER
115 * The dimension of the array WORK.
116 * LWORK is local input and must be at least
117 * LWORK = LOCr(N+MOD(IA-1,MB_A))*NB_A. WORK is used to keep a
118 * copy of at most an entire column block of sub( A ).
119 *
120 * If LWORK = -1, then LWORK is global input and a workspace
121 * query is assumed; the routine only calculates the minimum
122 * and optimal size for all work arrays. Each of these
123 * values is returned in the first entry of the corresponding
124 * work array, and no error message is issued by PXERBLA.
125 *
126 * IWORK (local workspace/local output) INTEGER array,
127 * dimension (LIWORK)
128 * On exit, IWORK(1) returns the minimal and optimal LIWORK.
129 *
130 * LIWORK (local or global input) INTEGER
131 * The dimension of the array IWORK used as workspace for
132 * physically transposing the pivots.
133 * LIWORK is local input and must be at least
134 * if NPROW == NPCOL then
135 * LIWORK = LOCc( N_A + MOD(JA-1, NB_A) ) + NB_A,
136 * else
137 * LIWORK = LOCc( N_A + MOD(JA-1, NB_A) ) +
138 * MAX( CEIL(CEIL(LOCr(M_A)/MB_A)/(LCM/NPROW)),
139 * NB_A )
140 * where LCM is the least common multiple of process
141 * rows and columns (NPROW and NPCOL).
142 * end if
143 *
144 * If LIWORK = -1, then LIWORK is global input and a workspace
145 * query is assumed; the routine only calculates the minimum
146 * and optimal size for all work arrays. Each of these
147 * values is returned in the first entry of the corresponding
148 * work array, and no error message is issued by PXERBLA.
149 *
150 * INFO (global output) INTEGER
151 * = 0: successful exit
152 * < 0: If the i-th argument is an array and the j-entry had
153 * an illegal value, then INFO = -(i*100+j), if the i-th
154 * argument is a scalar and had an illegal value, then
155 * INFO = -i.
156 * > 0: If INFO = K, U(IA+K-1,IA+K-1) is exactly zero; the
157 * matrix is singular and its inverse could not be
158 * computed.
159 *
160 * =====================================================================
161 *
162 * .. Parameters ..
163  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
164  $ lld_, mb_, m_, nb_, n_, rsrc_
165  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
166  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
167  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
168  DOUBLE PRECISION ZERO, ONE
169  parameter( zero = 0.0d+0, one = 1.0d+0 )
170 * ..
171 * .. Local Scalars ..
172  LOGICAL LQUERY
173  INTEGER I, IACOL, IAROW, ICOFF, ICTXT, IROFF, IW, J,
174  $ jb, jn, lcm, liwmin, lwmin, mp, mycol, myrow,
175  $ nn, np, npcol, nprow, nq
176 * ..
177 * .. Local Arrays ..
178  INTEGER DESCW( DLEN_ ), IDUM1( 2 ), IDUM2( 2 )
179 * ..
180 * .. External Subroutines ..
181  EXTERNAL blacs_gridinfo, chk1mat, descset, pchk1mat,
182  $ pdgemm, pdlacpy, pdlaset, pdlapiv,
183  $ pdtrsm, pdtrtri, pxerbla
184 * ..
185 * .. External Functions ..
186  INTEGER ICEIL, ILCM, INDXG2P, NUMROC
187  EXTERNAL iceil, ilcm, indxg2p, numroc
188 * ..
189 * .. Intrinsic Functions ..
190  INTRINSIC dble, max, min, mod
191 * ..
192 * .. Executable Statements ..
193 *
194 * Get grid parameters
195 *
196  ictxt = desca( ctxt_ )
197  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
198 *
199 * Test the input parameters
200 *
201  info = 0
202  IF( nprow.EQ.-1 ) THEN
203  info = -(500+ctxt_)
204  ELSE
205  CALL chk1mat( n, 1, n, 1, ia, ja, desca, 5, info )
206  IF( info.EQ.0 ) THEN
207  iroff = mod( ia-1, desca( mb_ ) )
208  icoff = mod( ja-1, desca( nb_ ) )
209  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
210  $ nprow )
211  np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
212  lwmin = np * desca( nb_ )
213 *
214  mp = numroc( desca( m_ ), desca( mb_ ), myrow,
215  $ desca( rsrc_ ), nprow )
216  nq = numroc( desca( n_ ), desca( nb_ ), mycol,
217  $ desca( csrc_ ), npcol )
218  IF( nprow.EQ.npcol ) THEN
219  liwmin = nq + desca( nb_ )
220  ELSE
221 *
222 * Use the formula for the workspace given in PxLAPIV
223 * to compute the minimum size LIWORK for IWORK
224 *
225 * The formula in PxLAPIV is
226 * LDW = LOCc( M_P + MOD(IP-1, MB_P) ) +
227 * MB_P * CEIL( CEIL(LOCr(M_P)/MB_P) / (LCM/NPROW) )
228 *
229 * where
230 * M_P is the global length of the pivot vector
231 * MP = DESCA( M_ ) + DESCA( MB_ ) * NPROW
232 * I_P is IA
233 * I_P = IA
234 * MB_P is the block size use for the block cyclic distribution of the
235 * pivot vector
236 * MB_P = DESCA (MB_ )
237 * LOCc ( . )
238 * NUMROC ( . , DESCA ( NB_ ), MYCOL, DESCA ( CSRC_ ), NPCOL )
239 * LOCr ( . )
240 * NUMROC ( . , DESCA ( MB_ ), MYROW, DESCA ( RSRC_ ), NPROW )
241 * CEIL ( X / Y )
242 * ICEIL( X, Y )
243 * LCM
244 * LCM = ILCM( NPROW, NPCOL )
245 *
246  lcm = ilcm( nprow, npcol )
247  liwmin = numroc( desca( m_ ) + desca( mb_ ) * nprow
248  $ + mod( ia - 1, desca( mb_ ) ), desca( nb_ ),
249  $ mycol, desca( csrc_ ), npcol ) +
250  $ max( desca( mb_ ) * iceil( iceil(
251  $ numroc( desca( m_ ) + desca( mb_ ) * nprow,
252  $ desca( mb_ ), myrow, desca( rsrc_ ), nprow ),
253  $ desca( mb_ ) ), lcm / nprow ), desca( nb_ ) )
254 *
255  END IF
256 *
257  work( 1 ) = dble( lwmin )
258  iwork( 1 ) = liwmin
259  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
260  IF( iroff.NE.icoff .OR. iroff.NE.0 ) THEN
261  info = -4
262  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
263  info = -(500+nb_)
264  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
265  info = -8
266  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
267  info = -10
268  END IF
269  END IF
270  IF( lwork.EQ.-1 ) THEN
271  idum1( 1 ) = -1
272  ELSE
273  idum1( 1 ) = 1
274  END IF
275  idum2( 1 ) = 8
276  IF( liwork.EQ.-1 ) THEN
277  idum1( 2 ) = -1
278  ELSE
279  idum1( 2 ) = 1
280  END IF
281  idum2( 2 ) = 10
282  CALL pchk1mat( n, 1, n, 1, ia, ja, desca, 5, 2, idum1, idum2,
283  $ info )
284  END IF
285 *
286  IF( info.NE.0 ) THEN
287  CALL pxerbla( ictxt, 'PDGETRI', -info )
288  RETURN
289  ELSE IF( lquery ) THEN
290  RETURN
291  END IF
292 *
293 * Quick return if possible
294 *
295  IF( n.EQ.0 )
296  $ RETURN
297 *
298 * Form inv(U). If INFO > 0 from PDTRTRI, then U is singular,
299 * and the inverse is not computed.
300 *
301  CALL pdtrtri( 'Upper', 'Non-unit', n, a, ia, ja, desca, info )
302  IF( info.GT.0 )
303  $ RETURN
304 *
305 * Define array descriptor for working array WORK
306 *
307  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
308  nn = ( ( ja+n-2 ) / desca( nb_ ) ) * desca( nb_ ) + 1
309  iacol = indxg2p( nn, desca( nb_ ), mycol, desca( csrc_ ), npcol )
310  CALL descset( descw, n+iroff, desca( nb_ ), desca( mb_ ),
311  $ desca( nb_ ), iarow, iacol, ictxt, max( 1, np ) )
312  iw = iroff + 1
313 *
314 * Solve the equation inv(A)*L=inv(U) for inv(A) using blocked code.
315 *
316  DO 10 j = nn, jn+1, -desca( nb_ )
317  jb = min( desca( nb_ ), ja+n-j )
318  i = ia + j - ja
319 *
320 * Copy current block column of L to WORK and replace with zeros.
321 *
322  CALL pdlacpy( 'Lower', ja+n-1-j, jb, a, i+1, j, desca,
323  $ work, iw+j-ja+1, 1, descw )
324  CALL pdlaset( 'Lower', ja+n-1-j, jb, zero, zero, a, i+1, j,
325  $ desca )
326 *
327 * Compute current block column of inv(A).
328 *
329  IF( j+jb.LE.ja+n-1 )
330  $ CALL pdgemm( 'No transpose', 'No transpose', n, jb,
331  $ ja+n-j-jb, -one, a, ia, j+jb, desca, work,
332  $ iw+j+jb-ja, 1, descw, one, a, ia, j, desca )
333  CALL pdtrsm( 'Right', 'Lower', 'No transpose', 'Unit', n, jb,
334  $ one, work, iw+j-ja, 1, descw, a, ia, j, desca )
335  descw( csrc_ ) = mod( descw( csrc_ ) + npcol - 1, npcol )
336 *
337  10 CONTINUE
338 *
339 * Handle the last block of columns separately
340 *
341  jb = jn-ja+1
342 *
343 * Copy current block column of L to WORK and replace with zeros.
344 *
345  CALL pdlacpy( 'Lower', n-1, jb, a, ia+1, ja, desca, work, iw+1,
346  $ 1, descw )
347  CALL pdlaset( 'Lower', n-1, jb, zero, zero, a, ia+1, ja, desca )
348 *
349 * Compute current block column of inv(A).
350 *
351  IF( ja+jb.LE.ja+n-1 )
352  $ CALL pdgemm( 'No transpose', 'No transpose', n, jb,
353  $ n-jb, -one, a, ia, ja+jb, desca, work, iw+jb, 1,
354  $ descw, one, a, ia, ja, desca )
355  CALL pdtrsm( 'Right', 'Lower', 'No transpose', 'Unit', n, jb,
356  $ one, work, iw, 1, descw, a, ia, ja, desca )
357 *
358 * Use the row pivots and apply them to the columns of the global
359 * matrix.
360 *
361 *
362 * JL: I do not get why the size of the PIVOT vector is DESCA( M_ ) + DESCA( MB_ )*NPROW
363 * should be DESCA( M_ ) + DESCA( MB_ ) no?
364 *
365  CALL descset( descw, desca( m_ ) + desca( mb_ )*nprow, 1,
366  $ desca( mb_ ), 1, desca( rsrc_ ), mycol, ictxt,
367  $ mp+desca( mb_ ) )
368  CALL pdlapiv( 'Backward', 'Columns', 'Column', n, n, a, ia,
369  $ ja, desca, ipiv, ia, 1, descw, iwork )
370 *
371  work( 1 ) = dble( lwmin )
372  iwork( 1 ) = liwmin
373 *
374  RETURN
375 *
376 * End of PDGETRI
377 *
378  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdgetri
subroutine pdgetri(N, A, IA, JA, DESCA, IPIV, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdgetri.f:3
pdlapiv
subroutine pdlapiv(DIREC, ROWCOL, PIVROC, M, N, A, IA, JA, DESCA, IPIV, IP, JP, DESCIP, IWORK)
Definition: pdlapiv.f:3
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
pdlaset
subroutine pdlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pdblastst.f:6862
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pdlacpy
subroutine pdlacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pdlacpy.f:3
min
#define min(A, B)
Definition: pcgemr.c:181
pdtrtri
subroutine pdtrtri(UPLO, DIAG, N, A, IA, JA, DESCA, INFO)
Definition: pdtrtri.f:2