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