ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
psscaex.f
Go to the documentation of this file.
1  PROGRAM psscaex
2 *
3 * -- ScaLAPACK example code --
4 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5 * and University of California, Berkeley.
6 *
7 * Written by Antoine Petitet, August 1995 (petitet@cs.utk.edu)
8 *
9 * This program solves a linear system by calling the ScaLAPACK
10 * routine PSGESV. The input matrix and right-and-sides are
11 * read from a file. The solution is written to a file.
12 *
13 * .. Parameters ..
14  INTEGER realsz, intgsz, memsiz, totmem
15  parameter( realsz = 4, intgsz = 4, totmem = 2000000,
16  $ memsiz = totmem / realsz )
17  INTEGER block_cyclic_2d, csrc_, ctxt_, dlen_, dt_,
18  $ lld_, mb_, m_, nb_, n_, rsrc_
19  parameter( block_cyclic_2d = 1, dlen_ = 9, dt_ = 1,
20  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
21  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
22  REAL one
23  parameter( one = 1.0e+0 )
24 * ..
25 * .. Local Scalars ..
26  CHARACTER*80 outfile
27  INTEGER iam, ictxt, info, ipa, ipacpy, ipb, ippiv, ipx,
28  $ ipw, lipiv, mycol, myrow, n, nb, nout, npcol,
29  $ nprocs, nprow, np, nq, nqrhs, nrhs, worksiz
30  REAL anorm, bnorm, eps, xnorm, resid
31 * ..
32 * .. Local Arrays ..
33  INTEGER desca( dlen_ ), descb( dlen_ ), descx( dlen_ )
34  REAL mem( memsiz )
35 * ..
36 * .. External Subroutines ..
37  EXTERNAL blacs_exit, blacs_get, blacs_gridexit,
38  $ blacs_gridinfo, blacs_gridinit, blacs_pinfo,
39  $ descinit, igsum2d, pdscaexinfo, psgesv,
40  $ psgemm, pslacpy, pslaprnt, pslaread, pslawrite
41 * ..
42 * .. External Functions ..
43  INTEGER iceil, numroc
44  REAL pslamch, pslange
45  EXTERNAL iceil, numroc, pslamch, pslange
46 * ..
47 * .. Intrinsic Functions ..
48  INTRINSIC dble, max
49 * ..
50 * .. Executable Statements ..
51 *
52 * Get starting information
53 *
54  CALL blacs_pinfo( iam, nprocs )
55  CALL pdscaexinfo( outfile, nout, n, nrhs, nb, nprow, npcol, mem,
56  $ iam, nprocs )
57 *
58 * Define process grid
59 *
60  CALL blacs_get( -1, 0, ictxt )
61  CALL blacs_gridinit( ictxt, 'Row-major', nprow, npcol )
62  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
63 *
64 * Go to bottom of process grid loop if this case doesn't use my
65 * process
66 *
67  IF( myrow.GE.nprow .OR. mycol.GE.npcol )
68  $ GO TO 20
69 *
70  np = numroc( n, nb, myrow, 0, nprow )
71  nq = numroc( n, nb, mycol, 0, npcol )
72  nqrhs = numroc( nrhs, nb, mycol, 0, npcol )
73 *
74 * Initialize the array descriptor for the matrix A and B
75 *
76  CALL descinit( desca, n, n, nb, nb, 0, 0, ictxt, max( 1, np ),
77  $ info )
78  CALL descinit( descb, n, nrhs, nb, nb, 0, 0, ictxt, max( 1, np ),
79  $ info )
80  CALL descinit( descx, n, nrhs, nb, nb, 0, 0, ictxt, max( 1, np ),
81  $ info )
82 *
83 * Assign pointers into MEM for SCALAPACK arrays, A is
84 * allocated starting at position MEM( 1 )
85 *
86  ipa = 1
87  ipacpy = ipa + desca( lld_ )*nq
88  ipb = ipacpy + desca( lld_ )*nq
89  ipx = ipb + descb( lld_ )*nqrhs
90  ippiv = ipx + descb( lld_ )*nqrhs
91  lipiv = iceil( intgsz*( np+nb ), realsz )
92  ipw = ippiv + max( np, lipiv )
93 *
94  worksiz = nb
95 *
96 * Check for adequate memory for problem size
97 *
98  info = 0
99  IF( ipw+worksiz.GT.memsiz ) THEN
100  IF( iam.EQ.0 )
101  $ WRITE( nout, fmt = 9998 ) 'test', ( ipw+worksiz )*realsz
102  info = 1
103  END IF
104 *
105 * Check all processes for an error
106 *
107  CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, 0 )
108  IF( info.GT.0 ) THEN
109  IF( iam.EQ.0 )
110  $ WRITE( nout, fmt = 9999 ) 'MEMORY'
111  GO TO 10
112  END IF
113 *
114 * Read from file and distribute matrices A and B
115 *
116  CALL pslaread( 'SSCAEXMAT.dat', mem( ipa ), desca, 0, 0,
117  $ mem( ipw ) )
118  CALL pslaread( 'SSCAEXRHS.dat', mem( ipb ), descb, 0, 0,
119  $ mem( ipw ) )
120 *
121 * Make a copy of A and the rhs for checking purposes
122 *
123  CALL pslacpy( 'All', n, n, mem( ipa ), 1, 1, desca,
124  $ mem( ipacpy ), 1, 1, desca )
125  CALL pslacpy( 'All', n, nrhs, mem( ipb ), 1, 1, descb,
126  $ mem( ipx ), 1, 1, descx )
127 *
128 **********************************************************************
129 * Call ScaLAPACK PSGESV routine
130 **********************************************************************
131 *
132  IF( iam.EQ.0 ) THEN
133  WRITE( nout, fmt = * )
134  WRITE( nout, fmt = * )
135  $ '***********************************************'
136  WRITE( nout, fmt = * )
137  $ 'Example of ScaLAPACK routine call: (PSGESV)'
138  WRITE( nout, fmt = * )
139  $ '***********************************************'
140  WRITE( nout, fmt = * )
141  WRITE( nout, fmt = * ) 'A * X = B, Matrix A:'
142  WRITE( nout, fmt = * )
143  END IF
144  CALL pslaprnt( n, n, mem( ipa ), 1, 1, desca, 0, 0,
145  $ 'A', nout, mem( ipw ) )
146  IF( iam.EQ.0 ) THEN
147  WRITE( nout, fmt = * )
148  WRITE( nout, fmt = * ) 'Matrix B:'
149  WRITE( nout, fmt = * )
150  END IF
151  CALL pslaprnt( n, nrhs, mem( ipb ), 1, 1, descb, 0, 0,
152  $ 'B', nout, mem( ipw ) )
153 *
154  CALL psgesv( n, nrhs, mem( ipa ), 1, 1, desca, mem( ippiv ),
155  $ mem( ipb ), 1, 1, descb, info )
156 *
157  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
158  WRITE( nout, fmt = * )
159  WRITE( nout, fmt = * ) 'INFO code returned by PSGESV = ', info
160  WRITE( nout, fmt = * )
161  WRITE( nout, fmt = * ) 'Matrix X = A^{-1} * B'
162  WRITE( nout, fmt = * )
163  END IF
164  CALL pslaprnt( n, nrhs, mem( ipb ), 1, 1, descb, 0, 0, 'X', nout,
165  $ mem( ipw ) )
166  CALL pslawrite( 'SSCAEXSOL.dat', n, nrhs, mem( ipb ), 1, 1, descb,
167  $ 0, 0, mem( ipw ) )
168 *
169 * Compute residual ||A * X - B|| / ( ||X|| * ||A|| * eps * N )
170 *
171  eps = pslamch( ictxt, 'Epsilon' )
172  anorm = pslange( 'I', n, n, mem( ipa ), 1, 1, desca, mem( ipw ) )
173  bnorm = pslange( 'I', n, nrhs, mem( ipb ), 1, 1, descb,
174  $ mem( ipw ) )
175  CALL psgemm( 'No transpose', 'No transpose', n, nrhs, n, one,
176  $ mem( ipacpy ), 1, 1, desca, mem( ipb ), 1, 1, descb,
177  $ -one, mem( ipx ), 1, 1, descx )
178  xnorm = pslange( 'I', n, nrhs, mem( ipx ), 1, 1, descx,
179  $ mem( ipw ) )
180  resid = xnorm / ( anorm * bnorm * eps * dble( n ) )
181 *
182  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
183  WRITE( nout, fmt = * )
184  WRITE( nout, fmt = * )
185  $ '||A * X - B|| / ( ||X|| * ||A|| * eps * N ) = ', resid
186  WRITE( nout, fmt = * )
187  IF( resid.LT.10.0d+0 ) THEN
188  WRITE( nout, fmt = * ) 'The answer is correct.'
189  ELSE
190  WRITE( nout, fmt = * ) 'The answer is suspicious.'
191  END IF
192  END IF
193 *
194  10 CONTINUE
195 *
196  CALL blacs_gridexit( ictxt )
197 *
198  20 CONTINUE
199 *
200 * Print ending messages and close output file
201 *
202  IF( iam.EQ.0 ) THEN
203  WRITE( nout, fmt = * )
204  WRITE( nout, fmt = * )
205  WRITE( nout, fmt = 9997 )
206  WRITE( nout, fmt = * )
207  IF( nout.NE.6 .AND. nout.NE.0 )
208  $ CLOSE ( nout )
209  END IF
210 *
211  CALL blacs_exit( 0 )
212 *
213  9999 FORMAT( 'Bad ', a6, ' parameters: going on to next test case.' )
214  9998 FORMAT( 'Unable to perform ', a, ': need TOTMEM of at least',
215  $ i11 )
216  9997 FORMAT( 'END OF TESTS.' )
217 *
218  stop
219 *
220 * End of PSSCAEX
221 *
222  END
pslamch
real function pslamch(ICTXT, CMACH)
Definition: pcblastst.f:7455
max
#define max(A, B)
Definition: pcgemr.c:180
pslange
real function pslange(NORM, M, N, A, IA, JA, DESCA, WORK)
Definition: pslange.f:3
pslaprnt
subroutine pslaprnt(M, N, A, IA, JA, DESCA, IRPRNT, ICPRNT, CMATNM, NOUT, WORK)
Definition: pslaprnt.f:3
pslaread
subroutine pslaread(FILNAM, A, DESCA, IRREAD, ICREAD, WORK)
Definition: pslaread.f:2
pslacpy
subroutine pslacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pslacpy.f:3
pslawrite
subroutine pslawrite(FILNAM, M, N, A, IA, JA, DESCA, IRWRIT, ICWRIT, WORK)
Definition: pslawrite.f:3
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
psgesv
subroutine psgesv(N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, INFO)
Definition: psgesv.f:3
pdscaexinfo
subroutine pdscaexinfo(SUMMRY, NOUT, N, NRHS, NB, NPROW, NPCOL, WORK, IAM, NPROCS)
Definition: pdscaexinfo.f:3
psscaex
program psscaex
Definition: psscaex.f:1
iceil
integer function iceil(INUM, IDENOM)
Definition: iceil.f:2