ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
psinvchk.f
Go to the documentation of this file.
1  SUBROUTINE psinvchk( MATTYP, N, A, IA, JA, DESCA, IASEED, ANORM,
2  $ FRESID, RCOND, WORK )
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  REAL ANORM, FRESID, RCOND
12 * ..
13 * .. Array Arguments ..
14  CHARACTER*3 MATTYP
15  INTEGER DESCA( * )
16  REAL A( * ), WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * PSINVCHK 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) symmetric Positive Definite,
92 * if MATTYP = 'LPD' then (Lower) symmetric 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) REAL 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 PSGETRI, PSPOTRI or PSTRTRI.
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) REAL
118 * The 1-norm of the original matrix sub( A ).
119 *
120 * FRESID (global output) REAL
121 * The inversion residual.
122 *
123 * RCOND (global output) REAL
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) REAL 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 PSLANGE, PSLANSY and PSLANTR.
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  REAL ZERO, ONE
142  parameter( zero = 0.0e+0, one = 1.0e+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  REAL AUXNORM, EPS, NRMINVAXA, TEMP
150 * ..
151 * .. Local Arrays ..
152  INTEGER DESCW( DLEN_ )
153 * ..
154 * .. External Subroutines ..
155  EXTERNAL blacs_gridinfo, descset, infog2l, psgemm,
156  $ pslaset, psmatgen, pssymm, pstrmm
157 * ..
158 * .. External Functions ..
159  LOGICAL LSAMEN
160  INTEGER ICEIL, NUMROC
161  REAL PSLAMCH, PSLANGE, PSLANSY, PSLANTR
162  EXTERNAL iceil, lsamen, numroc, pslamch, pslange,
163  $ pslansy, pslantr
164 * ..
165 * .. Intrinsic Functions ..
166  INTRINSIC max, min, mod
167 * ..
168 * .. Executable Statements ..
169 *
170  eps = pslamch( 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 = pslange( '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 = pslantr( '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 = 'S'
200  diag = 'D'
201  auxnorm = pslansy( '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 psmatgen( 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 pslaset( 'Lower', n-1, jb, zero, zero, work, iw+1,
241  $ 1, descw )
242  ELSE
243  CALL pslaset( 'Upper', jb-1, jb-1, zero, zero, work, iw,
244  $ 2, descw )
245  END IF
246  ELSE
247  CALL psmatgen( 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 psgemm( '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 pstrmm( '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 pssymm( '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 = pslange( '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 psmatgen( 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 pslaset( 'Lower', ja+n-j-1, jb, zero, zero,
312  $ work, iw+j-ja+1, 1, descw )
313  ELSE
314  CALL pslaset( 'All', j-ja, jb, zero, zero, work, iw,
315  $ 1, descw )
316  CALL pslaset( 'Upper', jb-1, jb-1, zero, zero,
317  $ work, iw+j-ja, 2, descw )
318  END IF
319  ELSE
320  CALL psmatgen( 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 psgemm( '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 pstrmm( '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 pssymm( '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 = pslange( '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 PSINVCHK
382 *
383  END
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
psinvchk
subroutine psinvchk(MATTYP, N, A, IA, JA, DESCA, IASEED, ANORM, FRESID, RCOND, WORK)
Definition: psinvchk.f:3
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
psmatgen
subroutine psmatgen(ICTXT, AFORM, DIAG, M, N, MB, NB, A, LDA, IAROW, IACOL, ISEED, IROFF, IRNUM, ICOFF, ICNUM, MYROW, MYCOL, NPROW, NPCOL)
Definition: psmatgen.f:4
pslaset
subroutine pslaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: psblastst.f:6863
min
#define min(A, B)
Definition: pcgemr.c:181