1 SUBROUTINE pslaschk( 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
16 INTEGER DESCA( * ), DESCX( * )
17 REAL 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.0e+0, zero = 0.0e+0 )
158 INTEGER IACOL, IAROW, IB, ICOFF, ICTXT, ICURCOL, IDUMM,
159 $ II, IIA, IIX, IOFFX, IPA, IPB, IPW, IPX, IROFF,
160 $ ixcol, ixrow, j, jbrhs, jj, jja, jjx, ldx,
161 $ mycol, myrow, np, npcol, nprow, nq
162 REAL BETA, DIVISOR, EPS, RESID1
166 $ sgamx2d, sgebr2d, sgebs2d, sgemm,
167 $ sgerv2d, sgesd2d, sgsum2d, slaset
170 INTEGER ISAMAX, NUMROC
172 EXTERNAL isamax, numroc, pslamch
175 INTRINSIC abs,
max,
min, mod, real
181 ictxt = desca( ctxt_ )
182 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
184 eps = pslamch( ictxt,
'eps' )
186 divisor = anorm * eps * real( n )
188 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
190 CALL infog2l( ix, jx, descx, nprow, npcol, myrow, mycol, iix, jjx,
192 iroff = mod( ia-1, desca( mb_ ) )
193 icoff = mod( ja-1, desca( nb_ ) )
194 np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
195 nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
199 ipx = ipb + np * descx( nb_ )
200 ipa = ipx + nq * descx( nb_ )
211 DO 40 j = 1, nrhs, descx( nb_ )
212 jbrhs =
min( descx( nb_ ), nrhs-j+1 )
216 ioffx = iix + ( jjx - 1 ) * descx( lld_ )
217 CALL pbstran( ictxt,
'Column',
'Transpose', n, jbrhs,
218 $ descx( mb_ ), x( ioffx ), descx( lld_ ), zero,
219 $ work( ipx ), jbrhs, ixrow, icurcol, -1, iacol,
224 IF( mycol.EQ.icurcol )
THEN
225 CALL psmatgen( ictxt,
'N',
'N', descx( m_ ), descx( n_ ),
226 $ descx( mb_ ), descx( nb_ ), work( ipb ), ldx,
227 $ ixrow, ixcol, ibseed, iix-1, np, jjx-1,
228 $ jbrhs, myrow, mycol, nprow, npcol )
235 DO 10 ii = iia, iia+np-1, desca( mb_ )
236 ib =
min( desca( mb_ ), iia+np-ii )
240 CALL psmatgen( ictxt, symm, diag, desca( m_ ),
241 $ desca( n_ ), desca( mb_ ), desca( nb_ ),
242 $ work( ipa ), ib, desca( rsrc_ ),
243 $ desca( csrc_ ), iaseed, ii-1, ib,
244 $ jja-1, nq, myrow, mycol, nprow, npcol )
248 CALL sgemm(
'No transpose',
'Transpose', ib, jbrhs, nq,
249 $ -one, work( ipa ), ib, work( ipx ), jbrhs,
250 $ beta, work( ipb+ii-iia ), ldx )
254 ELSE IF( mycol.NE.icurcol )
THEN
256 CALL slaset(
'All', np, jbrhs, zero, zero, work( ipb ),
263 CALL sgsum2d( ictxt,
'Row',
' ', np, jbrhs, work( ipb ), ldx,
266 IF( mycol.EQ.icurcol )
THEN
271 DO 20 jj = 0, jbrhs - 1
273 ii = isamax( np, work( ipb+jj*ldx ), 1 )
274 work( ipa+jj ) = abs( work( ipb+ii-1+jj*ldx ) )
275 work( ipw+jj ) = abs( x( ioffx + isamax( np,
276 $ x( ioffx + jj*descx( lld_ ) ), 1 )-1+jj*
279 work( ipa+jj ) = zero
280 work( ipw+jj ) = zero
288 CALL sgamx2d( ictxt,
'Column',
' ', 1, 2*jbrhs,
289 $ work( ipa ), 1, idumm, idumm, -1, 0, icurcol )
293 IF( myrow.EQ.0 )
THEN
294 DO 30 jj = 0, jbrhs - 1
295 resid1 = work( ipa+jj ) / ( work( ipw+jj )*divisor )
296 IF( resid.LT.resid1 )
300 $
CALL sgesd2d( ictxt, 1, 1, resid, 1, 0, 0 )
303 ELSE IF( myrow.EQ.0 .AND. mycol.EQ.0 )
THEN
305 CALL sgerv2d( ictxt, 1, 1, resid1, 1, 0, icurcol )
306 IF( resid.LT.resid1 )
311 IF( mycol.EQ.icurcol )
313 icurcol = mod( icurcol+1, npcol )
317 IF( myrow.EQ.0 .AND. mycol.EQ.0 )
THEN
318 CALL sgebs2d( ictxt,
'All',
' ', 1, 1, resid, 1 )
320 CALL sgebr2d( ictxt,
'All',
' ', 1, 1, resid, 1, 0, 0 )