1 SUBROUTINE pzlaschk( SYMM, DIAG, N, NRHS, X, IX, JX, DESCX,
2 $ IASEED, IA, JA, DESCA, IBSEED, ANORM, RESID,
12 INTEGER IA, IASEED, IBSEED, IX, JA, JX, N, NRHS
13 DOUBLE PRECISION ANORM, RESID
16 INTEGER DESCA( * ), DESCX( * )
17 COMPLEX*16 WORK( * ), X( * )
149 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
150 $ LLD_, MB_, M_, NB_, N_, RSRC_
151 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
152 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
153 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
155 PARAMETER ( ONE = ( 1.0d+0, 0.0d+0 ),
156 $ zero = ( 0.0d+0, 0.0d+0 ) )
159 INTEGER IACOL, IAROW, IB, ICOFF, ICTXT, ICURCOL, IDUMM,
160 $ II, IIA, IIX, IOFFX, IPA, IPB, IPW, IPX, IROFF,
161 $ ixcol, ixrow, j, jbrhs, jj, jja, jjx, ldx,
162 $ mycol, myrow, np, npcol, nprow, nq
163 DOUBLE PRECISION DIVISOR, EPS, RESID1
167 EXTERNAL blacs_gridinfo, dgebr2d, dgebs2d,
169 $
pzmatgen, zgamx2d, zgemm, zgsum2d,
173 INTEGER IZAMAX, NUMROC
174 DOUBLE PRECISION PDLAMCH
175 EXTERNAL izamax, numroc, pdlamch
178 INTRINSIC abs, dble,
max,
min, mod
184 ictxt = desca( ctxt_ )
185 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
187 eps = pdlamch( ictxt,
'eps' )
189 divisor = anorm * eps * dble( n )
191 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
193 CALL infog2l( ix, jx, descx, nprow, npcol, myrow, mycol, iix, jjx,
195 iroff = mod( ia-1, desca( mb_ ) )
196 icoff = mod( ja-1, desca( nb_ ) )
197 np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
198 nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
202 ipx = ipb + np * descx( nb_ )
203 ipa = ipx + nq * descx( nb_ )
214 DO 40 j = 1, nrhs, descx( nb_ )
215 jbrhs =
min( descx( nb_ ), nrhs-j+1 )
219 ioffx = iix + ( jjx - 1 ) * descx( lld_ )
220 CALL pbztran( ictxt,
'Column',
'Transpose', n, jbrhs,
221 $ descx( mb_ ), x( ioffx ), descx( lld_ ), zero,
222 $ work( ipx ), jbrhs, ixrow, icurcol, -1, iacol,
227 IF( mycol.EQ.icurcol )
THEN
228 CALL pzmatgen( ictxt,
'N',
'N', descx( m_ ), descx( n_ ),
229 $ descx( mb_ ), descx( nb_ ), work( ipb ), ldx,
230 $ ixrow, ixcol, ibseed, iix-1, np, jjx-1,
231 $ jbrhs, myrow, mycol, nprow, npcol )
238 DO 10 ii = iia, iia+np-1, desca( mb_ )
239 ib =
min( desca( mb_ ), iia+np-ii )
243 CALL pzmatgen( ictxt, symm, diag, desca( m_ ),
244 $ desca( n_ ), desca( mb_ ), desca( nb_ ),
245 $ work( ipa ), ib, desca( rsrc_ ),
246 $ desca( csrc_ ), iaseed, ii-1, ib,
247 $ jja-1, nq, myrow, mycol, nprow, npcol )
251 CALL zgemm(
'No transpose',
'Transpose', ib, jbrhs, nq,
252 $ -one, work( ipa ), ib, work( ipx ), jbrhs,
253 $ beta, work( ipb+ii-iia ), ldx )
257 ELSE IF( mycol.NE.icurcol )
THEN
259 CALL zlaset(
'All', np, jbrhs, zero, zero, work( ipb ),
266 CALL zgsum2d( ictxt,
'Row',
' ', np, jbrhs, work( ipb ), ldx,
269 IF( mycol.EQ.icurcol )
THEN
274 DO 20 jj = 0, jbrhs - 1
276 ii = izamax( np, work( ipb+jj*ldx ), 1 )
277 work( ipa+jj ) = abs( work( ipb+ii-1+jj*ldx ) )
278 work( ipw+jj ) = abs( x( ioffx + izamax( np,
279 $ x( ioffx + jj*descx( lld_ ) ), 1 )-1+jj*
282 work( ipa+jj ) = zero
283 work( ipw+jj ) = zero
291 CALL zgamx2d( ictxt,
'Column',
' ', 1, 2*jbrhs,
292 $ work( ipa ), 1, idumm, idumm, -1, 0, icurcol )
296 IF( myrow.EQ.0 )
THEN
297 DO 30 jj = 0, jbrhs - 1
298 resid1 = dble( work( ipa+jj ) ) /
299 $ ( dble( work( ipw+jj ) )*divisor )
300 IF( resid.LT.resid1 )
304 $
CALL dgesd2d( ictxt, 1, 1, resid, 1, 0, 0 )
307 ELSE IF( myrow.EQ.0 .AND. mycol.EQ.0 )
THEN
309 CALL dgerv2d( ictxt, 1, 1, resid1, 1, 0, icurcol )
310 IF( resid.LT.resid1 )
315 IF( mycol.EQ.icurcol )
317 icurcol = mod( icurcol+1, npcol )
321 IF( myrow.EQ.0 .AND. mycol.EQ.0 )
THEN
322 CALL dgebs2d( ictxt,
'All',
' ', 1, 1, resid, 1 )
324 CALL dgebr2d( ictxt,
'All',
' ', 1, 1, resid, 1, 0, 0 )