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