131
132
133
134
135
136
137 CHARACTER DIAG, TRANS, UPLO
138 INTEGER IMAT, INFO, N
139
140
141 INTEGER ISEED( 4 )
142 REAL RWORK( * )
143 COMPLEX AP( * ), B( * ), WORK( * )
144
145
146
147
148
149 REAL ONE, TWO, ZERO
150 parameter( one = 1.0e+0, two = 2.0e+0, zero = 0.0e+0 )
151
152
153 LOGICAL UPPER
154 CHARACTER DIST, PACKIT, TYPE
155 CHARACTER*3 PATH
156 INTEGER I, IY, J, JC, JCNEXT, JCOUNT, JJ, JL, JR, JX,
157 $ KL, KU, MODE
158 REAL ANORM, BIGNUM, BNORM, BSCAL, C, CNDNUM, REXP,
159 $ SFAC, SMLNUM, T, TEXP, TLEFT, TSCAL, ULP, UNFL,
160 $ X, Y, Z
161 COMPLEX CTEMP, PLUS1, PLUS2, RA, RB, S, STAR1
162
163
164 LOGICAL LSAME
165 INTEGER ICAMAX
166 REAL SLAMCH
167 COMPLEX CLARND
169
170
173
174
175 INTRINSIC abs, cmplx, conjg, max, real, sqrt
176
177
178
179 path( 1: 1 ) = 'Complex precision'
180 path( 2: 3 ) = 'TP'
181 unfl =
slamch(
'Safe minimum' )
183 smlnum = unfl
184 bignum = ( one-ulp ) / smlnum
185 IF( ( imat.GE.7 .AND. imat.LE.10 ) .OR. imat.EQ.18 ) THEN
186 diag = 'U'
187 ELSE
188 diag = 'N'
189 END IF
190 info = 0
191
192
193
194 IF( n.LE.0 )
195 $ RETURN
196
197
198
199 upper =
lsame( uplo,
'U' )
200 IF( upper ) THEN
201 CALL clatb4( path, imat, n, n,
TYPE, KL, KU, ANORM, MODE,
202 $ CNDNUM, DIST )
203 packit = 'C'
204 ELSE
205 CALL clatb4( path, -imat, n, n,
TYPE, KL, KU, ANORM, MODE,
206 $ CNDNUM, DIST )
207 packit = 'R'
208 END IF
209
210
211
212 IF( imat.LE.6 ) THEN
213 CALL clatms( n, n, dist, iseed,
TYPE, RWORK, MODE, CNDNUM,
214 $ ANORM, KL, KU, PACKIT, AP, N, WORK, INFO )
215
216
217
218
219
220
221 ELSE IF( imat.EQ.7 ) THEN
222 IF( upper ) THEN
223 jc = 1
224 DO 20 j = 1, n
225 DO 10 i = 1, j - 1
226 ap( jc+i-1 ) = zero
227 10 CONTINUE
228 ap( jc+j-1 ) = j
229 jc = jc + j
230 20 CONTINUE
231 ELSE
232 jc = 1
233 DO 40 j = 1, n
234 ap( jc ) = j
235 DO 30 i = j + 1, n
236 ap( jc+i-j ) = zero
237 30 CONTINUE
238 jc = jc + n - j + 1
239 40 CONTINUE
240 END IF
241
242
243
244
245
246
247
248 ELSE IF( imat.LE.10 ) THEN
249 IF( upper ) THEN
250 jc = 0
251 DO 60 j = 1, n
252 DO 50 i = 1, j - 1
253 ap( jc+i ) = zero
254 50 CONTINUE
255 ap( jc+j ) = j
256 jc = jc + j
257 60 CONTINUE
258 ELSE
259 jc = 1
260 DO 80 j = 1, n
261 ap( jc ) = j
262 DO 70 i = j + 1, n
263 ap( jc+i-j ) = zero
264 70 CONTINUE
265 jc = jc + n - j + 1
266 80 CONTINUE
267 END IF
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327 star1 = 0.25*
clarnd( 5, iseed )
328 sfac = 0.5
329 plus1 = sfac*
clarnd( 5, iseed )
330 DO 90 j = 1, n, 2
331 plus2 = star1 / plus1
332 work( j ) = plus1
333 work( n+j ) = star1
334 IF( j+1.LE.n ) THEN
335 work( j+1 ) = plus2
336 work( n+j+1 ) = zero
337 plus1 = star1 / plus2
338 rexp = real(
clarnd( 2, iseed ) )
339 IF( rexp.LT.zero ) THEN
340 star1 = -sfac**( one-rexp )*
clarnd( 5, iseed )
341 ELSE
342 star1 = sfac**( one+rexp )*
clarnd( 5, iseed )
343 END IF
344 END IF
345 90 CONTINUE
346
347 x = sqrt( cndnum ) - one / sqrt( cndnum )
348 IF( n.GT.2 ) THEN
349 y = sqrt( two / real( n-2 ) )*x
350 ELSE
351 y = zero
352 END IF
353 z = x*x
354
355 IF( upper ) THEN
356
357
358
359
360 jc = 1
361 DO 100 j = 2, n
362 ap( jc+1 ) = y
363 IF( j.GT.2 )
364 $ ap( jc+j-1 ) = work( j-2 )
365 IF( j.GT.3 )
366 $ ap( jc+j-2 ) = work( n+j-3 )
367 jc = jc + j
368 100 CONTINUE
369 jc = jc - n
370 ap( jc+1 ) = z
371 DO 110 j = 2, n - 1
372 ap( jc+j ) = y
373 110 CONTINUE
374 ELSE
375
376
377
378
379 DO 120 i = 2, n - 1
380 ap( i ) = y
381 120 CONTINUE
382 ap( n ) = z
383 jc = n + 1
384 DO 130 j = 2, n - 1
385 ap( jc+1 ) = work( j-1 )
386 IF( j.LT.n-1 )
387 $ ap( jc+2 ) = work( n+j-1 )
388 ap( jc+n-j ) = y
389 jc = jc + n - j + 1
390 130 CONTINUE
391 END IF
392
393
394
395 IF( upper ) THEN
396 jc = 1
397 DO 150 j = 1, n - 1
398 jcnext = jc + j
399 ra = ap( jcnext+j-1 )
400 rb = two
401 CALL crotg( ra, rb, c, s )
402
403
404
405 IF( n.GT.j+1 ) THEN
406 jx = jcnext + j
407 DO 140 i = j + 2, n
408 ctemp = c*ap( jx+j ) + s*ap( jx+j+1 )
409 ap( jx+j+1 ) = -conjg( s )*ap( jx+j ) +
410 $ c*ap( jx+j+1 )
411 ap( jx+j ) = ctemp
412 jx = jx + i
413 140 CONTINUE
414 END IF
415
416
417
418 IF( j.GT.1 )
419 $
CALL crot( j-1, ap( jcnext ), 1, ap( jc ), 1, -c, -s )
420
421
422
423 ap( jcnext+j-1 ) = -ap( jcnext+j-1 )
424 jc = jcnext
425 150 CONTINUE
426 ELSE
427 jc = 1
428 DO 170 j = 1, n - 1
429 jcnext = jc + n - j + 1
430 ra = ap( jc+1 )
431 rb = two
432 CALL crotg( ra, rb, c, s )
433 s = conjg( s )
434
435
436
437 IF( n.GT.j+1 )
438 $
CALL crot( n-j-1, ap( jcnext+1 ), 1, ap( jc+2 ), 1, c,
439 $ -s )
440
441
442
443 IF( j.GT.1 ) THEN
444 jx = 1
445 DO 160 i = 1, j - 1
446 ctemp = -c*ap( jx+j-i ) + s*ap( jx+j-i+1 )
447 ap( jx+j-i+1 ) = -conjg( s )*ap( jx+j-i ) -
448 $ c*ap( jx+j-i+1 )
449 ap( jx+j-i ) = ctemp
450 jx = jx + n - i + 1
451 160 CONTINUE
452 END IF
453
454
455
456 ap( jc+1 ) = -ap( jc+1 )
457 jc = jcnext
458 170 CONTINUE
459 END IF
460
461
462
463
464
465 ELSE IF( imat.EQ.11 ) THEN
466
467
468
469
470
471 IF( upper ) THEN
472 jc = 1
473 DO 180 j = 1, n
474 CALL clarnv( 4, iseed, j-1, ap( jc ) )
475 ap( jc+j-1 ) =
clarnd( 5, iseed )*two
476 jc = jc + j
477 180 CONTINUE
478 ELSE
479 jc = 1
480 DO 190 j = 1, n
481 IF( j.LT.n )
482 $
CALL clarnv( 4, iseed, n-j, ap( jc+1 ) )
483 ap( jc ) =
clarnd( 5, iseed )*two
484 jc = jc + n - j + 1
485 190 CONTINUE
486 END IF
487
488
489
490 CALL clarnv( 2, iseed, n, b )
492 bnorm = abs( b( iy ) )
493 bscal = bignum / max( one, bnorm )
494 CALL csscal( n, bscal, b, 1 )
495
496 ELSE IF( imat.EQ.12 ) THEN
497
498
499
500
501
502 CALL clarnv( 2, iseed, n, b )
503 tscal = one / max( one, real( n-1 ) )
504 IF( upper ) THEN
505 jc = 1
506 DO 200 j = 1, n
507 CALL clarnv( 4, iseed, j-1, ap( jc ) )
508 CALL csscal( j-1, tscal, ap( jc ), 1 )
509 ap( jc+j-1 ) =
clarnd( 5, iseed )
510 jc = jc + j
511 200 CONTINUE
512 ap( n*( n+1 ) / 2 ) = smlnum*ap( n*( n+1 ) / 2 )
513 ELSE
514 jc = 1
515 DO 210 j = 1, n
516 CALL clarnv( 2, iseed, n-j, ap( jc+1 ) )
517 CALL csscal( n-j, tscal, ap( jc+1 ), 1 )
518 ap( jc ) =
clarnd( 5, iseed )
519 jc = jc + n - j + 1
520 210 CONTINUE
521 ap( 1 ) = smlnum*ap( 1 )
522 END IF
523
524 ELSE IF( imat.EQ.13 ) THEN
525
526
527
528
529
530 CALL clarnv( 2, iseed, n, b )
531 IF( upper ) THEN
532 jc = 1
533 DO 220 j = 1, n
534 CALL clarnv( 4, iseed, j-1, ap( jc ) )
535 ap( jc+j-1 ) =
clarnd( 5, iseed )
536 jc = jc + j
537 220 CONTINUE
538 ap( n*( n+1 ) / 2 ) = smlnum*ap( n*( n+1 ) / 2 )
539 ELSE
540 jc = 1
541 DO 230 j = 1, n
542 CALL clarnv( 4, iseed, n-j, ap( jc+1 ) )
543 ap( jc ) =
clarnd( 5, iseed )
544 jc = jc + n - j + 1
545 230 CONTINUE
546 ap( 1 ) = smlnum*ap( 1 )
547 END IF
548
549 ELSE IF( imat.EQ.14 ) THEN
550
551
552
553
554
555 IF( upper ) THEN
556 jcount = 1
557 jc = ( n-1 )*n / 2 + 1
558 DO 250 j = n, 1, -1
559 DO 240 i = 1, j - 1
560 ap( jc+i-1 ) = zero
561 240 CONTINUE
562 IF( jcount.LE.2 ) THEN
563 ap( jc+j-1 ) = smlnum*
clarnd( 5, iseed )
564 ELSE
565 ap( jc+j-1 ) =
clarnd( 5, iseed )
566 END IF
567 jcount = jcount + 1
568 IF( jcount.GT.4 )
569 $ jcount = 1
570 jc = jc - j + 1
571 250 CONTINUE
572 ELSE
573 jcount = 1
574 jc = 1
575 DO 270 j = 1, n
576 DO 260 i = j + 1, n
577 ap( jc+i-j ) = zero
578 260 CONTINUE
579 IF( jcount.LE.2 ) THEN
580 ap( jc ) = smlnum*
clarnd( 5, iseed )
581 ELSE
582 ap( jc ) =
clarnd( 5, iseed )
583 END IF
584 jcount = jcount + 1
585 IF( jcount.GT.4 )
586 $ jcount = 1
587 jc = jc + n - j + 1
588 270 CONTINUE
589 END IF
590
591
592
593 IF( upper ) THEN
594 b( 1 ) = zero
595 DO 280 i = n, 2, -2
596 b( i ) = zero
597 b( i-1 ) = smlnum*
clarnd( 5, iseed )
598 280 CONTINUE
599 ELSE
600 b( n ) = zero
601 DO 290 i = 1, n - 1, 2
602 b( i ) = zero
603 b( i+1 ) = smlnum*
clarnd( 5, iseed )
604 290 CONTINUE
605 END IF
606
607 ELSE IF( imat.EQ.15 ) THEN
608
609
610
611
612
613 texp = one / max( one, real( n-1 ) )
614 tscal = smlnum**texp
615 CALL clarnv( 4, iseed, n, b )
616 IF( upper ) THEN
617 jc = 1
618 DO 310 j = 1, n
619 DO 300 i = 1, j - 2
620 ap( jc+i-1 ) = zero
621 300 CONTINUE
622 IF( j.GT.1 )
623 $ ap( jc+j-2 ) = cmplx( -one, -one )
624 ap( jc+j-1 ) = tscal*
clarnd( 5, iseed )
625 jc = jc + j
626 310 CONTINUE
627 b( n ) = cmplx( one, one )
628 ELSE
629 jc = 1
630 DO 330 j = 1, n
631 DO 320 i = j + 2, n
632 ap( jc+i-j ) = zero
633 320 CONTINUE
634 IF( j.LT.n )
635 $ ap( jc+1 ) = cmplx( -one, -one )
636 ap( jc ) = tscal*
clarnd( 5, iseed )
637 jc = jc + n - j + 1
638 330 CONTINUE
639 b( 1 ) = cmplx( one, one )
640 END IF
641
642 ELSE IF( imat.EQ.16 ) THEN
643
644
645
646 iy = n / 2 + 1
647 IF( upper ) THEN
648 jc = 1
649 DO 340 j = 1, n
650 CALL clarnv( 4, iseed, j, ap( jc ) )
651 IF( j.NE.iy ) THEN
652 ap( jc+j-1 ) =
clarnd( 5, iseed )*two
653 ELSE
654 ap( jc+j-1 ) = zero
655 END IF
656 jc = jc + j
657 340 CONTINUE
658 ELSE
659 jc = 1
660 DO 350 j = 1, n
661 CALL clarnv( 4, iseed, n-j+1, ap( jc ) )
662 IF( j.NE.iy ) THEN
663 ap( jc ) =
clarnd( 5, iseed )*two
664 ELSE
665 ap( jc ) = zero
666 END IF
667 jc = jc + n - j + 1
668 350 CONTINUE
669 END IF
670 CALL clarnv( 2, iseed, n, b )
671 CALL csscal( n, two, b, 1 )
672
673 ELSE IF( imat.EQ.17 ) THEN
674
675
676
677
678
679
680 tscal = unfl / ulp
681 tscal = ( one-ulp ) / tscal
682 DO 360 j = 1, n*( n+1 ) / 2
683 ap( j ) = zero
684 360 CONTINUE
685 texp = one
686 IF( upper ) THEN
687 jc = ( n-1 )*n / 2 + 1
688 DO 370 j = n, 2, -2
689 ap( jc ) = -tscal / real( n+1 )
690 ap( jc+j-1 ) = one
691 b( j ) = texp*( one-ulp )
692 jc = jc - j + 1
693 ap( jc ) = -( tscal / real( n+1 ) ) / real( n+2 )
694 ap( jc+j-2 ) = one
695 b( j-1 ) = texp*real( n*n+n-1 )
696 texp = texp*two
697 jc = jc - j + 2
698 370 CONTINUE
699 b( 1 ) = ( real( n+1 ) / real( n+2 ) )*tscal
700 ELSE
701 jc = 1
702 DO 380 j = 1, n - 1, 2
703 ap( jc+n-j ) = -tscal / real( n+1 )
704 ap( jc ) = one
705 b( j ) = texp*( one-ulp )
706 jc = jc + n - j + 1
707 ap( jc+n-j-1 ) = -( tscal / real( n+1 ) ) / real( n+2 )
708 ap( jc ) = one
709 b( j+1 ) = texp*real( n*n+n-1 )
710 texp = texp*two
711 jc = jc + n - j
712 380 CONTINUE
713 b( n ) = ( real( n+1 ) / real( n+2 ) )*tscal
714 END IF
715
716 ELSE IF( imat.EQ.18 ) THEN
717
718
719
720
721
722 IF( upper ) THEN
723 jc = 1
724 DO 390 j = 1, n
725 CALL clarnv( 4, iseed, j-1, ap( jc ) )
726 ap( jc+j-1 ) = zero
727 jc = jc + j
728 390 CONTINUE
729 ELSE
730 jc = 1
731 DO 400 j = 1, n
732 IF( j.LT.n )
733 $
CALL clarnv( 4, iseed, n-j, ap( jc+1 ) )
734 ap( jc ) = zero
735 jc = jc + n - j + 1
736 400 CONTINUE
737 END IF
738
739
740
741 CALL clarnv( 2, iseed, n, b )
743 bnorm = abs( b( iy ) )
744 bscal = bignum / max( one, bnorm )
745 CALL csscal( n, bscal, b, 1 )
746
747 ELSE IF( imat.EQ.19 ) THEN
748
749
750
751
752
753
754 tleft = bignum / max( one, real( n-1 ) )
755 tscal = bignum*( real( n-1 ) / max( one, real( n ) ) )
756 IF( upper ) THEN
757 jc = 1
758 DO 420 j = 1, n
759 CALL clarnv( 5, iseed, j, ap( jc ) )
760 CALL slarnv( 1, iseed, j, rwork )
761 DO 410 i = 1, j
762 ap( jc+i-1 ) = ap( jc+i-1 )*( tleft+rwork( i )*tscal )
763 410 CONTINUE
764 jc = jc + j
765 420 CONTINUE
766 ELSE
767 jc = 1
768 DO 440 j = 1, n
769 CALL clarnv( 5, iseed, n-j+1, ap( jc ) )
770 CALL slarnv( 1, iseed, n-j+1, rwork )
771 DO 430 i = j, n
772 ap( jc+i-j ) = ap( jc+i-j )*
773 $ ( tleft+rwork( i-j+1 )*tscal )
774 430 CONTINUE
775 jc = jc + n - j + 1
776 440 CONTINUE
777 END IF
778 CALL clarnv( 2, iseed, n, b )
779 CALL csscal( n, two, b, 1 )
780 END IF
781
782
783
784
785 IF( .NOT.
lsame( trans,
'N' ) )
THEN
786 IF( upper ) THEN
787 jj = 1
788 jr = n*( n+1 ) / 2
789 DO 460 j = 1, n / 2
790 jl = jj
791 DO 450 i = j, n - j
792 t = real( ap( jr-i+j ) )
793 ap( jr-i+j ) = ap( jl )
794 ap( jl ) = t
795 jl = jl + i
796 450 CONTINUE
797 jj = jj + j + 1
798 jr = jr - ( n-j+1 )
799 460 CONTINUE
800 ELSE
801 jl = 1
802 jj = n*( n+1 ) / 2
803 DO 480 j = 1, n / 2
804 jr = jj
805 DO 470 i = j, n - j
806 t = real( ap( jl+i-j ) )
807 ap( jl+i-j ) = ap( jr )
808 ap( jr ) = t
809 jr = jr - i
810 470 CONTINUE
811 jl = jl + n - j + 1
812 jj = jj - j - 1
813 480 CONTINUE
814 END IF
815 END IF
816
817 RETURN
818
819
820
complex function clarnd(idist, iseed)
CLARND
subroutine clatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
CLATB4
subroutine clatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
CLATMS
integer function icamax(n, cx, incx)
ICAMAX
real function slamch(cmach)
SLAMCH
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
logical function lsame(ca, cb)
LSAME
subroutine crot(n, cx, incx, cy, incy, c, s)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
subroutine crotg(a, b, c, s)
CROTG generates a Givens rotation with real cosine and complex sine.
subroutine csscal(n, sa, cx, incx)
CSSCAL