3
4
5
6
7
8
9
10 CHARACTER AFORM, DIAG
11 INTEGER IA, IASEED, JA, M, N
12 DOUBLE PRECISION ANORM, FRESID
13
14
15 INTEGER DESCA( * )
16 COMPLEX*16 A( * ), WORK( * )
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
149 $ LLD_, MB_, M_, NB_, N_, RSRC_
150 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
151 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
152 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
153 COMPLEX*16 ONE
154 parameter( one = (1.0d+0, 0.0d+0) )
155
156
157 INTEGER IACOL, IAROW, ICOFF, ICTXT, ICURCOL, ICURROW,
158 $ II, IIA, IOFFA, IROFF, JB, JJ, JJA, JN, KK,
159 $ LDA, LDW, LDWP1, MP, MYCOL, MYROW, NPCOL,
160 $ NPROW, NQ
161 DOUBLE PRECISION EPS
162
163
165
166
167 LOGICAL LSAME
168 INTEGER ICEIL, NUMROC
169 DOUBLE PRECISION PDLAMCH, PZLANGE
171
172
173 INTRINSIC dble,
max,
min, mod
174
175
176
177 ictxt = desca( ctxt_ )
178 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
180 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
181 $ iarow, iacol )
182
183
184
185 iroff = mod( ia-1, desca( mb_ ) )
186 icoff = mod( ja-1, desca( nb_ ) )
187 mp =
numroc( m+iroff, desca( mb_ ), myrow, iarow, nprow )
188 nq =
numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
189 IF( myrow.EQ.iarow )
190 $ mp = mp-iroff
191 IF( mycol.EQ.iacol )
192 $ nq = nq-icoff
193 jn =
min(
iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
195 ldwp1 = ldw + 1
196 lda = desca( lld_ )
197 ioffa = iia + ( jja - 1 )*lda
198
199 IF(
lsame( aform,
'H' ) )
THEN
200
201
202
203 ii = 1
204 icurrow = iarow
205 icurcol = iacol
206 jb = jn - ja + 1
207
208 IF( mycol.EQ.icurcol ) THEN
209 CALL pzmatgen( ictxt, aform, diag, desca( m_ ), desca( n_ ),
210 $ desca( mb_ ), desca( nb_ ), work, ldw,
211 $ desca( rsrc_ ), desca( csrc_ ), iaseed,
212 $ iia-1, mp, jja-1, jb, myrow, mycol, nprow,
213 $ npcol )
214 IF( myrow.EQ.icurrow ) THEN
215 DO 10, kk = 0, jb-1
216 work( ii+kk*ldwp1 ) = dble( work( ii+kk*ldwp1 ) )
217 10 CONTINUE
218 END IF
219 CALL zmatadd( mp, jb, -one, work, ldw, one, a( ioffa ),
220 $ lda )
221 jja = jja + jb
222 ioffa = ioffa + jb*lda
223 END IF
224
225 IF( myrow.EQ.icurrow )
226 $ ii = ii + jb
227 icurrow = mod( icurrow+1, nprow )
228 icurcol = mod( icurcol+1, npcol )
229
230 DO 30, jj = jn+1, ja+n-1, desca( nb_ )
231 jb =
min( ja+n-jj, desca( nb_ ) )
232
233 IF( mycol.EQ.icurcol ) THEN
234 CALL pzmatgen( ictxt, aform, diag, desca( m_ ),
235 $ desca( n_ ), desca( mb_ ), desca( nb_ ),
236 $ work, ldw, desca( rsrc_ ), desca( csrc_ ),
237 $ iaseed, iia-1, mp, jja-1, jb, myrow,
238 $ mycol, nprow, npcol )
239 IF( myrow.EQ.icurrow ) THEN
240 DO 20, kk = 0, jb-1
241 work( ii+kk*ldwp1 ) = dble( work( ii+kk*ldwp1 ) )
242 20 CONTINUE
243 END IF
244 CALL zmatadd( mp, jb, -one, work, ldw, one, a( ioffa ),
245 $ lda )
246 jja = jja + jb
247 ioffa = ioffa + jb*lda
248 END IF
249 IF( myrow.EQ.icurrow )
250 $ ii = ii + jb
251 icurrow = mod( icurrow+1, nprow )
252 icurcol = mod( icurcol+1, npcol )
253 30 CONTINUE
254
255 ELSE
256
257
258
259 IF( mycol.EQ.iacol ) THEN
260 jb = jn-ja+1
261 CALL pzmatgen( ictxt, aform, diag, desca( m_ ), desca( n_ ),
262 $ desca( mb_ ), desca( nb_ ), work, ldw,
263 $ desca( rsrc_ ), desca( csrc_ ), iaseed,
264 $ iia-1, mp, jja-1, jb, myrow, mycol, nprow,
265 $ npcol )
266 CALL zmatadd( mp, jb, -one, work, ldw, one, a( ioffa ),
267 $ lda )
268 jja = jja + jb
269 nq = nq - jb
270 ioffa = ioffa + jb * lda
271 END IF
272
273
274
275 DO 40 jj = jja, jja+nq-1, desca( nb_ )
276 jb =
min( desca( nb_ ), jja+nq-jj )
277 ioffa = iia + ( jj - 1 )*lda
278 CALL pzmatgen( ictxt, aform, diag, desca( m_ ), desca( n_ ),
279 $ desca( mb_ ), desca( nb_ ), work, ldw,
280 $ desca( rsrc_ ), desca( csrc_ ), iaseed,
281 $ iia-1, mp, jj-1, jb, myrow, mycol, nprow,
282 $ npcol )
283 CALL zmatadd( mp, jb, -one, work, ldw, one, a( ioffa ),
284 $ lda )
285 40 CONTINUE
286
287 END IF
288
289
290
291 fresid =
pzlange(
'I', m, n, a, ia, ja, desca, work ) /
292 $ (
max( m, n ) * eps * anorm )
293
294 RETURN
295
296
297
subroutine pzmatgen(ictxt, aform, diag, m, n, mb, nb, a, lda, iarow, iacol, iseed, iroff, irnum, icoff, icnum, myrow, mycol, nprow, npcol)
integer function iceil(inum, idenom)
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
integer function numroc(n, nb, iproc, isrcproc, nprocs)
double precision function pdlamch(ictxt, cmach)
double precision function pzlange(norm, m, n, a, ia, ja, desca, work)
subroutine zmatadd(m, n, alpha, a, lda, beta, c, ldc)