184 IMPLICIT NONE
185
186
187
188
189
190
191 INTEGER NM, NN, NNB, NNS, NOUT
192 DOUBLE PRECISION THRESH
193
194
195 LOGICAL DOTYPE( * )
196 INTEGER IWORK( * ), NBVAL( * ), MVAL( * ), NVAL( * ),
197 $ NSVAL( * ), NXVAL( * )
198 DOUBLE PRECISION A( * ), COPYA( * ), B( * ), COPYB( * ),
199 $ S( * ), TAU( * ), WORK( * )
200
201
202
203
204
205 INTEGER NTYPES
206 parameter( ntypes = 19 )
207 INTEGER NTESTS
208 parameter( ntests = 5 )
209 DOUBLE PRECISION ONE, ZERO, BIGNUM
210 parameter( one = 1.0d+0, zero = 0.0d+0,
211 $ bignum = 1.0d+38 )
212
213
214 CHARACTER DIST, TYPE
215 CHARACTER*3 PATH
216 INTEGER I, IHIGH, ILOW, IM, IMAT, IN, INC_ZERO,
217 $ INB, IND_OFFSET_GEN,
218 $ IND_IN, IND_OUT, INS, INFO,
219 $ ISTEP, J, J_INC, J_FIRST_NZ, JB_ZERO,
220 $ KFACT, KL, KMAX, KU, LDA, LW, LWORK,
221 $ LWORK_MQR, M, MINMN, MINMNB_GEN, MODE, N,
222 $ NB, NB_ZERO, NERRS, NFAIL, NB_GEN, NRHS,
223 $ NRUN, NX, T
224 DOUBLE PRECISION ANORM, CNDNUM, EPS, ABSTOL, RELTOL,
225 $ DTEMP, MAXC2NRMK, RELMAXC2NRMK
226
227
228 INTEGER ISEED( 4 ), ISEEDY( 4 )
229 DOUBLE PRECISION RESULT( NTESTS ), RDUMMY( 1 )
230
231
232 DOUBLE PRECISION DLAMCH, DQPT01, DQRT11, DQRT12, DLANGE,
233 $ DLAPY2
235
236
240
241
242 INTRINSIC abs, dble, max, min, mod
243
244
245 LOGICAL LERR, OK
246 CHARACTER*32 SRNAMT
247 INTEGER INFOT, IOUNIT
248
249
250 COMMON / infoc / infot, iounit, ok, lerr
251 COMMON / srnamc / srnamt
252
253
254 DATA iseedy / 1988, 1989, 1990, 1991 /
255
256
257
258
259
260 path( 1: 1 ) = 'Double precision'
261 path( 2: 3 ) = 'QK'
262 nrun = 0
263 nfail = 0
264 nerrs = 0
265 DO i = 1, 4
266 iseed( i ) = iseedy( i )
267 END DO
269 infot = 0
270
271 DO im = 1, nm
272
273
274
275 m = mval( im )
276 lda = max( 1, m )
277
278 DO in = 1, nn
279
280
281
282 n = nval( in )
283 minmn = min( m, n )
284 lwork = max( 1, m*max( m, n )+4*minmn+max( m, n ),
285 $ m*n + 2*minmn + 4*n )
286
287 DO ins = 1, nns
288 nrhs = nsval( ins )
289
290
291
292
293
294
295
296 CALL dlatb4( path, 14, m, nrhs,
TYPE, KL, KU, ANORM,
297 $ MODE, CNDNUM, DIST )
298
299 srnamt = 'DLATMS'
300 CALL dlatms( m, nrhs, dist, iseed,
TYPE, S, MODE,
301 $ CNDNUM, ANORM, KL, KU, 'No packing',
302 $ COPYB, LDA, WORK, INFO )
303
304
305
306
307
308 IF( info.NE.0 ) THEN
309 CALL alaerh( path,
'DLATMS', info, 0,
' ', m,
310 $ nrhs, -1, -1, -1, 6, nfail, nerrs,
311 $ nout )
312 cycle
313 END IF
314
315 DO imat = 1, ntypes
316
317
318
319 IF( .NOT.dotype( imat ) )
320 $ cycle
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349 IF( imat.EQ.1 ) THEN
350
351
352
353 CALL dlaset(
'Full', m, n, zero, zero, copya, lda )
354 DO i = 1, minmn
355 s( i ) = zero
356 END DO
357
358 ELSE IF( (imat.GE.2 .AND. imat.LE.4 )
359 $ .OR. (imat.GE.14 .AND. imat.LE.19 ) ) THEN
360
361
362
363
364
365
366 CALL dlatb4( path, imat, m, n,
TYPE, KL, KU, ANORM,
367 $ MODE, CNDNUM, DIST )
368
369 srnamt = 'DLATMS'
370 CALL dlatms( m, n, dist, iseed,
TYPE, S, MODE,
371 $ CNDNUM, ANORM, KL, KU, 'No packing',
372 $ COPYA, LDA, WORK, INFO )
373
374
375
376 IF( info.NE.0 ) THEN
377 CALL alaerh( path,
'DLATMS', info, 0,
' ', m, n,
378 $ -1, -1, -1, imat, nfail, nerrs,
379 $ nout )
380 cycle
381 END IF
382
383 CALL dlaord(
'Decreasing', minmn, s, 1 )
384
385 ELSE IF( minmn.GE.2
386 $ .AND. imat.GE.5 .AND. imat.LE.13 ) THEN
387
388
389
390
391
392
393
394
395
396
397
398
399
400 IF( imat.EQ.5 ) THEN
401
402
403
404 jb_zero = 1
405 nb_zero = 1
406 nb_gen = n - nb_zero
407
408 ELSE IF( imat.EQ.6 ) THEN
409
410
411
412 jb_zero = minmn
413 nb_zero = 1
414 nb_gen = n - nb_zero
415
416 ELSE IF( imat.EQ.7 ) THEN
417
418
419
420 jb_zero = n
421 nb_zero = 1
422 nb_gen = n - nb_zero
423
424 ELSE IF( imat.EQ.8 ) THEN
425
426
427
428 jb_zero = minmn / 2 + 1
429 nb_zero = 1
430 nb_gen = n - nb_zero
431
432 ELSE IF( imat.EQ.9 ) THEN
433
434
435
436 jb_zero = 1
437 nb_zero = minmn / 2
438 nb_gen = n - nb_zero
439
440 ELSE IF( imat.EQ.10 ) THEN
441
442
443
444
445 jb_zero = minmn / 2 + 1
446 nb_zero = n - jb_zero + 1
447 nb_gen = n - nb_zero
448
449 ELSE IF( imat.EQ.11 ) THEN
450
451
452
453
454
455 jb_zero = minmn / 2 - (minmn / 2) / 2 + 1
456 nb_zero = minmn / 2
457 nb_gen = n - nb_zero
458
459 ELSE IF( imat.EQ.12 ) THEN
460
461
462
463 nb_gen = n / 2
464 nb_zero = n - nb_gen
465 j_inc = 2
466 j_first_nz = 2
467
468 ELSE IF( imat.EQ.13 ) THEN
469
470
471
472 nb_zero = n / 2
473 nb_gen = n - nb_zero
474 j_inc = 2
475 j_first_nz = 1
476
477 END IF
478
479
480
481
482
483 CALL dlaset(
'Full', m, nb_zero, zero, zero,
484 $ copya, lda )
485
486
487
488
489
490 CALL dlatb4( path, imat, m, nb_gen,
TYPE, KL, KU,
491 $ ANORM, MODE, CNDNUM, DIST )
492
493 srnamt = 'DLATMS'
494
495 ind_offset_gen = nb_zero * lda
496
497 CALL dlatms( m, nb_gen, dist, iseed,
TYPE, S, MODE,
498 $ CNDNUM, ANORM, KL, KU, 'No packing',
499 $ COPYA( IND_OFFSET_GEN + 1 ), LDA,
500 $ WORK, INFO )
501
502
503
504 IF( info.NE.0 ) THEN
505 CALL alaerh( path,
'DLATMS', info, 0,
' ', m,
506 $ nb_gen, -1, -1, -1, imat, nfail,
507 $ nerrs, nout )
508 cycle
509 END IF
510
511
512
513
514
515 IF( imat.EQ.6
516 $ .OR. imat.EQ.7
517 $ .OR. imat.EQ.8
518 $ .OR. imat.EQ.10
519 $ .OR. imat.EQ.11 ) THEN
520
521
522
523
524
525
526 DO j = 1, jb_zero-1, 1
528 $ copya( ( nb_zero+j-1)*lda+1), 1,
529 $ copya( (j-1)*lda + 1 ), 1 )
530 END DO
531
532 ELSE IF( imat.EQ.12 .OR. imat.EQ.13 ) THEN
533
534
535
536
537
538
539
540
541
542
543
544 DO j = 1, nb_gen, 1
545 ind_out = ( nb_zero+j-1 )*lda + 1
546 ind_in = ( j_inc*(j-1)+(j_first_nz-1) )*lda
547 $ + 1
549 $ copya( ind_out ), 1,
550 $ copya( ind_in), 1 )
551 END DO
552
553 END IF
554
555
556
557
558
559
560 minmnb_gen = min( m, nb_gen )
561
562 DO i = minmnb_gen+1, minmn
563 s( i ) = zero
564 END DO
565
566 ELSE
567
568
569
570 cycle
571 END IF
572
573
574
575 DO i = 1, n
576 iwork( i ) = 0
577 END DO
578
579 DO inb = 1, nnb
580
581
582
583 nb = nbval( inb )
585 nx = nxval( inb )
587
588
589
590
591
592 DO kmax = 0, min(m,n)+1
593
594
595
596
597
598
599
600
601
602 CALL dlacpy(
'All', m, n, copya, lda, a, lda )
603 CALL dlacpy(
'All', m, nrhs, copyb, lda,
604 $ a( lda*n + 1 ), lda )
605 CALL dlacpy(
'All', m, nrhs, copyb, lda,
606 $ b, lda )
607 CALL icopy( n, iwork( 1 ), 1, iwork( n+1 ), 1 )
608 DO i = 1, ntests
609 result( i ) = zero
610 END DO
611
612 abstol = -1.0
613 reltol = -1.0
614
615
616
617 lw = max( 1, max( 2*n + nb*( n+nrhs+1 ),
618 $ 3*n + nrhs - 1 ) )
619
620
621
622 srnamt = 'DGEQP3RK'
623 CALL dgeqp3rk( m, n, nrhs, kmax, abstol, reltol,
624 $ a, lda, kfact, maxc2nrmk,
625 $ relmaxc2nrmk, iwork( n+1 ), tau,
626 $ work, lw, iwork( 2*n+1 ), info )
627
628
629
630 IF( info.LT.0 )
631 $
CALL alaerh( path,
'DGEQP3RK', info, 0,
' ',
632 $ m, n, nx, -1, nb, imat,
633 $ nfail, nerrs, nout )
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649 IF( kfact.EQ.minmn ) THEN
650
651 result( 1 ) =
dqrt12( m, n, a, lda, s, work,
652 $ lwork )
653
654 nrun = nrun + 1
655
656
657
658 END IF
659
660
661
662
663
664
665
666 result( 2 ) =
dqpt01( m, n, kfact, copya, a, lda, tau,
667 $ iwork( n+1 ), work, lwork )
668
669
670
671
672
673
674
675 result( 3 ) =
dqrt11( m, kfact, a, lda, tau, work,
676 $ lwork )
677
678 nrun = nrun + 2
679
680
681
682
683
684
685
686
687
688
689
690
691
692 IF( min(kfact, minmn).GE.2 ) THEN
693
694 DO j = 1, kfact-1, 1
695
696 dtemp = (( abs( a( (j-1)*lda+j ) ) -
697 $ abs( a( (j)*lda+j+1 ) ) ) /
698 $ abs( a(1) ) )
699
700 IF( dtemp.LT.zero ) THEN
701 result( 4 ) = bignum
702 END IF
703
704 END DO
705
706 nrun = nrun + 1
707
708
709
710 END IF
711
712
713
714
715
716
717
718
719
720
721
722
723 IF( minmn.GT.0 ) THEN
724
725 lwork_mqr = max(1, nrhs)
726 CALL dormqr(
'Left',
'Transpose',
727 $ m, nrhs, kfact, a, lda, tau, b, lda,
728 $ work, lwork_mqr, info )
729
730 DO i = 1, nrhs
731
732
733
734 CALL daxpy( m, -one, a( ( n+i-1 )*lda+1 ), 1,
735 $ b( ( i-1 )*lda+1 ), 1 )
736 END DO
737
738 result( 5 ) = abs(
739 $
dlange(
'One-norm', m, nrhs, b, lda, rdummy ) /
740 $ ( dble( m )*
dlamch(
'Epsilon' ) ) )
741
742 nrun = nrun + 1
743
744
745
746 END IF
747
748
749
750
751 DO t = 1, ntests
752 IF( result( t ).GE.thresh ) THEN
753 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
754 $
CALL alahd( nout, path )
755 WRITE( nout, fmt = 9999 ) 'DGEQP3RK', m, n,
756 $ nrhs, kmax, abstol, reltol, nb, nx,
757 $ imat, t, result( t )
758 nfail = nfail + 1
759 END IF
760 END DO
761
762
763
764 END DO
765
766
767
768 END DO
769
770
771
772 END DO
773
774
775
776 END DO
777
778
779
780 END DO
781
782
783
784 END DO
785
786
787
788 CALL alasum( path, nout, nfail, nrun, nerrs )
789
790 9999 FORMAT( 1x, a, ' M =', i5, ', N =', i5, ', NRHS =', i5,
791 $ ', KMAX =', i5, ', ABSTOL =', g12.5,
792 $ ', RELTOL =', g12.5, ', NB =', i4, ', NX =', i4,
793 $ ', type ', i2, ', test ', i2, ', ratio =', g12.5 )
794
795
796
subroutine alasum(type, nout, nfail, nrun, nerrs)
ALASUM
subroutine xlaenv(ispec, nvalue)
XLAENV
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
subroutine alahd(iounit, path)
ALAHD
subroutine dgeqp3rk(m, n, nrhs, kmax, abstol, reltol, a, lda, k, maxc2nrmk, relmaxc2nrmk, jpiv, tau, work, lwork, iwork, info)
DGEQP3RK computes a truncated Householder QR factorization with column pivoting of a real m-by-n matr...
subroutine dlaord(job, n, x, incx)
DLAORD
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
double precision function dqpt01(m, n, k, a, af, lda, tau, jpvt, work, lwork)
DQPT01
double precision function dqrt11(m, k, a, lda, tau, work, lwork)
DQRT11
double precision function dqrt12(m, n, a, lda, s, work, lwork)
DQRT12
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
subroutine icopy(n, sx, incx, sy, incy)
ICOPY