SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ pzinvchk()

subroutine pzinvchk ( character*3  mattyp,
integer  n,
complex*16, dimension( * )  a,
integer  ia,
integer  ja,
integer, dimension( * )  desca,
integer  iaseed,
double precision  anorm,
double precision  fresid,
double precision  rcond,
complex*16, dimension( * )  work 
)

Definition at line 1 of file pzinvchk.f.

3*
4* -- ScaLAPACK routine (version 1.7) --
5* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6* and University of California, Berkeley.
7* May 28, 2001
8*
9* .. Scalar Arguments ..
10 INTEGER IA, IASEED, JA, N
11 DOUBLE PRECISION ANORM, FRESID, RCOND
12* ..
13* .. Array Arguments ..
14 CHARACTER*3 MATTYP
15 INTEGER DESCA( * )
16 COMPLEX*16 A( * ), WORK( * )
17* ..
18*
19* Purpose
20* =======
21*
22* PZINVCHK computes the scaled residual
23*
24* || sub( A ) * inv( sub( A ) ) - I || / ( || sub( A ) || * N * eps ),
25*
26* where sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1). to check the result
27* returned by the matrix inversion routines.
28*
29* Notes
30* =====
31*
32* Each global data object is described by an associated description
33* vector. This vector stores the information required to establish
34* the mapping between an object element and its corresponding process
35* and memory location.
36*
37* Let A be a generic term for any 2D block cyclicly distributed array.
38* Such a global array has an associated description vector DESCA.
39* In the following comments, the character _ should be read as
40* "of the global array".
41*
42* NOTATION STORED IN EXPLANATION
43* --------------- -------------- --------------------------------------
44* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
45* DTYPE_A = 1.
46* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
47* the BLACS process grid A is distribu-
48* ted over. The context itself is glo-
49* bal, but the handle (the integer
50* value) may vary.
51* M_A (global) DESCA( M_ ) The number of rows in the global
52* array A.
53* N_A (global) DESCA( N_ ) The number of columns in the global
54* array A.
55* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
56* the rows of the array.
57* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
58* the columns of the array.
59* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
60* row of the array A is distributed.
61* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
62* first column of the array A is
63* distributed.
64* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
65* array. LLD_A >= MAX(1,LOCr(M_A)).
66*
67* Let K be the number of rows or columns of a distributed matrix,
68* and assume that its process grid has dimension p x q.
69* LOCr( K ) denotes the number of elements of K that a process
70* would receive if K were distributed over the p processes of its
71* process column.
72* Similarly, LOCc( K ) denotes the number of elements of K that a
73* process would receive if K were distributed over the q processes of
74* its process row.
75* The values of LOCr() and LOCc() may be determined via a call to the
76* ScaLAPACK tool function, NUMROC:
77* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
78* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
79* An upper bound for these quantities may be computed by:
80* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
81* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
82*
83* Arguments
84* =========
85*
86* MATTYP (global input) CHARACTER*3
87* The type of the distributed matrix to be generated:
88* if MATTYP = 'GEN' then GENeral matrix,
89* if MATTYP = 'UTR' then Upper TRiangular matrix,
90* if MATTYP = 'LTR' then Lower TRiangular matrix,
91* if MATTYP = 'UPD' then (Upper) Hermitian Positive Definite,
92* if MATTYP = 'LPD' then (Lower) Hermitian Positive Definite.
93*
94* N (global input) INTEGER
95* The number of rows and columns to be operated on, i.e. the
96* order of the distributed submatrix sub( A ). N >= 0.
97*
98* A (local input) COMPLEX*16 pointer into the local memory
99* to an array of local dimension (LLD_A, LOCc(JA+N-1)). On
100* entry, sub( A ) contains the distributed matrix inverse
101* computed by PZGETRI, PZPOTRI or PZTRTRI.
102*
103* IA (global input) INTEGER
104* The row index in the global array A indicating the first
105* row of sub( A ).
106*
107* JA (global input) INTEGER
108* The column index in the global array A indicating the
109* first column of sub( A ).
110*
111* DESCA (global and local input) INTEGER array of dimension DLEN_.
112* The array descriptor for the distributed matrix A.
113*
114* IASEED (global input) INTEGER
115* Seed for the random generation of sub( A ).
116*
117* ANORM (global input) DOUBLE PRECISION
118* The 1-norm of the original matrix sub( A ).
119*
120* FRESID (global output) DOUBLE PRECISION
121* The inversion residual.
122*
123* RCOND (global output) DOUBLE PRECISION
124* The condition number of the original distributed matrix.
125* RCOND = || sub( A ) ||.|| sub( A )^{-1} || where ||A||
126* denotes the 1-norm of A.
127*
128* WORK (local workspace) COMPLEX*16 array, dimension
129* MAX(2*LOCr(N_A+MOD(IA-1,MB_A))*MB_A, LDW)
130* where LDW is the workspace requirement for the norm computa-
131* tions, see PZLANGE, PZLANHE, PZLANSY and PZLANTR.
132*
133* =====================================================================
134*
135* .. Parameters ..
136 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
137 $ LLD_, MB_, M_, NB_, N_, RSRC_
138 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
139 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
140 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
141 COMPLEX*16 ZERO, ONE
142 parameter( zero = 0.0d+0, one = 1.0d+0 )
143* ..
144* .. Local Scalars ..
145 CHARACTER AFORM, DIAG, UPLO
146 INTEGER ICTXT, ICURCOL, ICURROW, II, IIA, IPW, IROFF,
147 $ IW, J, JB, JJA, JN, KK, MYCOL, MYROW, NP,
148 $ NPCOL, NPROW
149 DOUBLE PRECISION AUXNORM, EPS, NRMINVAXA, TEMP
150* ..
151* .. Local Arrays ..
152 INTEGER DESCW( DLEN_ )
153* ..
154* .. External Subroutines ..
155 EXTERNAL blacs_gridinfo, descset, infog2l, pzgemm,
156 $ pzhemm, pzlaset, pzmatgen, pztrmm
157* ..
158* .. External Functions ..
159 LOGICAL LSAMEN
160 INTEGER ICEIL, NUMROC
161 DOUBLE PRECISION PDLAMCH, PZLANGE, PZLANHE, PZLANTR
162 EXTERNAL iceil, lsamen, numroc, pdlamch, pzlange,
163 $ pzlanhe, pzlansy, pzlantr
164* ..
165* .. Intrinsic Functions ..
166 INTRINSIC max, min, mod
167* ..
168* .. Executable Statements ..
169*
170 eps = pdlamch( desca( ctxt_ ), 'eps' )
171*
172* Get grid parameters
173*
174 ictxt = desca( ctxt_ )
175 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
176*
177* Compute the condition number
178*
179 IF( lsamen( 1, mattyp( 1:1 ), 'U' ) ) THEN
180 uplo = 'U'
181 ELSE
182 uplo = 'L'
183 END IF
184*
185 IF( lsamen( 3, mattyp, 'GEN' ) ) THEN
186*
187 aform = 'N'
188 diag = 'D'
189 auxnorm = pzlange( '1', n, n, a, ia, ja, desca, work )
190*
191 ELSE IF( lsamen( 2, mattyp( 2:3 ), 'TR' ) ) THEN
192*
193 aform = 'N'
194 diag = 'D'
195 auxnorm = pzlantr( '1', uplo, 'Non unit', n, n, a, ia, ja,
196 $ desca, work )
197 ELSE IF( lsamen( 2, mattyp( 2:3 ), 'PD' ) ) THEN
198*
199 aform = 'H'
200 diag = 'D'
201 auxnorm = pzlanhe( '1', uplo, n, a, ia, ja, desca, work )
202*
203 END IF
204 rcond = anorm*auxnorm
205*
206* Compute inv(A)*A
207*
208 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
209 $ icurrow, icurcol )
210*
211* Define array descriptor for working array WORK
212*
213 iroff = mod( ia-1, desca( mb_ ) )
214 np = numroc( n+iroff, desca( mb_ ), myrow, icurrow, nprow )
215 CALL descset( descw, n+iroff, desca( nb_ ), desca( mb_ ),
216 $ desca( nb_ ), icurrow, icurcol, desca( ctxt_ ),
217 $ max( 1, np ) )
218 ipw = descw( lld_ ) * descw( nb_ ) + 1
219*
220 IF( myrow.EQ.icurrow ) THEN
221 ii = iroff + 1
222 np = np - iroff
223 ELSE
224 ii = 1
225 END IF
226 jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
227 jb = jn - ja + 1
228*
229* Handle first block separately, regenerate a block of columns of A
230*
231 iw = iroff + 1
232 IF( mycol.EQ.icurcol ) THEN
233 IF( lsamen( 2, mattyp( 2:3 ), 'TR' ) ) THEN
234 CALL pzmatgen( ictxt, aform, diag, desca( m_ ), desca( n_ ),
235 $ descw( mb_ ), descw( nb_ ), work,
236 $ descw( lld_ ), desca( rsrc_ ),
237 $ desca( csrc_ ), iaseed, iia-1, np,
238 $ jja-1, jb, myrow, mycol, nprow, npcol )
239 IF( lsamen( 3, mattyp, 'UTR' ) ) THEN
240 CALL pzlaset( 'Lower', n-1, jb, zero, zero, work, iw+1,
241 $ 1, descw )
242 ELSE
243 CALL pzlaset( 'Upper', jb-1, jb-1, zero, zero, work, iw,
244 $ 2, descw )
245 END IF
246 ELSE
247 CALL pzmatgen( ictxt, aform, diag, desca( m_ ), desca( n_ ),
248 $ descw( mb_ ), descw( nb_ ), work( ipw ),
249 $ descw( lld_ ), desca( rsrc_ ),
250 $ desca( csrc_ ), iaseed,
251 $ iia-1, np, jja-1, jb, myrow, mycol, nprow,
252 $ npcol )
253 END IF
254 END IF
255*
256* Multiply A^{-1}*A
257*
258 IF( lsamen( 3, mattyp, 'GEN' ) ) THEN
259*
260 CALL pzgemm( 'No tranpose', 'No transpose', n, jb, n, one, a,
261 $ ia, ja, desca, work( ipw ), iw, 1, descw, zero,
262 $ work, iw, 1, descw )
263*
264 ELSE IF( lsamen( 2, mattyp( 2:3 ), 'TR' ) ) THEN
265*
266 CALL pztrmm( 'Left', uplo, 'No tranpose', 'Non unit', n, jb,
267 $ one, a, ia, ja, desca, work, iw, 1, descw )
268*
269 ELSE IF( lsamen( 2, mattyp( 2:3 ), 'PD' ) ) THEN
270*
271 CALL pzhemm( 'Left', uplo, n, jb, one, a, ia, ja, desca,
272 $ work(ipw), iw, 1, descw, zero, work, iw, 1,
273 $ descw )
274*
275 END IF
276*
277* subtract the identity matrix to the diagonal block of these cols
278*
279 IF( myrow.EQ.icurrow .AND. mycol.EQ.icurcol ) THEN
280 DO 10 kk = 0, jb-1
281 work( ii+kk*(descw(lld_)+1) ) =
282 $ work( ii+kk*(descw( lld_ )+1) )-one
283 10 CONTINUE
284 END IF
285*
286 nrminvaxa = pzlange( '1', n, jb, work, iw, 1, descw, work( ipw ) )
287*
288 IF( myrow.EQ.icurrow )
289 $ ii = ii + jb
290 IF( mycol.EQ.icurcol )
291 $ jja = jja + jb
292 icurrow = mod( icurrow+1, nprow )
293 icurcol = mod( icurcol+1, npcol )
294 descw( csrc_ ) = icurcol
295*
296 DO 30 j = jn+1, ja+n-1, desca( nb_ )
297*
298 jb = min( n-j+ja, desca( nb_ ) )
299*
300* regenerate a block of columns of A
301*
302 IF( mycol.EQ.icurcol ) THEN
303 IF( lsamen( 2, mattyp( 2:3 ), 'TR' ) ) THEN
304 CALL pzmatgen( ictxt, aform, diag, desca( m_ ),
305 $ desca( n_ ), descw( mb_ ), descw( nb_ ),
306 $ work, descw( lld_ ), desca( rsrc_ ),
307 $ desca( csrc_ ),
308 $ iaseed, iia-1, np, jja-1, jb, myrow,
309 $ mycol, nprow, npcol )
310 IF( lsamen( 3, mattyp, 'UTR' ) ) THEN
311 CALL pzlaset( 'Lower', ja+n-j-1, jb, zero, zero,
312 $ work, iw+j-ja+1, 1, descw )
313 ELSE
314 CALL pzlaset( 'All', j-ja, jb, zero, zero, work, iw,
315 $ 1, descw )
316 CALL pzlaset( 'Upper', jb-1, jb-1, zero, zero,
317 $ work, iw+j-ja, 2, descw )
318 END IF
319 ELSE
320 CALL pzmatgen( ictxt, aform, diag, desca( m_ ),
321 $ desca( n_ ), descw( mb_ ), descw( nb_ ),
322 $ work( ipw ), descw( lld_ ),
323 $ desca( rsrc_ ), desca( csrc_ ), iaseed,
324 $ iia-1, np,
325 $ jja-1, jb, myrow, mycol, nprow, npcol )
326 END IF
327 END IF
328*
329* Multiply A^{-1}*A
330*
331 IF( lsamen( 3, mattyp, 'GEN' ) ) THEN
332*
333 CALL pzgemm( 'No tranpose', 'No transpose', n, jb, n, one,
334 $ a, ia, ja, desca, work( ipw ), iw, 1, descw,
335 $ zero, work, iw, 1, descw )
336*
337 ELSE IF( lsamen( 2, mattyp(2:3), 'TR' ) ) THEN
338*
339 CALL pztrmm( 'Left', uplo, 'No tranpose', 'Non unit', n, jb,
340 $ one, a, ia, ja, desca, work, iw, 1, descw )
341*
342 ELSE IF( lsamen( 2, mattyp( 2:3 ), 'PD' ) ) THEN
343*
344 CALL pzhemm( 'Left', uplo, n, jb, one, a, ia, ja, desca,
345 $ work(ipw), iw, 1, descw, zero, work, iw, 1,
346 $ descw )
347*
348 END IF
349*
350* subtract the identity matrix to the diagonal block of these
351* cols
352*
353 IF( myrow.EQ.icurrow .AND. mycol.EQ.icurcol ) THEN
354 DO 20 kk = 0, jb-1
355 work( ii+kk*(descw( lld_ )+1) ) =
356 $ work( ii+kk*(descw( lld_ )+1) ) - one
357 20 CONTINUE
358 END IF
359*
360* Compute the 1-norm of these JB cols
361*
362 temp = pzlange( '1', n, jb, work, iw, 1, descw, work( ipw ) )
363 nrminvaxa = max( temp, nrminvaxa )
364*
365 IF( myrow.EQ.icurrow )
366 $ ii = ii + jb
367 IF( mycol.EQ.icurcol )
368 $ jja = jja + jb
369 icurrow = mod( icurrow+1, nprow )
370 icurcol = mod( icurcol+1, npcol )
371 descw( csrc_ ) = icurcol
372*
373 30 CONTINUE
374*
375* Compute the scaled residual
376*
377 fresid = nrminvaxa / ( n * eps * anorm )
378*
379 RETURN
380*
381* End of PZINVCHK
382*
pzmatgen
subroutine pzmatgen(ictxt, aform, diag, m, n, mb, nb, a, lda, iarow, iacol, iseed, iroff, irnum, icoff, icnum, myrow, mycol, nprow, npcol)
Definition pzmatgen.f:4
descset
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
Definition descset.f:3
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
lsamen
logical function lsamen(n, ca, cb)
Definition pblastst.f:1457
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
pdlamch
double precision function pdlamch(ictxt, cmach)
Definition pdblastst.f:6769
pzlaset
subroutine pzlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pzblastst.f:7509
pzlange
double precision function pzlange(norm, m, n, a, ia, ja, desca, work)
Definition pzlange.f:3
pzlanhe
double precision function pzlanhe(norm, uplo, n, a, ia, ja, desca, work)
Definition pzlanhe.f:3
pzlansy
double precision function pzlansy(norm, uplo, n, a, ia, ja, desca, work)
Definition pzlansy.f:3
pzlantr
double precision function pzlantr(norm, uplo, diag, m, n, a, ia, ja, desca, work)
Definition pzlantr.f:3
Here is the call graph for this function:
Here is the caller graph for this function: