3
4
5
6
7
8
9
10 CHARACTER TRANS
11 INTEGER IA, IB, INFO, JA, JB, LWORK, M, N, NRHS
12
13
14 INTEGER DESCA( * ), DESCB( * )
15 COMPLEX*16 A( * ), B( * ), WORK( * )
16
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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
230 $ LLD_, MB_, M_, NB_, N_, RSRC_
231 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
232 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
233 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
234 DOUBLE PRECISION ZERO, ONE
235 parameter( zero = 0.0d+0, one = 1.0d+0 )
236 COMPLEX*16 CZERO, CONE
237 parameter( czero = ( 0.0d+0, 0.0d+0 ),
238 $ cone = ( 1.0d+0, 0.0d+0 ) )
239
240
241 LOGICAL LQUERY, TPSD
242 INTEGER BROW, IACOL, IAROW, IASCL, IBCOL, IBROW, IBSCL,
243 $ ICOFFA, ICOFFB, ICTXT, IPW, IROFFA, IROFFB,
244 $ LCM, LCMP, LTAU, LWF, LWMIN, LWS, MPA0, MPB0,
245 $ MYCOL, MYROW, NPB0, NPCOL, NPROW, NQA0,
246 $ NRHSQB0, SCLLEN
247 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
248
249
250 INTEGER IDUM1( 2 ), IDUM2( 2 )
251 DOUBLE PRECISION RWORK( 1 )
252
253
254 LOGICAL LSAME
255 INTEGER ILCM
256 INTEGER INDXG2P, NUMROC
257 DOUBLE PRECISION PDLAMCH, PZLANGE
260
261
265
266
267 INTRINSIC dble, dcmplx, ichar,
max,
min, mod
268
269
270
271
272
273 ictxt = desca( ctxt_ )
274 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
275
276
277
278 info = 0
279 IF( nprow.EQ.-1 ) THEN
280 info = -( 800 + ctxt_ )
281 ELSE
282 CALL chk1mat( m, 2, n, 3, ia, ja, desca, 8, info )
283 IF ( m .GE. n ) THEN
284 CALL chk1mat( m, 2, nrhs, 4, ib, jb, descb, 12, info )
285 ELSE
286 CALL chk1mat( n, 3, nrhs, 4, ib, jb, descb, 12, info )
287 ENDIF
288 IF( info.EQ.0 ) THEN
289 iroffa = mod( ia-1, desca( mb_ ) )
290 icoffa = mod( ja-1, desca( nb_ ) )
291 iroffb = mod( ib-1, descb( mb_ ) )
292 icoffb = mod( jb-1, descb( nb_ ) )
293 iarow =
indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
294 $ nprow )
295 iacol =
indxg2p( ia, desca( nb_ ), mycol, desca( csrc_ ),
296 $ npcol )
297 mpa0 =
numroc( m+iroffa, desca( mb_ ), myrow, iarow, nprow )
298 nqa0 =
numroc( n+icoffa, desca( nb_ ), mycol, iacol, npcol )
299
300 ibrow =
indxg2p( ib, descb( mb_ ), myrow, descb( rsrc_ ),
301 $ nprow )
302 ibcol =
indxg2p( ib, descb( nb_ ), mycol, descb( csrc_ ),
303 $ npcol )
304 nrhsqb0 =
numroc( nrhs+icoffb, descb( nb_ ), mycol, ibcol,
305 $ npcol )
306 IF( m.GE.n ) THEN
307 mpb0 =
numroc( m+iroffb, descb( mb_ ), myrow, ibrow,
308 $ nprow )
309 ltau =
numroc( ja+
min(m,n)-1, desca( nb_ ), mycol,
310 $ desca( csrc_ ), npcol )
311 lwf = desca( nb_ ) * ( mpa0 + nqa0 + desca( nb_ ) )
312 lws =
max( ( desca( nb_ )*( desca( nb_ ) - 1 ) ) / 2,
313 $ ( mpb0 + nrhsqb0 ) * desca( nb_ ) ) +
314 $ desca( nb_ )*desca( nb_ )
315 ELSE
316 lcm =
ilcm( nprow, npcol )
317 lcmp = lcm / nprow
318 npb0 =
numroc( n+iroffb, descb( mb_ ), myrow, ibrow,
319 $ nprow )
320 ltau =
numroc( ia+
min(m,n)-1, desca( mb_ ), myrow,
321 $ desca( rsrc_ ), nprow )
322 lwf = desca( mb_ ) * ( mpa0 + nqa0 + desca( mb_ ) )
323 lws =
max( ( desca( mb_ )*( desca( mb_ ) - 1 ) ) / 2,
325 $ desca( mb_ ), 0, 0, nprow ), desca( mb_ ), 0, 0,
326 $ lcmp ), nrhsqb0 ) )*desca( mb_ ) ) +
327 $ desca( mb_ ) * desca( mb_ )
328 END IF
329 lwmin = ltau +
max( lwf, lws )
330 work( 1 ) = dcmplx( dble( lwmin ) )
331 lquery = ( lwork.EQ.-1 )
332
333 tpsd = .true.
334 IF(
lsame( trans,
'N' ) )
335 $ tpsd = .false.
336
337 IF( .NOT.(
lsame( trans,
'N' ) .OR.
338 $
lsame( trans,
'C' ) ) )
THEN
339 info = -1
340 ELSE IF( m.LT.0 ) THEN
341 info = -2
342 ELSE IF( n.LT.0 ) THEN
343 info = -3
344 ELSE IF( nrhs.LT.0 ) THEN
345 info = -4
346 ELSE IF( m.GE.n .AND. iroffa.NE.iroffb ) THEN
347 info = -10
348 ELSE IF( m.GE.n .AND. iarow.NE.ibrow ) THEN
349 info = -10
350 ELSE IF( m.LT.n .AND. icoffa.NE.iroffb ) THEN
351 info = -10
352 ELSE IF( m.GE.n .AND. desca( mb_ ).NE.descb( mb_ ) ) THEN
353 info = -( 1200 + mb_ )
354 ELSE IF( m.LT.n .AND. desca( nb_ ).NE.descb( mb_ ) ) THEN
355 info = -( 1200 + mb_ )
356 ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
357 info = -( 1200 + ctxt_ )
358 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
359 info = -14
360 END IF
361 END IF
362
363 IF( .NOT.tpsd ) THEN
364 idum1( 1 ) = ichar( 'N' )
365 ELSE
366 idum1( 1 ) = ichar( 'C' )
367 END IF
368 idum2( 1 ) = 1
369 IF( lwork.EQ.-1 ) THEN
370 idum1( 2 ) = -1
371 ELSE
372 idum1( 2 ) = 1
373 END IF
374 idum2( 2 ) = 14
375 CALL pchk2mat( m, 2, n, 3, ia, ja, desca, 8, n, 3, nrhs, 4,
376 $ ib, jb, descb, 12, 2, idum1, idum2, info )
377 END IF
378
379 IF( info.NE.0 ) THEN
380 CALL pxerbla( ictxt,
'PZGELS', -info )
381 RETURN
382 ELSE IF( lquery ) THEN
383 RETURN
384 END IF
385
386
387
388 IF(
min( m, n, nrhs ).EQ.0 )
THEN
389 CALL pzlaset(
'Full',
max( m, n ), nrhs, czero, czero, b,
390 $ ib, jb, descb )
391 RETURN
392 END IF
393
394
395
397 smlnum = smlnum /
pdlamch( ictxt,
'P' )
398 bignum = one / smlnum
399 CALL pdlabad( ictxt, smlnum, bignum )
400
401
402
403 anrm =
pzlange(
'M', m, n, a, ia, ja, desca, rwork )
404 iascl = 0
405 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
406
407
408
409 CALL pzlascl(
'G', anrm, smlnum, m, n, a, ia, ja, desca,
410 $ info )
411 iascl = 1
412 ELSE IF( anrm.GT.bignum ) THEN
413
414
415
416 CALL pzlascl(
'G', anrm, bignum, m, n, a, ia, ja, desca,
417 $ info )
418 iascl = 2
419 ELSE IF( anrm.EQ.zero ) THEN
420
421
422
423 CALL pzlaset(
'F',
max( m, n ), nrhs, czero, czero, b, ib,
424 $ jb, descb )
425 GO TO 10
426 END IF
427
428 brow = m
429 IF( tpsd )
430 $ brow = n
431
432 bnrm =
pzlange(
'M', brow, nrhs, b, ib, jb, descb, rwork )
433
434 ibscl = 0
435 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
436
437
438
439 CALL pzlascl(
'G', bnrm, smlnum, brow, nrhs, b, ib, jb,
440 $ descb, info )
441 ibscl = 1
442 ELSE IF( bnrm.GT.bignum ) THEN
443
444
445
446 CALL pzlascl(
'G', bnrm, bignum, brow, nrhs, b, ib, jb,
447 $ descb, info )
448 ibscl = 2
449 END IF
450
451 ipw = ltau + 1
452
453 IF( m.GE.n ) THEN
454
455
456
457 CALL pzgeqrf( m, n, a, ia, ja, desca, work, work( ipw ),
458 $ lwork-ltau, info )
459
460
461
462 IF( .NOT.tpsd ) THEN
463
464
465
466
467
468 CALL pzunmqr(
'Left',
'Conjugate transpose', m, nrhs, n, a,
469 $ ia, ja, desca, work, b, ib, jb, descb,
470 $ work( ipw ), lwork-ltau, info )
471
472
473
474
475
476
477 CALL pztrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', n,
478 $ nrhs, cone, a, ia, ja, desca, b, ib, jb,
479 $ descb )
480
481 scllen = n
482
483 ELSE
484
485
486
487
488
489 CALL pztrsm( 'Left', 'Upper', 'Conjugate transpose',
490 $ 'Non-unit', n, nrhs, cone, a, ia, ja, desca,
491 $ b, ib, jb, descb )
492
493
494
495 CALL pzlaset(
'All', m-n, nrhs, czero, czero, b, ib+n, jb,
496 $ descb )
497
498
499
500
501 CALL pzunmqr(
'Left',
'No transpose', m, nrhs, n, a, ia, ja,
502 $ desca, work, b, ib, jb, descb, work( ipw ),
503 $ lwork-ltau, info )
504
505
506
507 scllen = m
508
509 END IF
510
511 ELSE
512
513
514
515 CALL pzgelqf( m, n, a, ia, ja, desca, work, work( ipw ),
516 $ lwork-ltau, info )
517
518
519
520 IF( .NOT.tpsd ) THEN
521
522
523
524
525
526
527 CALL pztrsm( 'Left', 'Lower', 'No transpose', 'Non-unit', m,
528 $ nrhs, cone, a, ia, ja, desca, b, ib, jb,
529 $ descb )
530
531
532
533 CALL pzlaset(
'All', n-m, nrhs, czero, czero, b, ib+m, jb,
534 $ descb )
535
536
537
538
539 CALL pzunmlq(
'Left',
'Conjugate transpose', n, nrhs, m, a,
540 $ ia, ja, desca, work, b, ib, jb, descb,
541 $ work( ipw ), lwork-ltau, info )
542
543
544
545 scllen = n
546
547 ELSE
548
549
550
551
552
553 CALL pzunmlq(
'Left',
'No transpose', n, nrhs, m, a, ia, ja,
554 $ desca, work, b, ib, jb, descb, work( ipw ),
555 $ lwork-ltau, info )
556
557
558
559
560
561
562 CALL pztrsm( 'Left', 'Lower', 'Conjugate transpose',
563 $ 'Non-unit', m, nrhs, cone, a, ia, ja, desca,
564 $ b, ib, jb, descb )
565
566 scllen = m
567
568 END IF
569
570 END IF
571
572
573
574 IF( iascl.EQ.1 ) THEN
575 CALL pzlascl(
'G', anrm, smlnum, scllen, nrhs, b, ib, jb,
576 $ descb, info )
577 ELSE IF( iascl.EQ.2 ) THEN
578 CALL pzlascl(
'G', anrm, bignum, scllen, nrhs, b, ib, jb,
579 $ descb, info )
580 END IF
581 IF( ibscl.EQ.1 ) THEN
582 CALL pzlascl(
'G', smlnum, bnrm, scllen, nrhs, b, ib, jb,
583 $ descb, info )
584 ELSE IF( ibscl.EQ.2 ) THEN
585 CALL pzlascl(
'G', bignum, bnrm, scllen, nrhs, b, ib, jb,
586 $ descb, info )
587 END IF
588
589 10 CONTINUE
590
591 work( 1 ) = dcmplx( dble( lwmin ) )
592
593 RETURN
594
595
596
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
integer function ilcm(m, n)
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
integer function numroc(n, nb, iproc, isrcproc, nprocs)
subroutine pchk2mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, mb, mbpos0, nb, nbpos0, ib, jb, descb, descbpos0, nextra, ex, expos, info)
double precision function pdlamch(ictxt, cmach)
subroutine pdlabad(ictxt, small, large)
subroutine pxerbla(ictxt, srname, info)
subroutine pzlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
subroutine pzgelqf(m, n, a, ia, ja, desca, tau, work, lwork, info)
subroutine pzgeqrf(m, n, a, ia, ja, desca, tau, work, lwork, info)
double precision function pzlange(norm, m, n, a, ia, ja, desca, work)
subroutine pzlascl(type, cfrom, cto, m, n, a, ia, ja, desca, info)
subroutine pzunmlq(side, trans, m, n, k, a, ia, ja, desca, tau, c, ic, jc, descc, work, lwork, info)
subroutine pzunmqr(side, trans, m, n, k, a, ia, ja, desca, tau, c, ic, jc, descc, work, lwork, info)