From a50de3c8cc87505adbfdb759def537a44da129a3 Mon Sep 17 00:00:00 2001 From: Julie Date: Sun, 9 Oct 2016 10:07:56 -0700 Subject: Adding Assen for SV, TRF, SV from @iyamazak@icl.utk.edu Contribution by Ichitaro Yamazaki, University of Tennessee Note: need to add corresponding LAPACKE routines --- TESTING/LIN/CMakeLists.txt | 24 +- TESTING/LIN/Makefile | 24 +- TESTING/LIN/aladhd.f | 26 +- TESTING/LIN/alaerh.f | 3 + TESTING/LIN/alahd.f | 24 ++ TESTING/LIN/cchkaa.f | 36 ++- TESTING/LIN/cchkhe_aasen.f | 578 ++++++++++++++++++++++++++++++++++++++++++++ TESTING/LIN/cdrvhe_aasen.f | 529 +++++++++++++++++++++++++++++++++++++++++ TESTING/LIN/cerrvx.f | 21 +- TESTING/LIN/chet01_aasen.f | 267 +++++++++++++++++++++ TESTING/LIN/dchkaa.f | 29 +++ TESTING/LIN/dchksy_aasen.f | 577 ++++++++++++++++++++++++++++++++++++++++++++ TESTING/LIN/ddrvsy_aasen.f | 517 ++++++++++++++++++++++++++++++++++++++++ TESTING/LIN/derrvx.f | 21 +- TESTING/LIN/dsyt01_aasen.f | 263 ++++++++++++++++++++ TESTING/LIN/schkaa.f | 27 +++ TESTING/LIN/schksy_aasen.f | 578 ++++++++++++++++++++++++++++++++++++++++++++ TESTING/LIN/sdrvsy_aasen.f | 517 ++++++++++++++++++++++++++++++++++++++++ TESTING/LIN/serrvx.f | 20 +- TESTING/LIN/ssyt01_aasen.f | 262 ++++++++++++++++++++ TESTING/LIN/zchkaa.f | 34 ++- TESTING/LIN/zchkhe_aasen.f | 579 +++++++++++++++++++++++++++++++++++++++++++++ TESTING/LIN/zdrvhe_aasen.f | 525 ++++++++++++++++++++++++++++++++++++++++ TESTING/LIN/zerrvx.f | 21 +- TESTING/LIN/zhet01_aasen.f | 267 +++++++++++++++++++++ 25 files changed, 5732 insertions(+), 37 deletions(-) create mode 100644 TESTING/LIN/cchkhe_aasen.f create mode 100644 TESTING/LIN/cdrvhe_aasen.f create mode 100644 TESTING/LIN/chet01_aasen.f create mode 100644 TESTING/LIN/dchksy_aasen.f create mode 100644 TESTING/LIN/ddrvsy_aasen.f create mode 100644 TESTING/LIN/dsyt01_aasen.f create mode 100644 TESTING/LIN/schksy_aasen.f create mode 100644 TESTING/LIN/sdrvsy_aasen.f create mode 100644 TESTING/LIN/ssyt01_aasen.f create mode 100644 TESTING/LIN/zchkhe_aasen.f create mode 100644 TESTING/LIN/zdrvhe_aasen.f create mode 100644 TESTING/LIN/zhet01_aasen.f (limited to 'TESTING/LIN') diff --git a/TESTING/LIN/CMakeLists.txt b/TESTING/LIN/CMakeLists.txt index a924b12f..b8320977 100644 --- a/TESTING/LIN/CMakeLists.txt +++ b/TESTING/LIN/CMakeLists.txt @@ -10,10 +10,10 @@ set(SLINTST schkaa.f schkeq.f schkgb.f schkge.f schkgt.f schklq.f schkpb.f schkpo.f schkps.f schkpp.f schkpt.f schkq3.f schkql.f schkqr.f schkrq.f - schksp.f schksy.f schksy_rook.f schktb.f schktp.f schktr.f + schksp.f schksy.f schksy_rook.f schksy_aasen.f schktb.f schktp.f schktr.f schktz.f sdrvgt.f sdrvls.f sdrvpb.f - sdrvpp.f sdrvpt.f sdrvsp.f sdrvsy.f sdrvsy_rook.f + sdrvpp.f sdrvpt.f sdrvsp.f sdrvsy.f sdrvsy_rook.f sdrvsy_aasen.f serrgt.f serrlq.f serrls.f serrpo.f serrps.f serrql.f serrqp.f serrqr.f serrrq.f serrsy.f serrtr.f serrtz.f serrvx.f @@ -29,7 +29,7 @@ set(SLINTST schkaa.f sqrt01.f sqrt01p.f sqrt02.f sqrt03.f sqrt11.f sqrt12.f sqrt13.f sqrt14.f sqrt15.f sqrt16.f sqrt17.f srqt01.f srqt02.f srqt03.f srzt01.f srzt02.f - sspt01.f ssyt01.f ssyt01_rook.f + sspt01.f ssyt01.f ssyt01_rook.f ssyt01_aasen.f stbt02.f stbt03.f stbt05.f stbt06.f stpt01.f stpt02.f stpt03.f stpt05.f stpt06.f strt01.f strt02.f strt03.f strt05.f strt06.f @@ -44,11 +44,11 @@ endif() set(CLINTST cchkaa.f cchkeq.f cchkgb.f cchkge.f cchkgt.f - cchkhe.f cchkhe_rook.f cchkhp.f cchklq.f cchkpb.f + cchkhe.f cchkhe_rook.f cchkhe_aasen.f cchkhp.f cchklq.f cchkpb.f cchkpo.f cchkps.f cchkpp.f cchkpt.f cchkq3.f cchkql.f cchkqr.f cchkrq.f cchksp.f cchksy.f cchksy_rook.f cchktb.f cchktp.f cchktr.f cchktz.f - cdrvgt.f cdrvhe.f cdrvhe_rook.f cdrvhp.f + cdrvgt.f cdrvhe.f cdrvhe_rook.f cdrvhe_aasen.f cdrvhp.f cdrvls.f cdrvpb.f cdrvpp.f cdrvpt.f cdrvsp.f cdrvsy.f cdrvsy_rook.f cerrgt.f cerrhe.f cerrlq.f @@ -58,7 +58,7 @@ set(CLINTST cchkaa.f cgbt01.f cgbt02.f cgbt05.f cgelqs.f cgeqls.f cgeqrs.f cgerqs.f cget01.f cget02.f cget03.f cget04.f cget07.f cgtt01.f cgtt02.f - cgtt05.f chet01.f chet01_rook.f chpt01.f claipd.f claptm.f clarhs.f clatb4.f clatb5.f + cgtt05.f chet01.f chet01_rook.f chet01_aasen.f chpt01.f claipd.f claptm.f clarhs.f clatb4.f clatb5.f clatsp.f clatsy.f clattb.f clattp.f clattr.f clavhe.f clavhe_rook.f clavhp.f clavsp.f clavsy.f clavsy_rook.f clqt01.f clqt02.f clqt03.f cpbt01.f cpbt02.f cpbt05.f @@ -87,10 +87,10 @@ set(DLINTST dchkaa.f dchkeq.f dchkgb.f dchkge.f dchkgt.f dchklq.f dchkpb.f dchkpo.f dchkps.f dchkpp.f dchkpt.f dchkq3.f dchkql.f dchkqr.f dchkrq.f - dchksp.f dchksy.f dchksy_rook.f dchktb.f dchktp.f dchktr.f + dchksp.f dchksy.f dchksy_rook.f dchksy_aasen.f dchktb.f dchktp.f dchktr.f dchktz.f ddrvgt.f ddrvls.f ddrvpb.f - ddrvpp.f ddrvpt.f ddrvsp.f ddrvsy.f ddrvsy_rook.f + ddrvpp.f ddrvpt.f ddrvsp.f ddrvsy.f ddrvsy_rook.f ddrvsy_aasen.f derrgt.f derrlq.f derrls.f derrps.f derrql.f derrqp.f derrqr.f derrrq.f derrsy.f derrtr.f derrtz.f derrvx.f @@ -106,7 +106,7 @@ set(DLINTST dchkaa.f dqrt01.f dqrt01p.f dqrt02.f dqrt03.f dqrt11.f dqrt12.f dqrt13.f dqrt14.f dqrt15.f dqrt16.f dqrt17.f drqt01.f drqt02.f drqt03.f drzt01.f drzt02.f - dspt01.f dsyt01.f dsyt01_rook.f + dspt01.f dsyt01.f dsyt01_rook.f dsyt01_aasen.f dtbt02.f dtbt03.f dtbt05.f dtbt06.f dtpt01.f dtpt02.f dtpt03.f dtpt05.f dtpt06.f dtrt01.f dtrt02.f dtrt03.f dtrt05.f dtrt06.f @@ -123,11 +123,11 @@ endif() set(ZLINTST zchkaa.f zchkeq.f zchkgb.f zchkge.f zchkgt.f - zchkhe.f zchkhe_rook.f zchkhp.f zchklq.f zchkpb.f + zchkhe.f zchkhe_rook.f zchkhe_aasen.f zchkhp.f zchklq.f zchkpb.f zchkpo.f zchkps.f zchkpp.f zchkpt.f zchkq3.f zchkql.f zchkqr.f zchkrq.f zchksp.f zchksy.f zchksy_rook.f zchktb.f zchktp.f zchktr.f zchktz.f - zdrvgt.f zdrvhe.f zdrvhe_rook.f zdrvhp.f + zdrvgt.f zdrvhe.f zdrvhe_rook.f zdrvhe_aasen.f zdrvhp.f zdrvls.f zdrvpb.f zdrvpp.f zdrvpt.f zdrvsp.f zdrvsy.f zdrvsy_rook.f zerrgt.f zerrhe.f zerrlq.f @@ -137,7 +137,7 @@ set(ZLINTST zchkaa.f zgbt01.f zgbt02.f zgbt05.f zgelqs.f zgeqls.f zgeqrs.f zgerqs.f zget01.f zget02.f zget03.f zget04.f zget07.f zgtt01.f zgtt02.f - zgtt05.f zhet01.f zhet01.f zhet01_rook.f zhpt01.f zlaipd.f zlaptm.f zlarhs.f zlatb4.f zlatb5.f + zgtt05.f zhet01.f zhet01.f zhet01_rook.f zhet01_aasen.f zhpt01.f zlaipd.f zlaptm.f zlarhs.f zlatb4.f zlatb5.f zlatsp.f zlatsy.f zlattb.f zlattp.f zlattr.f zlavhe.f zlavhe_rook.f zlavhp.f zlavsp.f zlavsy.f zlavsy_rook.f zlqt01.f zlqt02.f zlqt03.f zpbt01.f zpbt02.f zpbt05.f diff --git a/TESTING/LIN/Makefile b/TESTING/LIN/Makefile index 57d4fe25..3232f0fd 100644 --- a/TESTING/LIN/Makefile +++ b/TESTING/LIN/Makefile @@ -51,10 +51,10 @@ SLINTST = schkaa.o \ schkeq.o schkgb.o schkge.o schkgt.o \ schklq.o schkpb.o schkpo.o schkps.o schkpp.o \ schkpt.o schkq3.o schkql.o schkqr.o schkrq.o \ - schksp.o schksy.o schksy_rook.o schktb.o schktp.o schktr.o \ + schksp.o schksy.o schksy_rook.o schksy_aasen.o schktb.o schktp.o schktr.o \ schktz.o \ sdrvgt.o sdrvls.o sdrvpb.o \ - sdrvpp.o sdrvpt.o sdrvsp.o sdrvsy_rook.o\ + sdrvpp.o sdrvpt.o sdrvsp.o sdrvsy_rook.o sdrvsy_aasen.o\ serrgt.o serrlq.o serrls.o \ serrps.o serrql.o serrqp.o serrqr.o \ serrrq.o serrtr.o serrtz.o \ @@ -70,7 +70,7 @@ SLINTST = schkaa.o \ sqrt01.o sqrt01p.o sqrt02.o sqrt03.o sqrt11.o sqrt12.o \ sqrt13.o sqrt14.o sqrt15.o sqrt16.o sqrt17.o \ srqt01.o srqt02.o srqt03.o srzt01.o srzt02.o \ - sspt01.o ssyt01.o ssyt01_rook.o \ + sspt01.o ssyt01.o ssyt01_rook.o ssyt01_aasen.o\ stbt02.o stbt03.o stbt05.o stbt06.o stpt01.o \ stpt02.o stpt03.o stpt05.o stpt06.o strt01.o \ strt02.o strt03.o strt05.o strt06.o \ @@ -86,11 +86,11 @@ endif CLINTST = cchkaa.o \ cchkeq.o cchkgb.o cchkge.o cchkgt.o \ - cchkhe.o cchkhe_rook.o cchkhp.o cchklq.o cchkpb.o \ + cchkhe.o cchkhe_rook.o cchkhe_aasen.o cchkhp.o cchklq.o cchkpb.o \ cchkpo.o cchkps.o cchkpp.o cchkpt.o cchkq3.o cchkql.o \ cchkqr.o cchkrq.o cchksp.o cchksy.o cchksy_rook.o cchktb.o \ cchktp.o cchktr.o cchktz.o \ - cdrvgt.o cdrvhe_rook.o cdrvhp.o \ + cdrvgt.o cdrvhe_rook.o cdrvhe_aasen.o cdrvhp.o \ cdrvls.o cdrvpb.o cdrvpp.o cdrvpt.o \ cdrvsp.o cdrvsy_rook.o \ cerrgt.o cerrlq.o \ @@ -99,7 +99,7 @@ CLINTST = cchkaa.o \ cgbt01.o cgbt02.o cgbt05.o cgelqs.o cgeqls.o cgeqrs.o \ cgerqs.o cget01.o cget02.o \ cget03.o cget04.o cget07.o cgtt01.o cgtt02.o \ - cgtt05.o chet01.o chet01_rook.o chpt01.o claipd.o claptm.o clarhs.o clatb4.o clatb5.o \ + cgtt05.o chet01.o chet01_rook.o chet01_aasen.o chpt01.o claipd.o claptm.o clarhs.o clatb4.o clatb5.o \ clatsp.o clatsy.o clattb.o clattp.o clattr.o \ clavhe.o clavhe_rook.o clavhp.o clavsp.o clavsy.o clavsy_rook.o clqt01.o \ clqt02.o clqt03.o cpbt01.o cpbt02.o cpbt05.o \ @@ -129,10 +129,10 @@ DLINTST = dchkaa.o \ dchkeq.o dchkgb.o dchkge.o dchkgt.o \ dchklq.o dchkpb.o dchkpo.o dchkps.o dchkpp.o \ dchkpt.o dchkq3.o dchkql.o dchkqr.o dchkrq.o \ - dchksp.o dchksy.o dchksy_rook.o dchktb.o dchktp.o dchktr.o \ + dchksp.o dchksy.o dchksy_rook.o dchksy_aasen.o dchktb.o dchktp.o dchktr.o \ dchktz.o \ ddrvgt.o ddrvls.o ddrvpb.o \ - ddrvpp.o ddrvpt.o ddrvsp.o ddrvsy_rook.o \ + ddrvpp.o ddrvpt.o ddrvsp.o ddrvsy_rook.o ddrvsy_aasen.o\ derrgt.o derrlq.o derrls.o \ derrps.o derrql.o derrqp.o derrqr.o \ derrrq.o derrtr.o derrtz.o \ @@ -148,7 +148,7 @@ DLINTST = dchkaa.o \ dqrt01.o dqrt01p.o dqrt02.o dqrt03.o dqrt11.o dqrt12.o \ dqrt13.o dqrt14.o dqrt15.o dqrt16.o dqrt17.o \ drqt01.o drqt02.o drqt03.o drzt01.o drzt02.o \ - dspt01.o dsyt01.o dsyt01_rook.o \ + dspt01.o dsyt01.o dsyt01_rook.o dsyt01_aasen.o\ dtbt02.o dtbt03.o dtbt05.o dtbt06.o dtpt01.o \ dtpt02.o dtpt03.o dtpt05.o dtpt06.o dtrt01.o \ dtrt02.o dtrt03.o dtrt05.o dtrt06.o \ @@ -165,11 +165,11 @@ endif ZLINTST = zchkaa.o \ zchkeq.o zchkgb.o zchkge.o zchkgt.o \ - zchkhe.o zchkhe_rook.o zchkhp.o zchklq.o zchkpb.o \ + zchkhe.o zchkhe_rook.o zchkhe_aasen.o zchkhp.o zchklq.o zchkpb.o \ zchkpo.o zchkps.o zchkpp.o zchkpt.o zchkq3.o zchkql.o \ zchkqr.o zchkrq.o zchksp.o zchksy.o zchksy_rook.o zchktb.o \ zchktp.o zchktr.o zchktz.o \ - zdrvgt.o zdrvhe_rook.o zdrvhp.o \ + zdrvgt.o zdrvhe_rook.o zdrvhe_aasen.o zdrvhp.o \ zdrvls.o zdrvpb.o zdrvpp.o zdrvpt.o \ zdrvsp.o zdrvsy_rook.o \ zerrgt.o zerrlq.o \ @@ -178,7 +178,7 @@ ZLINTST = zchkaa.o \ zgbt01.o zgbt02.o zgbt05.o zgelqs.o zgeqls.o zgeqrs.o \ zgerqs.o zget01.o zget02.o \ zget03.o zget04.o zget07.o zgtt01.o zgtt02.o \ - zgtt05.o zhet01.o zhet01_rook.o zhpt01.o zlaipd.o zlaptm.o zlarhs.o zlatb4.o zlatb5.o \ + zgtt05.o zhet01.o zhet01_rook.o zhet01_aasen.o zhpt01.o zlaipd.o zlaptm.o zlarhs.o zlatb4.o zlatb5.o \ zlatsp.o zlatsy.o zlattb.o zlattp.o zlattr.o \ zlavhe.o zlavhe_rook.o zlavhp.o zlavsp.o zlavsy.o zlavsy_rook.o zlqt01.o \ zlqt02.o zlqt03.o zpbt01.o zpbt02.o zpbt05.o \ diff --git a/TESTING/LIN/aladhd.f b/TESTING/LIN/aladhd.f index b2bb2e03..3a53e0bd 100644 --- a/TESTING/LIN/aladhd.f +++ b/TESTING/LIN/aladhd.f @@ -53,6 +53,8 @@ *> with "rook" (bounded Bunch-Kaufman) pivoting *> _SP: Symmetric indefinite packed, *> with partial (Bunch-Kaufman) pivoting +*> _HA: (complex) Hermitian , +*> Assen Algorithm *> _HE: (complex) Hermitian indefinite, *> with partial (Bunch-Kaufman) pivoting *> _HR: (complex) Hermitian indefinite, @@ -275,7 +277,27 @@ WRITE( IOUNIT, FMT = 9979 )3 WRITE( IOUNIT, FMT = '( '' Messages:'' )' ) * - ELSE IF( LSAMEN( 2, P2, 'HE' ) .OR. LSAMEN( 2, P2, 'HP' ) ) THEN + ELSE IF( LSAMEN( 2, P2, 'HA' ) ) THEN +* +* HA: Hermitian +* Aasen algorithm + WRITE( IOUNIT, FMT = 9971 )PATH, 'Hermitian' +* + WRITE( IOUNIT, FMT = '( '' Matrix types:'' )' ) + WRITE( IOUNIT, FMT = 9983 ) +* + WRITE( IOUNIT, FMT = '( '' Test ratios:'' )' ) + WRITE( IOUNIT, FMT = 9974 )1 + WRITE( IOUNIT, FMT = 9980 )2 + WRITE( IOUNIT, FMT = 9979 )3 + WRITE( IOUNIT, FMT = 9977 )4 + WRITE( IOUNIT, FMT = 9978 )5 + WRITE( IOUNIT, FMT = 9976 )6 + WRITE( IOUNIT, FMT = '( '' Messages:'' )' ) + + + ELSE IF( LSAMEN( 2, P2, 'HE' ) .OR. + $ LSAMEN( 2, P2, 'HP' ) ) THEN * * HE: Hermitian indefinite full * with partial (Bunch-Kaufman) pivoting algorithm @@ -336,6 +358,8 @@ $ ' positive definite band matrices' ) 9993 FORMAT( / 1X, A3, ' drivers: ', A9, $ ' positive definite tridiagonal' ) + 9971 FORMAT( / 1X, A3, ' drivers: ', A9, ' indefinite matrices', + $ ', "Aasen" Algorithm' ) 9992 FORMAT( / 1X, A3, ' drivers: ', A9, ' indefinite matrices', $ ', "rook" (bounded Bunch-Kaufman) pivoting' ) 9991 FORMAT( / 1X, A3, ' drivers: ', A9, diff --git a/TESTING/LIN/alaerh.f b/TESTING/LIN/alaerh.f index 2f45b1d2..4fec4522 100644 --- a/TESTING/LIN/alaerh.f +++ b/TESTING/LIN/alaerh.f @@ -490,6 +490,7 @@ ELSE IF( LSAMEN( 2, P2, 'SY' ) $ .OR. LSAMEN( 2, P2, 'SR' ) $ .OR. LSAMEN( 2, P2, 'HE' ) + $ .OR. LSAMEN( 2, P2, 'HA' ) $ .OR. LSAMEN( 2, P2, 'HR' ) ) THEN * * xSY: symmetric indefinite matrices @@ -498,6 +499,8 @@ * with rook (bounded Bunch-Kaufman) pivoting; * xHE: Hermitian indefinite matrices * with partial (Bunch-Kaufman) pivoting. +* xHA: Hermitian matrices +* Aasen Algorithm * xHR: Hermitian indefinite matrices * with rook (bounded Bunch-Kaufman) pivoting; * diff --git a/TESTING/LIN/alahd.f b/TESTING/LIN/alahd.f index 8d56c798..995a5444 100644 --- a/TESTING/LIN/alahd.f +++ b/TESTING/LIN/alahd.f @@ -53,6 +53,8 @@ *> with "rook" (bounded Bunch-Kaufman) pivoting *> _SP: Symmetric indefinite packed, *> with partial (Bunch-Kaufman) pivoting +*> _HA: (complex) Hermitian , +*> with Aasen Algorithm *> _HE: (complex) Hermitian indefinite, *> with partial (Bunch-Kaufman) pivoting *> _HR: Symmetric indefinite, @@ -354,6 +356,28 @@ WRITE( IOUNIT, FMT = 9957 )7 WRITE( IOUNIT, FMT = 9955 )8 WRITE( IOUNIT, FMT = '( '' Messages:'' )' ) +* + ELSE IF( LSAMEN( 2, P2, 'HA' ) ) THEN +* +* HA: Hermitian, +* with Assen Algorithm +* + WRITE( IOUNIT, FMT = 9992 )PATH, 'Hermitian' +* + WRITE( IOUNIT, FMT = '( '' Matrix types:'' )' ) + WRITE( IOUNIT, FMT = 9972 ) +* + WRITE( IOUNIT, FMT = '( '' Test ratios:'' )' ) + WRITE( IOUNIT, FMT = 9953 )1 + WRITE( IOUNIT, FMT = 9961 )2 + WRITE( IOUNIT, FMT = 9960 )3 + WRITE( IOUNIT, FMT = 9960 )4 + WRITE( IOUNIT, FMT = 9959 )5 + WRITE( IOUNIT, FMT = 9958 )6 + WRITE( IOUNIT, FMT = 9956 )7 + WRITE( IOUNIT, FMT = 9957 )8 + WRITE( IOUNIT, FMT = 9955 )9 + WRITE( IOUNIT, FMT = '( '' Messages:'' )' ) * ELSE IF( LSAMEN( 2, P2, 'HE' ) ) THEN * diff --git a/TESTING/LIN/cchkaa.f b/TESTING/LIN/cchkaa.f index 8a97fbb5..1f8b2c65 100644 --- a/TESTING/LIN/cchkaa.f +++ b/TESTING/LIN/cchkaa.f @@ -51,6 +51,7 @@ *> CPT 12 List types on next line if 0 < NTYPES < 12 *> CHE 10 List types on next line if 0 < NTYPES < 10 *> CHR 10 List types on next line if 0 < NTYPES < 10 +*> CHA 10 List types on next line if 0 < NTYPES < 10 *> CHP 10 List types on next line if 0 < NTYPES < 10 *> CSY 11 List types on next line if 0 < NTYPES < 11 *> CSR 11 List types on next line if 0 < NTYPES < 11 @@ -102,17 +103,17 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date November 2015 +*> \date November 2016 * *> \ingroup complex_lin * * ===================================================================== PROGRAM CCHKAA * -* -- LAPACK test routine (version 3.6.0) -- +* -- LAPACK test routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2015 +* November 2016 * * ===================================================================== * @@ -640,6 +641,33 @@ ELSE WRITE( NOUT, FMT = 9988 )PATH END IF +* + ELSE IF( LSAMEN( 2, C2, 'HA' ) ) THEN +* +* HA: Hermitian matrices, +* Aasen Algorithm +* + NTYPES = 10 + CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) +* + IF( TSTCHK ) THEN + CALL CCHKHE_AASEN( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS, + $ NSVAL, THRESH, TSTERR, LDA, + $ A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), + $ WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9989 )PATH + END IF +* + IF( TSTDRV ) THEN + CALL CDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), + $ WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9988 )PATH + END IF * ELSE IF( LSAMEN( 2, C2, 'HR' ) ) THEN * @@ -867,7 +895,6 @@ ELSE WRITE( NOUT, FMT = 9989 )PATH END IF - * ELSE IF( LSAMEN( 2, C2, 'RQ' ) ) THEN * @@ -926,7 +953,6 @@ ELSE WRITE( NOUT, FMT = 9989 )PATH END IF - * ELSE IF( LSAMEN( 2, C2, 'LS' ) ) THEN * diff --git a/TESTING/LIN/cchkhe_aasen.f b/TESTING/LIN/cchkhe_aasen.f new file mode 100644 index 00000000..355dd25b --- /dev/null +++ b/TESTING/LIN/cchkhe_aasen.f @@ -0,0 +1,578 @@ +*> \brief \b CCHKHE_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE CCHKHE_AASEN( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, +* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, +* XACT, WORK, RWORK, IWORK, NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NN, NNB, NNS, NOUT +* REAL THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) +* REAL RWORK( * ) +* COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), +* $ WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CCHKHE_AASEN tests CHETRF_AASEN, -TRS_AASEN. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NNB +*> \verbatim +*> NNB is INTEGER +*> The number of values of NB contained in the vector NBVAL. +*> \endverbatim +*> +*> \param[in] NBVAL +*> \verbatim +*> NBVAL is INTEGER array, dimension (NBVAL) +*> The values of the blocksize NB. +*> \endverbatim +*> +*> \param[in] NNS +*> \verbatim +*> NNS is INTEGER +*> The number of values of NRHS contained in the vector NSVAL. +*> \endverbatim +*> +*> \param[in] NSVAL +*> \verbatim +*> NSVAL is INTEGER array, dimension (NNS) +*> The values of the number of right hand sides NRHS. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is REAL +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is COMPLEX array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is COMPLEX array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is COMPLEX array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is COMPLEX array, dimension (NMAX*NSMAX) +*> where NSMAX is the largest entry in NSVAL. +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is COMPLEX array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is COMPLEX array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (NMAX*max(3,NSMAX)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (max(NMAX,2*NSMAX)) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* +*> \ingroup complex_lin +* +* ===================================================================== + SUBROUTINE CCHKHE_AASEN( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, + $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, + $ X, XACT, WORK, RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NMAX, NN, NNB, NNS, NOUT + REAL THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) + REAL RWORK( * ) + COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), + $ WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO + PARAMETER ( ZERO = 0.0E+0 ) + COMPLEX CZERO + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) + INTEGER NTYPES + PARAMETER ( NTYPES = 10 ) + INTEGER NTESTS + PARAMETER ( NTESTS = 9 ) +* .. +* .. Local Scalars .. + LOGICAL TRFCON, ZEROT + CHARACTER DIST, TYPE, UPLO, XTYPE + CHARACTER*3 PATH, MATPATH + INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, + $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, + $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT + REAL ANORM, CNDNUM, RCOND, RCONDC +* .. +* .. Local Arrays .. + CHARACTER UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + REAL RESULT( NTESTS ) +* .. +* .. External Functions .. + REAL DGET06, CLANHE + EXTERNAL DGET06, CLANHE +* .. +* .. External Subroutines .. + EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, CERRHE, CGET04, + $ ZHECON, CHERFS, CHET01, CHETRF_AASEN, ZHETRI2, + $ CHETRS_AASEN, CLACPY, CLAIPD, CLARHS, CLATB4, + $ CLATMS, CPOT02, ZPOT03, ZPOT05 +* .. +* .. Intrinsic Functions .. + INTRINSIC REAL, IMAG, MAX, MIN +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* +* Test path +* + PATH( 1: 1 ) = 'Complex precision' + PATH( 2: 3 ) = 'HA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Complex precision' + MATPATH( 2: 3 ) = 'HE' + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE +* +* Test the error exits +* + IF( TSTERR ) + $ CALL CERRHE( PATH, NOUT ) + INFOT = 0 +* +* Set the minimum block size for which the block routine should +* be used, which will be later returned by ILAENV +* + CALL XLAENV( 2, 2 ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + IF( N .GT. NMAX ) THEN + NFAIL = NFAIL + 1 + WRITE(NOUT, 9995) 'M ', N, NMAX + GO TO 180 + END IF + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + IZERO = 0 + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Set up parameters with CLATB4 for the matrix generator +* based on the type of matrix to be generated. +* + CALL CLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, + $ ANORM, MODE, CNDNUM, DIST ) +* +* Generate a matrix with CLATMS. +* + SRNAMT = 'CLATMS' + CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from CLATMS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) +* +* Skip all tests for this generated matrix +* + GO TO 160 + END IF +* +* For types 3-6, zero one or more rows and columns of +* the matrix to test that INFO is returned correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = CZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = CZERO + 50 CONTINUE + END IF + ELSE + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + IOFF = 0 + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = CZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + IOFF = 0 + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = CZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* +* End generate test matrix A. +* +* +* Set the imaginary part of the diagonals. +* + CALL CLAIPD( N, A, LDA+1, 0 ) +* +* Do for each value of NB in NBVAL +* + DO 150 INB = 1, NNB +* +* Set the optimal blocksize, which will be later +* returned by ILAENV. +* + NB = NBVAL( INB ) + CALL XLAENV( 1, NB ) +* +* Copy the test matrix A into matrix AFAC which +* will be factorized in place. This is needed to +* preserve the test matrix A for subsequent tests. +* + CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +* +* Compute the L*D*L**T or U*D*U**T factorization of the +* matrix. IWORK stores details of the interchanges and +* the block structure of D. AINV is a work array for +* block factorization, LWORK is the length of AINV. +* + LWORK = ( NB+1 )*LDA + SRNAMT = 'CHETRF_AASEN' + CALL CHETRF_AASEN( UPLO, N, AFAC, LDA, IWORK, AINV, + $ LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from CHETRF and handle error. +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'CHETRF_AASEN', INFO, K, UPLO, + $ N, N, -1, -1, NB, IMAT, NFAIL, NERRS, + $ NOUT ) + END IF +* +* Set the condition estimate flag if the INFO is not 0. +* + IF( INFO.NE.0 ) THEN + TRFCON = .TRUE. + ELSE + TRFCON = .FALSE. + END IF +* +*+ TEST 1 +* Reconstruct matrix from factors and compute residual. +* + CALL CHET01_AASEN( UPLO, N, A, LDA, AFAC, LDA, IWORK, + $ AINV, LDA, RWORK, RESULT( 1 ) ) + NT = 1 +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, + $ RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT +* +* Do only the condition estimate if INFO is not 0. +* + IF( TRFCON ) THEN + RCONDC = ZERO + GO TO 140 + END IF +* +* Do for each value of NRHS in NSVAL. +* + DO 130 IRHS = 1, NNS + NRHS = NSVAL( IRHS ) +* +*+ TEST 3 (Using TRS) +* Solve and compute residual for A * X = B. +* +* Choose a set of NRHS random solution vectors +* stored in XACT and set up the right hand side B +* + SRNAMT = 'CLARHS' + CALL CLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, + $ KL, KU, NRHS, A, LDA, XACT, LDA, + $ B, LDA, ISEED, INFO ) + CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* + SRNAMT = 'CHETRS_AASEN' + LWORK = 3*N-2 + CALL CHETRS_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK, + $ X, LDA, WORK, LWORK, INFO ) +* +* Check error code from CHETRS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'CHETRS_AASEN', INFO, 0, + $ UPLO, N, N, -1, -1, NRHS, IMAT, + $ NFAIL, NERRS, NOUT ) + END IF +* + CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) +* +* Compute the residual for the solution +* + CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* Print information about the tests that did not pass +* the threshold. +* + DO 120 K = 2, 2 + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, + $ IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 120 CONTINUE + NRUN = NRUN + 1 +* +* End do for each value of NRHS in NSVAL. +* + 130 CONTINUE + 140 CONTINUE + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ', + $ I2, ', test ', I2, ', ratio =', G12.5 ) + 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', + $ I2, ', test(', I2, ') =', G12.5 ) + 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=', + $ I6 ) + RETURN +* +* End of CCHKHE_AASEN +* + END diff --git a/TESTING/LIN/cdrvhe_aasen.f b/TESTING/LIN/cdrvhe_aasen.f new file mode 100644 index 00000000..617f307a --- /dev/null +++ b/TESTING/LIN/cdrvhe_aasen.f @@ -0,0 +1,529 @@ +*> \brief \b CDRVHE_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE CDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NOUT, NRHS +* REAL THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NVAL( * ) +* REAL RWORK( * ) +* COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), +* $ WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CDRVHE_AASEN tests the driver routine CHESV_AASEN. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand side vectors to be generated for +*> each linear system. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is REAL +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is COMPLEX array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is COMPLEX array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is COMPLEX array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is COMPLEX array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is COMPLEX array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is COMPLEX array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (NMAX*max(2,NRHS)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (NMAX+2*NRHS) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup complex_lin +* +* ===================================================================== + SUBROUTINE CDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, + $ RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NMAX, NN, NOUT, NRHS + REAL THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NVAL( * ) + REAL RWORK( * ) + COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), + $ WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) + INTEGER NTYPES, NTESTS + PARAMETER ( NTYPES = 10, NTESTS = 3 ) + INTEGER NFACT + PARAMETER ( NFACT = 2 ) +* .. +* .. Local Scalars .. + LOGICAL ZEROT + CHARACTER DIST, FACT, TYPE, UPLO, XTYPE + CHARACTER*3 MATPATH, PATH + INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, + $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N, + $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT + REAL AINVNM, ANORM, CNDNUM, RCOND, RCONDC +* .. +* .. Local Arrays .. + CHARACTER FACTS( NFACT ), UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + REAL RESULT( NTESTS ) +* .. +* .. External Functions .. + REAL CLANHE, SGET06 + EXTERNAL CLANHE, SGET06 +* .. +* .. External Subroutines .. + EXTERNAL ALADHD, ALAERH, ALASVM, XLAENV, CERRVX, + $ CGET04, CLACPY, CLARHS, CLATB4, CLATMS, + $ CHESV_AASEN, CHET01_AASEN, CPOT02, + $ CHETRF_AASEN +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Intrinsic Functions .. + INTRINSIC CMPLX, MAX, MIN +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Complex precision' + PATH( 2: 3 ) = 'HA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Complex precision' + MATPATH( 2: 3 ) = 'HE' +* + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE + LWORK = MAX( 2*NMAX, NMAX*NRHS ) +* +* Test the error exits +* + IF( TSTERR ) + $ CALL CERRVX( PATH, NOUT ) + INFOT = 0 +* +* Set the block size and minimum block size for testing. +* + NB = 1 + NBMIN = 2 + CALL XLAENV( 1, NB ) + CALL XLAENV( 2, NBMIN ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Begin generate the test matrix A. +* +* Set up parameters with CLATB4 for the matrix generator +* based on the type of matrix to be generated. +* + CALL CLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM, + $ MODE, CNDNUM, DIST ) +* +* Generate a matrix with CLATMS. +* + SRNAMT = 'CLATMS' + CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, + $ WORK, INFO ) +* +* Check error code from CLATMS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, + $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT ) + GO TO 160 + END IF +* +* For types 3-6, zero one or more rows and columns of +* the matrix to test that INFO is returned correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = ZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = ZERO + 50 CONTINUE + END IF + ELSE + IOFF = 0 + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = ZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the first IZERO rows and columns to zero. +* + IOFF = 0 + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = ZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* +* End generate the test matrix A. +* +* + DO 150 IFACT = 1, NFACT +* +* Do first for FACT = 'F', then for other values. +* + FACT = FACTS( IFACT ) +* +* Compute the condition number for comparison with +* the value returned by CHESVX. +* + IF( ZEROT ) THEN + IF( IFACT.EQ.1 ) + $ GO TO 150 + RCONDC = ZERO +* + ELSE IF( IFACT.EQ.1 ) THEN +* +* Compute the 1-norm of A. +* + ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK ) +* +* Factor the matrix A. +* +c CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +c SRNAMT = 'CHETRF_AASEN' +c CALL CHETRF_AASEN( UPLO, N, AFAC, LDA, IWORK, +c $ WORK, LWORK, INFO ) +* +* Compute inv(A) and take its norm. +* +c CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) +c LWORK = (N+NB+1)*(NB+3) +c SRNAMT = 'CHETRI2' +c CALL CHETRI2( UPLO, N, AINV, LDA, IWORK, WORK, +c $ LWORK, INFO ) +c AINVNM = CLANHE( '1', UPLO, N, AINV, LDA, RWORK ) +* +* Compute the 1-norm condition number of A. +* +c IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN +c RCONDC = ONE +c ELSE +c RCONDC = ( ONE / ANORM ) / AINVNM +c END IF + END IF +* +* Form an exact solution and set the right hand side. +* + SRNAMT = 'CLARHS' + CALL CLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU, + $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED, + $ INFO ) + XTYPE = 'C' +* +* --- Test CHESV_AASEN --- +* + IF( IFACT.EQ.2 ) THEN + CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) + CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* +* Factor the matrix and solve the system using CHESV_AASEN. +* + SRNAMT = 'CHESV_AASEN ' + CALL CHESV_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK, + $ X, LDA, WORK, LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from CHESV_AASEN . +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'CHESV_AASEN', INFO, K, + $ UPLO, N, N, -1, -1, NRHS, + $ IMAT, NFAIL, NERRS, NOUT ) + GO TO 120 + ELSE IF( INFO.NE.0 ) THEN + GO TO 120 + END IF +* +* Reconstruct matrix from factors and compute +* residual. +* + CALL CHET01_AASEN( UPLO, N, A, LDA, AFAC, LDA, + $ IWORK, AINV, LDA, RWORK, + $ RESULT( 1 ) ) +* +* Compute residual of the computed solution. +* + CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* Check solution from generated exact solution. +* + CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, + $ RESULT( 3 ) ) + NT = 3 +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALADHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )'CHESV_AASEN ', + $ UPLO, N, IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT + 120 CONTINUE + END IF +* + 150 CONTINUE +* + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2, + $ ', test ', I2, ', ratio =', G12.5 ) + RETURN +* +* End of CDRVHE_AASEN +* + END diff --git a/TESTING/LIN/cerrvx.f b/TESTING/LIN/cerrvx.f index 52ca890d..380dc699 100644 --- a/TESTING/LIN/cerrvx.f +++ b/TESTING/LIN/cerrvx.f @@ -93,7 +93,7 @@ $ CHESV, CHESV_ROOK, CHESVX, CHKXER, CHPSV, $ CHPSVX, CPBSV, CPBSVX, CPOSV, CPOSVX, CPPSV, $ CPPSVX, CPTSV, CPTSVX, CSPSV, CSPSVX, CSYSV, - $ CSYSV_ROOK, CSYSVX + $ CSYSV_AASEN, CSYSV_ROOK, CSYSVX * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -632,6 +632,25 @@ $ RCOND, R1, R2, W, 3, RW, INFO ) CALL CHKXER( 'CHESVX', INFOT, NOUT, LERR, OK ) * + ELSE IF( LSAMEN( 2, C2, 'HA' ) ) THEN +* +* CHESV_AASEN +* + SRNAMT = 'CHESV_AASEN' + INFOT = 1 + CALL CHESV_AASEN( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'CHESV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL CHESV_AASEN( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'CHESV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL CHESV_AASEN( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'CHESV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL CHESV_AASEN( 'U', 2, 0, A, 2, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'CHESV_AASEN', INFOT, NOUT, LERR, OK ) +* + ELSE IF( LSAMEN( 2, C2, 'HR' ) ) THEN * * CHESV_ROOK diff --git a/TESTING/LIN/chet01_aasen.f b/TESTING/LIN/chet01_aasen.f new file mode 100644 index 00000000..d87a6102 --- /dev/null +++ b/TESTING/LIN/chet01_aasen.f @@ -0,0 +1,267 @@ +*> \brief \b CHET01_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE CHET01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, +* C, LDC, RWORK, RESID ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER LDA, LDAFAC, LDC, N +* COMPLEX RESID +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), +* $ RWORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CHET01_AASEN reconstructs a hermitian indefinite matrix A from its +*> block L*D*L' or U*D*U' factorization and computes the residual +*> norm( C - A ) / ( N * norm(A) * EPS ), +*> where C is the reconstructed matrix and EPS is the machine epsilon. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the upper or lower triangular part of the +*> hermitian matrix A is stored: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of rows and columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> The original hermitian matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N) +*> \endverbatim +*> +*> \param[in] AFAC +*> \verbatim +*> AFAC is COMPLEX array, dimension (LDAFAC,N) +*> The factored form of the matrix A. AFAC contains the block +*> diagonal matrix D and the multipliers used to obtain the +*> factor L or U from the block L*D*L' or U*D*U' factorization +*> as computed by CHETRF. +*> \endverbatim +*> +*> \param[in] LDAFAC +*> \verbatim +*> LDAFAC is INTEGER +*> The leading dimension of the array AFAC. LDAFAC >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> The pivot indices from CHETRF. +*> \endverbatim +*> +*> \param[out] C +*> \verbatim +*> C is COMPLEX array, dimension (LDC,N) +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is COMPLEX array, dimension (N) +*> \endverbatim +*> +*> \param[out] RESID +*> \verbatim +*> RESID is COMPLEX +*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) +*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* +*> \ingroup complex_lin +* +* ===================================================================== + SUBROUTINE CHET01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, + $ LDC, RWORK, RESID ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDAFAC, LDC, N + REAL RESID +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), + $ RWORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX CZERO, CONE + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), + $ CONE = ( 1.0E+0, 0.0E+0 ) ) + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J + REAL ANORM, EPS +* .. +* .. External Functions .. + LOGICAL LSAME + REAL SLAMCH, CLANHE + EXTERNAL LSAME, SLAMCH, CLANHE +* .. +* .. External Subroutines .. + EXTERNAL CLASET, CLAVHE +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE +* .. +* .. Executable Statements .. +* +* Quick exit if N = 0. +* + IF( N.LE.0 ) THEN + RESID = ZERO + RETURN + END IF +* +* Determine EPS and the norm of A. +* + EPS = SLAMCH( 'Epsilon' ) + ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK ) +* +* Initialize C to the tridiagonal matrix T. +* + CALL CLASET( 'Full', N, N, CZERO, CZERO, C, LDC ) + CALL CLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 ) + IF( N.GT.1 ) THEN + IF( LSAME( UPLO, 'U' ) ) THEN + CALL CLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL CLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + CALL CLACGV( N-1, C( 2, 1 ), LDC+1 ) + ELSE + CALL CLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL CLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + CALL CLACGV( N-1, C( 1, 2 ), LDC+1 ) + ENDIF + ENDIF +* +* Call CTRMM to form the product U' * D (or L * D ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL CTRMM( 'Left', UPLO, 'Conjugate transpose', 'Unit', N-1, + $ N, CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) + ELSE + CALL CTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF +* +* Call CTRMM again to multiply by U (or L ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL CTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL CTRMM( 'Right', UPLO, 'Conjugate transpose', 'Unit', N, + $ N-1, CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) + END IF +* +* Apply hermitian pivots +* + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL CSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC ) + END DO + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL CSWAP( N, C( 1, J ), 1, C( 1, I ), 1 ) + END DO +* +* +* Compute the difference C - A . +* + IF( LSAME( UPLO, 'U' ) ) THEN + DO J = 1, N + DO I = 1, J + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + ELSE + DO J = 1, N + DO I = J, N + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + END IF +* +* Compute norm( C - A ) / ( N * norm(A) * EPS ) +* + RESID = CLANHE( '1', UPLO, N, C, LDC, RWORK ) +* + IF( ANORM.LE.ZERO ) THEN + IF( RESID.NE.ZERO ) + $ RESID = ONE / EPS + ELSE + RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS + END IF +* + RETURN +* +* End of CHET01_AASEN +* + END diff --git a/TESTING/LIN/dchkaa.f b/TESTING/LIN/dchkaa.f index 49478c93..9e7a14ab 100644 --- a/TESTING/LIN/dchkaa.f +++ b/TESTING/LIN/dchkaa.f @@ -49,6 +49,7 @@ *> DPP 9 List types on next line if 0 < NTYPES < 9 *> DPB 8 List types on next line if 0 < NTYPES < 8 *> DPT 12 List types on next line if 0 < NTYPES < 12 +*> DSA 10 List types on next line if 0 < NTYPES < 10 *> DSY 10 List types on next line if 0 < NTYPES < 10 *> DSR 10 List types on next line if 0 < NTYPES < 10 *> DSP 10 List types on next line if 0 < NTYPES < 10 @@ -663,6 +664,34 @@ ELSE WRITE( NOUT, FMT = 9988 )PATH END IF +* + ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN +* +* SY: symmetric indefinite matrices, +* with partial (Aasen's) pivoting algorithm +* + NTYPES = 10 + CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) +* + IF( TSTCHK ) THEN + CALL DCHKSY_AASEN( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS, + $ NSVAL, THRESH, TSTERR, LDA, + $ A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), + $ WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9989 )PATH + END IF +* + IF( TSTDRV ) THEN + CALL DDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), + $ WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9988 )PATH + END IF +* * ELSE IF( LSAMEN( 2, C2, 'SP' ) ) THEN * diff --git a/TESTING/LIN/dchksy_aasen.f b/TESTING/LIN/dchksy_aasen.f new file mode 100644 index 00000000..041ef754 --- /dev/null +++ b/TESTING/LIN/dchksy_aasen.f @@ -0,0 +1,577 @@ +*> \brief \b DCHKSY_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DCHKSY_AASEN( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, +* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, +* XACT, WORK, RWORK, IWORK, NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NNB, NNS, NOUT +* DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) +* DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ), +* $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DCHKSY_AASEN tests DSYTRF_AASEN, -TRS_AASEN. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NNB +*> \verbatim +*> NNB is INTEGER +*> The number of values of NB contained in the vector NBVAL. +*> \endverbatim +*> +*> \param[in] NBVAL +*> \verbatim +*> NBVAL is INTEGER array, dimension (NBVAL) +*> The values of the blocksize NB. +*> \endverbatim +*> +*> \param[in] NNS +*> \verbatim +*> NNS is INTEGER +*> The number of values of NRHS contained in the vector NSVAL. +*> \endverbatim +*> +*> \param[in] NSVAL +*> \verbatim +*> NSVAL is INTEGER array, dimension (NNS) +*> The values of the number of right hand sides NRHS. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is DOUBLE PRECISION +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is DOUBLE PRECISION array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is DOUBLE PRECISION array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is DOUBLE PRECISION array, dimension (NMAX*NSMAX) +*> where NSMAX is the largest entry in NSVAL. +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is DOUBLE PRECISION array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is DOUBLE PRECISION array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (NMAX*max(3,NSMAX)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (max(NMAX,2*NSMAX)) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* +*> \ingroup double_lin +* +* ===================================================================== + SUBROUTINE DCHKSY_AASEN( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, + $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, + $ X, XACT, WORK, RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NN, NNB, NNS, NMAX, NOUT + DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) + DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ), + $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + INTEGER NTYPES + PARAMETER ( NTYPES = 10 ) + INTEGER NTESTS + PARAMETER ( NTESTS = 9 ) +* .. +* .. Local Scalars .. + LOGICAL TRFCON, ZEROT + CHARACTER DIST, TYPE, UPLO, XTYPE + CHARACTER*3 PATH, MATPATH + INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, + $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, + $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT + DOUBLE PRECISION ANORM, CNDNUM, RCONDC +* .. +* .. Local Arrays .. + CHARACTER UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + DOUBLE PRECISION RESULT( NTESTS ) +* .. +* .. External Functions .. + DOUBLE PRECISION DGET06, DLANSY + EXTERNAL DGET06, DLANSY +* .. +* .. External Subroutines .. + EXTERNAL ALAERH, ALAHD, ALASUM, DERRSY, DGET04, DLACPY, + $ DLARHS, DLATB4, DLATMS, DPOT02, DPOT03, DPOT05, + $ DSYCON, DSYRFS, DSYT01, DSYTRF_AASEN, + $ DSYTRI2, DSYTRS_AASEN, XLAENV +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Double precision' + PATH( 2: 3 ) = 'SA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Double precision' + MATPATH( 2: 3 ) = 'SY' + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE +* +* Test the error exits +* + IF( TSTERR ) + $ CALL DERRSY( PATH, NOUT ) + INFOT = 0 +* +* Set the minimum block size for which the block routine should +* be used, which will be later returned by ILAENV +* + CALL XLAENV( 2, 2 ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + IF( N .GT. NMAX ) THEN + NFAIL = NFAIL + 1 + WRITE(NOUT, 9995) 'M ', N, NMAX + GO TO 180 + END IF + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + IZERO = 0 +* +* Do for each value of matrix type IMAT +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Begin generate the test matrix A. +* +* +* Set up parameters with DLATB4 for the matrix generator +* based on the type of matrix to be generated. +* + CALL DLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, + $ ANORM, MODE, CNDNUM, DIST ) +* +* Generate a matrix with DLATMS. +* + SRNAMT = 'DLATMS' + CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from DLATMS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) +* +* Skip all tests for this generated matrix +* + GO TO 160 + END IF +* +* For matrix types 3-6, zero one or more rows and +* columns of the matrix to test that INFO is returned +* correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = ZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = ZERO + 50 CONTINUE + END IF + ELSE + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + IOFF = 0 + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = ZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + IOFF = 0 + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = ZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* +* End generate the test matrix A. +* +* Do for each value of NB in NBVAL +* + DO 150 INB = 1, NNB +* +* Set the optimal blocksize, which will be later +* returned by ILAENV. +* + NB = NBVAL( INB ) + CALL XLAENV( 1, NB ) +* +* Copy the test matrix A into matrix AFAC which +* will be factorized in place. This is needed to +* preserve the test matrix A for subsequent tests. +* + CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +* +* Compute the L*D*L**T or U*D*U**T factorization of the +* matrix. IWORK stores details of the interchanges and +* the block structure of D. AINV is a work array for +* block factorization, LWORK is the length of AINV. +* + SRNAMT = 'DSYTRF_AASEN' + LWORK = N*NB + N + CALL DSYTRF_AASEN( UPLO, N, AFAC, LDA, IWORK, AINV, + $ LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from DSYTRF and handle error. +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'DSYTRF_AASEN', INFO, K, UPLO, + $ N, N, -1, -1, NB, IMAT, NFAIL, NERRS, + $ NOUT ) + END IF +* +* Set the condition estimate flag if the INFO is not 0. +* + IF( INFO.NE.0 ) THEN + TRFCON = .TRUE. + ELSE + TRFCON = .FALSE. + END IF +* +*+ TEST 1 +* Reconstruct matrix from factors and compute residual. +* + CALL DSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDA, IWORK, + $ AINV, LDA, RWORK, RESULT( 1 ) ) + NT = 1 +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, + $ RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT +* +* Do only the condition estimate if INFO is not 0. +* + IF( TRFCON ) THEN + RCONDC = ZERO + GO TO 140 + END IF +* +* Do for each value of NRHS in NSVAL. +* + DO 130 IRHS = 1, NNS + NRHS = NSVAL( IRHS ) +* +*+ TEST 3 ( Using TRS) +* Solve and compute residual for A * X = B. +* +* Choose a set of NRHS random solution vectors +* stored in XACT and set up the right hand side B +* + SRNAMT = 'DLARHS' + CALL DLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, + $ KL, KU, NRHS, A, LDA, XACT, LDA, + $ B, LDA, ISEED, INFO ) + CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* + SRNAMT = 'DSYTRS_AASEN' + LWORK = 3*N-2 + CALL DSYTRS_AASEN( UPLO, N, NRHS, AFAC, LDA, + $ IWORK, X, LDA, WORK, LWORK, + $ INFO ) +* +* Check error code from DSYTRS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'DSYTRS_AASEN', INFO, 0, + $ UPLO, N, N, -1, -1, NRHS, IMAT, + $ NFAIL, NERRS, NOUT ) + END IF +* + CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) +* +* Compute the residual for the solution +* + CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 120 K = 2, 2 + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, + $ IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 120 CONTINUE + NRUN = NRUN + 1 +* +* End do for each value of NRHS in NSVAL. +* + 130 CONTINUE + 140 CONTINUE + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ', + $ I2, ', test ', I2, ', ratio =', G12.5 ) + 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', + $ I2, ', test(', I2, ') =', G12.5 ) + 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=', + $ I6 ) + RETURN +* +* End of DCHKSY_AASEN +* + END diff --git a/TESTING/LIN/ddrvsy_aasen.f b/TESTING/LIN/ddrvsy_aasen.f new file mode 100644 index 00000000..a3520cb1 --- /dev/null +++ b/TESTING/LIN/ddrvsy_aasen.f @@ -0,0 +1,517 @@ +*> \brief \b DDRVSY_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NOUT, NRHS +* DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NVAL( * ) +* DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ), +* $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DDRVSY_AASEN tests the driver routine DSYSV_AASEN. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand side vectors to be generated for +*> each linear system. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is DOUBLE PRECISION +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is DOUBLE PRECISION array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is DOUBLE PRECISION array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is DOUBLE PRECISION array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is DOUBLE PRECISION array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is DOUBLE PRECISION array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (NMAX*max(2,NRHS)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (NMAX+2*NRHS) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup double_lin +* +* ===================================================================== + SUBROUTINE DDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, + $ RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NMAX, NN, NOUT, NRHS + DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NVAL( * ) + DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ), + $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + INTEGER NTYPES, NTESTS + PARAMETER ( NTYPES = 10, NTESTS = 3 ) + INTEGER NFACT + PARAMETER ( NFACT = 2 ) +* .. +* .. Local Scalars .. + LOGICAL ZEROT + CHARACTER DIST, FACT, TYPE, UPLO, XTYPE + CHARACTER*3 MATPATH, PATH + INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, + $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N, + $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT + DOUBLE PRECISION AINVNM, ANORM, CNDNUM, RCOND, RCONDC +* .. +* .. Local Arrays .. + CHARACTER FACTS( NFACT ), UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + DOUBLE PRECISION RESULT( NTESTS ) +* .. +* .. External Functions .. + DOUBLE PRECISION DGET06, DLANSY + EXTERNAL DGET06, DLANSY +* .. +* .. External Subroutines .. + EXTERNAL ALADHD, ALAERH, ALASVM, DERRVX, DGET04, DLACPY, + $ DLARHS, DLASET, DLATB4, DLATMS, DPOT02, DPOT05, + $ DSYSV_AASEN, DSYT01_AASEN, DSYTRF_AASEN, XLAENV +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Double precision' + PATH( 2: 3 ) = 'SA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Double precision' + MATPATH( 2: 3 ) = 'SY' +* + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE + LWORK = MAX( 2*NMAX, NMAX*NRHS ) +* +* Test the error exits +* + IF( TSTERR ) + $ CALL DERRVX( PATH, NOUT ) + INFOT = 0 +* +* Set the block size and minimum block size for testing. +* + NB = 1 + NBMIN = 2 + CALL XLAENV( 1, NB ) + CALL XLAENV( 2, NBMIN ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Set up parameters with DLATB4 and generate a test matrix +* with DLATMS. +* + CALL DLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM, + $ MODE, CNDNUM, DIST ) +* + SRNAMT = 'DLATMS' + CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from DLATMS. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) + GO TO 160 + END IF +* +* For types 3-6, zero one or more rows and columns of the +* matrix to test that INFO is returned correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = ZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = ZERO + 50 CONTINUE + END IF + ELSE + IOFF = 0 + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = ZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = ZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* + DO 150 IFACT = 1, NFACT +* +* Do first for FACT = 'F', then for other values. +* + FACT = FACTS( IFACT ) +* +* Compute the condition number for comparison with +* the value returned by DSYSVX. +* + IF( ZEROT ) THEN + IF( IFACT.EQ.1 ) + $ GO TO 150 + RCONDC = ZERO +* + ELSE IF( IFACT.EQ.1 ) THEN +* +* Compute the 1-norm of A. +* + ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) +* +* Factor the matrix A. +* +c CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +c CALL DSYTRF( UPLO, N, AFAC, LDA, IWORK, WORK, +c $ LWORK, INFO ) +* +* Compute inv(A) and take its norm. +* +c CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) +c LWORK = (N+NB+1)*(NB+3) +c SRNAMT = 'DSYTRI2' +c CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK, +c $ LWORK, INFO ) +c AINVNM = DLANSY( '1', UPLO, N, AINV, LDA, RWORK ) +* +* Compute the 1-norm condition number of A. +* +c IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN +c RCONDC = ONE +c ELSE +c RCONDC = ( ONE / ANORM ) / AINVNM +c END IF + END IF +* +* Form an exact solution and set the right hand side. +* + SRNAMT = 'DLARHS' + CALL DLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU, + $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED, + $ INFO ) + XTYPE = 'C' +* +* --- Test DSYSV_AASEN --- +* + IF( IFACT.EQ.2 ) THEN + CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) + CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* +* Factor the matrix and solve the system using DSYSV_AASEN. +* + SRNAMT = 'DSYSV_AASEN' + CALL DSYSV_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK, + $ X, LDA, WORK, LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from DSYSV_AASEN . +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'DSYSV_AASEN ', INFO, K, + $ UPLO, N, N, -1, -1, NRHS, + $ IMAT, NFAIL, NERRS, NOUT ) + GO TO 120 + ELSE IF( INFO.NE.0 ) THEN + GO TO 120 + END IF +* +* Reconstruct matrix from factors and compute +* residual. +* + CALL DSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDA, + $ IWORK, AINV, LDA, RWORK, + $ RESULT( 1 ) ) +* +* Compute residual of the computed solution. +* + CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* Check solution from generated exact solution. +* + CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, + $ RESULT( 3 ) ) + NT = 3 +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALADHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )'DSYSV_AASEN ', + $ UPLO, N, IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT + 120 CONTINUE + END IF +* + 150 CONTINUE +* + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2, + $ ', test ', I2, ', ratio =', G12.5 ) + RETURN +* +* End of DDRVSY_AASEN +* + END diff --git a/TESTING/LIN/derrvx.f b/TESTING/LIN/derrvx.f index 36eb8d6a..2affbe2f 100644 --- a/TESTING/LIN/derrvx.f +++ b/TESTING/LIN/derrvx.f @@ -91,7 +91,7 @@ EXTERNAL CHKXER, DGBSV, DGBSVX, DGESV, DGESVX, DGTSV, $ DGTSVX, DPBSV, DPBSVX, DPOSV, DPOSVX, DPPSV, $ DPPSVX, DPTSV, DPTSVX, DSPSV, DSPSVX, DSYSV, - $ DSYSV_ROOK, DSYSVX + $ DSYSV_AASEN, DSYSV_ROOK, DSYSVX * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -627,6 +627,25 @@ $ RCOND, R1, R2, W, 3, IW, INFO ) CALL CHKXER( 'DSYSVX', INFOT, NOUT, LERR, OK ) * + ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN +* +* DSYSV_AASEN +* + SRNAMT = 'DSYSV_AASEN' + INFOT = 1 + CALL DSYSV_AASEN( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'DSYSV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL DSYSV_AASEN( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'DSYSV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL DSYSV_AASEN( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'DSYSV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL DSYSV_AASEN( 'U', 2, 0, A, 2, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'DSYSV_AASEN', INFOT, NOUT, LERR, OK ) +* + ELSE IF( LSAMEN( 2, C2, 'SR' ) ) THEN * * DSYSV_ROOK diff --git a/TESTING/LIN/dsyt01_aasen.f b/TESTING/LIN/dsyt01_aasen.f new file mode 100644 index 00000000..f3756771 --- /dev/null +++ b/TESTING/LIN/dsyt01_aasen.f @@ -0,0 +1,263 @@ +*> \brief \b DSYT01 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DSYT01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC, +* RWORK, RESID ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER LDA, LDAFAC, LDC, N +* DOUBLE PRECISION RESID +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), +* $ RWORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DSYT01 reconstructs a symmetric indefinite matrix A from its +*> block L*D*L' or U*D*U' factorization and computes the residual +*> norm( C - A ) / ( N * norm(A) * EPS ), +*> where C is the reconstructed matrix and EPS is the machine epsilon. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the upper or lower triangular part of the +*> symmetric matrix A is stored: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of rows and columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> The original symmetric matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N) +*> \endverbatim +*> +*> \param[in] AFAC +*> \verbatim +*> AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N) +*> The factored form of the matrix A. AFAC contains the block +*> diagonal matrix D and the multipliers used to obtain the +*> factor L or U from the block L*D*L' or U*D*U' factorization +*> as computed by DSYTRF. +*> \endverbatim +*> +*> \param[in] LDAFAC +*> \verbatim +*> LDAFAC is INTEGER +*> The leading dimension of the array AFAC. LDAFAC >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> The pivot indices from DSYTRF. +*> \endverbatim +*> +*> \param[out] C +*> \verbatim +*> C is DOUBLE PRECISION array, dimension (LDC,N) +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (N) +*> \endverbatim +*> +*> \param[out] RESID +*> \verbatim +*> RESID is DOUBLE PRECISION +*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) +*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* @precisions fortran d -> s +* +*> \ingroup double_lin +* +* ===================================================================== + SUBROUTINE DSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, + $ LDC, RWORK, RESID ) +* +* -- LAPACK test routine (version 3.5.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDAFAC, LDC, N + DOUBLE PRECISION RESID +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), + $ RWORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J + DOUBLE PRECISION ANORM, EPS +* .. +* .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, DLANSY + EXTERNAL LSAME, DLAMCH, DLANSY +* .. +* .. External Subroutines .. + EXTERNAL DLASET, DLAVSY +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE +* .. +* .. Executable Statements .. +* +* Quick exit if N = 0. +* + IF( N.LE.0 ) THEN + RESID = ZERO + RETURN + END IF +* +* Determine EPS and the norm of A. +* + EPS = DLAMCH( 'Epsilon' ) + ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) +* +* Initialize C to the tridiagonal matrix T. +* + CALL DLASET( 'Full', N, N, ZERO, ZERO, C, LDC ) + CALL DLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 ) + IF( N.GT.1 ) THEN + IF( LSAME( UPLO, 'U' ) ) THEN + CALL DLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL DLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + ELSE + CALL DLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL DLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + ENDIF + ENDIF +* +* Call DTRMM to form the product U' * D (or L * D ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL DTRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N, + $ ONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) + ELSE + CALL DTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ ONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF +* +* Call DTRMM again to multiply by U (or L ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL DTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ ONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL DTRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1, + $ ONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) + END IF +* +* Apply symmetric pivots +* + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL DSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC ) + END DO + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL DSWAP( N, C( 1, J ), 1, C( 1, I ), 1 ) + END DO +* +* +* Compute the difference C - A . +* + IF( LSAME( UPLO, 'U' ) ) THEN + DO J = 1, N + DO I = 1, J + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + ELSE + DO J = 1, N + DO I = J, N + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + END IF +* +* Compute norm( C - A ) / ( N * norm(A) * EPS ) +* + RESID = DLANSY( '1', UPLO, N, C, LDC, RWORK ) +* + IF( ANORM.LE.ZERO ) THEN + IF( RESID.NE.ZERO ) + $ RESID = ONE / EPS + ELSE + RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS + END IF +* + RETURN +* +* End of DSYT01 +* + END diff --git a/TESTING/LIN/schkaa.f b/TESTING/LIN/schkaa.f index 821710ce..a6bfa101 100644 --- a/TESTING/LIN/schkaa.f +++ b/TESTING/LIN/schkaa.f @@ -663,6 +663,33 @@ ELSE WRITE( NOUT, FMT = 9988 )PATH END IF +* + ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN +* +* SY: symmetric indefinite matrices, +* with partial (Aasen's) pivoting algorithm +* + NTYPES = 10 + CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) +* + IF( TSTCHK ) THEN + CALL SCHKSY_AASEN( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS, + $ NSVAL, THRESH, TSTERR, LDA, + $ A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), + $ WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9989 )PATH + END IF +* + IF( TSTDRV ) THEN + CALL SDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), + $ WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9988 )PATH + END IF * ELSE IF( LSAMEN( 2, C2, 'SP' ) ) THEN * diff --git a/TESTING/LIN/schksy_aasen.f b/TESTING/LIN/schksy_aasen.f new file mode 100644 index 00000000..659e3fd3 --- /dev/null +++ b/TESTING/LIN/schksy_aasen.f @@ -0,0 +1,578 @@ +*> \brief \b SCHKSY_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE SCHKSY_AAEN( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, +* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, +* XACT, WORK, RWORK, IWORK, NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NNB, NNS, NOUT +* REAL THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) +* REAL A( * ), AFAC( * ), AINV( * ), B( * ), +* $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SCHKSY_AASEN tests SSYTRF_AASEN, -TRS_AASEN. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NNB +*> \verbatim +*> NNB is INTEGER +*> The number of values of NB contained in the vector NBVAL. +*> \endverbatim +*> +*> \param[in] NBVAL +*> \verbatim +*> NBVAL is INTEGER array, dimension (NBVAL) +*> The values of the blocksize NB. +*> \endverbatim +*> +*> \param[in] NNS +*> \verbatim +*> NNS is INTEGER +*> The number of values of NRHS contained in the vector NSVAL. +*> \endverbatim +*> +*> \param[in] NSVAL +*> \verbatim +*> NSVAL is INTEGER array, dimension (NNS) +*> The values of the number of right hand sides NRHS. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is REAL +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is REAL array, dimension (NMAX*NSMAX) +*> where NSMAX is the largest entry in NSVAL. +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is REAL array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is REAL array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension (NMAX*max(3,NSMAX)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (max(NMAX,2*NSMAX)) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* +*> \ingroup real_lin +* +* ===================================================================== + SUBROUTINE SCHKSY_AASEN( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, + $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, + $ X, XACT, WORK, RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NN, NNB, NNS, NMAX, NOUT + REAL THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) + REAL A( * ), AFAC( * ), AINV( * ), B( * ), + $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO + PARAMETER ( ZERO = 0.0E+0 ) + INTEGER NTYPES + PARAMETER ( NTYPES = 10 ) + INTEGER NTESTS + PARAMETER ( NTESTS = 9 ) +* .. +* .. Local Scalars .. + LOGICAL TRFCON, ZEROT + CHARACTER DIST, TYPE, UPLO, XTYPE + CHARACTER*3 PATH, MATPATH + INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, + $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, + $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT + REAL ANORM, CNDNUM, RCONDC +* .. +* .. Local Arrays .. + CHARACTER UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + REAL RESULT( NTESTS ) +* .. +* .. External Functions .. + REAL DGET06, SLANSY + EXTERNAL DGET06, SLANSY +* .. +* .. External Subroutines .. + EXTERNAL ALAERH, ALAHD, ALASUM, SERRSY, SGET04, SLACPY, + $ SLARHS, SLATB4, SLATMS, SPOT02, DPOT03, DPOT05, + $ DSYCON, SSYRFS, SSYT01_AASEN, SSYTRF_AASEN, + $ DSYTRI2, SSYTRS_AASEN, XLAENV +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* +* Test path +* + PATH( 1: 1 ) = 'Single precision' + PATH( 2: 3 ) = 'SA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Single precision' + MATPATH( 2: 3 ) = 'SY' + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE +* +* Test the error exits +* + IF( TSTERR ) + $ CALL SERRSY( PATH, NOUT ) + INFOT = 0 +* +* Set the minimum block size for which the block routine should +* be used, which will be later returned by ILAENV +* + CALL XLAENV( 2, 2 ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + IF( N .GT. NMAX ) THEN + NFAIL = NFAIL + 1 + WRITE(NOUT, 9995) 'M ', N, NMAX + GO TO 180 + END IF + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + IZERO = 0 +* +* Do for each value of matrix type IMAT +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Begin generate the test matrix A. +* +* +* Set up parameters with SLATB4 for the matrix generator +* based on the type of matrix to be generated. +* + CALL SLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, + $ ANORM, MODE, CNDNUM, DIST ) +* +* Generate a matrix with SLATMS. +* + SRNAMT = 'SLATMS' + CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from SLATMS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'SLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) +* +* Skip all tests for this generated matrix +* + GO TO 160 + END IF +* +* For matrix types 3-6, zero one or more rows and +* columns of the matrix to test that INFO is returned +* correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = ZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = ZERO + 50 CONTINUE + END IF + ELSE + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + IOFF = 0 + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = ZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + IOFF = 0 + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = ZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* +* End generate the test matrix A. +* +* Do for each value of NB in NBVAL +* + DO 150 INB = 1, NNB +* +* Set the optimal blocksize, which will be later +* returned by ILAENV. +* + NB = NBVAL( INB ) + CALL XLAENV( 1, NB ) +* +* Copy the test matrix A into matrix AFAC which +* will be factorized in place. This is needed to +* preserve the test matrix A for subsequent tests. +* + CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +* +* Compute the L*D*L**T or U*D*U**T factorization of the +* matrix. IWORK stores details of the interchanges and +* the block structure of D. AINV is a work array for +* block factorization, LWORK is the length of AINV. +* + SRNAMT = 'SSYTRF_AASEN' + LWORK = N*NB + N + CALL SSYTRF_AASEN( UPLO, N, AFAC, LDA, IWORK, AINV, + $ LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from SSYTRF and handle error. +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'SSYTRF_AASEN', INFO, K, UPLO, + $ N, N, -1, -1, NB, IMAT, NFAIL, NERRS, + $ NOUT ) + END IF +* +* Set the condition estimate flag if the INFO is not 0. +* + IF( INFO.NE.0 ) THEN + TRFCON = .TRUE. + ELSE + TRFCON = .FALSE. + END IF +* +*+ TEST 1 +* Reconstruct matrix from factors and compute residual. +* + CALL SSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDA, IWORK, + $ AINV, LDA, RWORK, RESULT( 1 ) ) + NT = 1 +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, + $ RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT +* +* Do only the condition estimate if INFO is not 0. +* + IF( TRFCON ) THEN + RCONDC = ZERO + GO TO 140 + END IF +* +* Do for each value of NRHS in NSVAL. +* + DO 130 IRHS = 1, NNS + NRHS = NSVAL( IRHS ) +* +*+ TEST 3 ( Using TRS) +* Solve and compute residual for A * X = B. +* +* Choose a set of NRHS random solution vectors +* stored in XACT and set up the right hand side B +* + SRNAMT = 'SLARHS' + CALL SLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, + $ KL, KU, NRHS, A, LDA, XACT, LDA, + $ B, LDA, ISEED, INFO ) + CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* + SRNAMT = 'SSYTRS_AASEN' + LWORK = 3*N-2 + CALL SSYTRS_AASEN( UPLO, N, NRHS, AFAC, LDA, + $ IWORK, X, LDA, WORK, LWORK, + $ INFO ) +* +* Check error code from SSYTRS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'SSYTRS_AASEN', INFO, 0, + $ UPLO, N, N, -1, -1, NRHS, IMAT, + $ NFAIL, NERRS, NOUT ) + END IF +* + CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) +* +* Compute the residual for the solution +* + CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 120 K = 2, 2 + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, + $ IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 120 CONTINUE + NRUN = NRUN + 1 +* +* End do for each value of NRHS in NSVAL. +* + 130 CONTINUE + 140 CONTINUE + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ', + $ I2, ', test ', I2, ', ratio =', G12.5 ) + 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', + $ I2, ', test(', I2, ') =', G12.5 ) + 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=', + $ I6 ) + RETURN +* +* End of SCHKSY_AASEN +* + END diff --git a/TESTING/LIN/sdrvsy_aasen.f b/TESTING/LIN/sdrvsy_aasen.f new file mode 100644 index 00000000..afa5dc49 --- /dev/null +++ b/TESTING/LIN/sdrvsy_aasen.f @@ -0,0 +1,517 @@ +*> \brief \b SDRVSY_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE SDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NOUT, NRHS +* REAL THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NVAL( * ) +* REAL A( * ), AFAC( * ), AINV( * ), B( * ), +* $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SDRVSY_AASEN tests the driver routine SSYSV_AASEN. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand side vectors to be generated for +*> each linear system. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is REAL +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is REAL array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is REAL array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is REAL array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension (NMAX*max(2,NRHS)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (NMAX+2*NRHS) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup real_lin +* +* ===================================================================== + SUBROUTINE SDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, + $ RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NMAX, NN, NOUT, NRHS + REAL THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NVAL( * ) + REAL A( * ), AFAC( * ), AINV( * ), B( * ), + $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) + INTEGER NTYPES, NTESTS + PARAMETER ( NTYPES = 10, NTESTS = 3 ) + INTEGER NFACT + PARAMETER ( NFACT = 2 ) +* .. +* .. Local Scalars .. + LOGICAL ZEROT + CHARACTER DIST, FACT, TYPE, UPLO, XTYPE + CHARACTER*3 MATPATH, PATH + INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, + $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N, + $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT + REAL AINVNM, ANORM, CNDNUM, RCOND, RCONDC +* .. +* .. Local Arrays .. + CHARACTER FACTS( NFACT ), UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + REAL RESULT( NTESTS ) +* .. +* .. External Functions .. + REAL DGET06, SLANSY + EXTERNAL DGET06, SLANSY +* .. +* .. External Subroutines .. + EXTERNAL ALADHD, ALAERH, ALASVM, SERRVX, SGET04, SLACPY, + $ SLARHS, SLASET, SLATB4, SLATMS, SPOT02, DPOT05, + $ SSYSV_AASEN, SSYT01_AASEN, SSYTRF_AASEN, XLAENV +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Single precision' + PATH( 2: 3 ) = 'SA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Single precision' + MATPATH( 2: 3 ) = 'SY' +* + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE + LWORK = MAX( 2*NMAX, NMAX*NRHS ) +* +* Test the error exits +* + IF( TSTERR ) + $ CALL SERRVX( PATH, NOUT ) + INFOT = 0 +* +* Set the block size and minimum block size for testing. +* + NB = 1 + NBMIN = 2 + CALL XLAENV( 1, NB ) + CALL XLAENV( 2, NBMIN ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Set up parameters with SLATB4 and generate a test matrix +* with SLATMS. +* + CALL SLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM, + $ MODE, CNDNUM, DIST ) +* + SRNAMT = 'SLATMS' + CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from SLATMS. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'SLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) + GO TO 160 + END IF +* +* For types 3-6, zero one or more rows and columns of the +* matrix to test that INFO is returned correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = ZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = ZERO + 50 CONTINUE + END IF + ELSE + IOFF = 0 + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = ZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = ZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* + DO 150 IFACT = 1, NFACT +* +* Do first for FACT = 'F', then for other values. +* + FACT = FACTS( IFACT ) +* +* Compute the condition number for comparison with +* the value returned by SSYSVX. +* + IF( ZEROT ) THEN + IF( IFACT.EQ.1 ) + $ GO TO 150 + RCONDC = ZERO +* + ELSE IF( IFACT.EQ.1 ) THEN +* +* Compute the 1-norm of A. +* + ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK ) +* +* Factor the matrix A. +* +c CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +c CALL SSYTRF( UPLO, N, AFAC, LDA, IWORK, WORK, +c $ LWORK, INFO ) +* +* Compute inv(A) and take its norm. +* +c CALL SLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) +c LWORK = (N+NB+1)*(NB+3) +c SRNAMT = 'DSYTRI2' +c CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK, +c $ LWORK, INFO ) +c AINVNM = SLANSY( '1', UPLO, N, AINV, LDA, RWORK ) +* +* Compute the 1-norm condition number of A. +* +c IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN +c RCONDC = ONE +c ELSE +c RCONDC = ( ONE / ANORM ) / AINVNM +c END IF + END IF +* +* Form an exact solution and set the right hand side. +* + SRNAMT = 'SLARHS' + CALL SLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU, + $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED, + $ INFO ) + XTYPE = 'C' +* +* --- Test SSYSV_AASEN --- +* + IF( IFACT.EQ.2 ) THEN + CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) + CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* +* Factor the matrix and solve the system using SSYSV_AASEN. +* + SRNAMT = 'SSYSV_AASEN' + CALL SSYSV_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK, + $ X, LDA, WORK, LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from SSYSV_AASEN . +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'SSYSV_AASEN ', INFO, K, + $ UPLO, N, N, -1, -1, NRHS, + $ IMAT, NFAIL, NERRS, NOUT ) + GO TO 120 + ELSE IF( INFO.NE.0 ) THEN + GO TO 120 + END IF +* +* Reconstruct matrix from factors and compute +* residual. +* + CALL SSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDA, + $ IWORK, AINV, LDA, RWORK, + $ RESULT( 1 ) ) +* +* Compute residual of the computed solution. +* + CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* Check solution from generated exact solution. +* + CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, + $ RESULT( 3 ) ) + NT = 3 +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALADHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )'SSYSV_AASEN ', + $ UPLO, N, IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT + 120 CONTINUE + END IF +* + 150 CONTINUE +* + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2, + $ ', test ', I2, ', ratio =', G12.5 ) + RETURN +* +* End of SDRVSY_AASEN +* + END diff --git a/TESTING/LIN/serrvx.f b/TESTING/LIN/serrvx.f index 7d999c4a..41334f5c 100644 --- a/TESTING/LIN/serrvx.f +++ b/TESTING/LIN/serrvx.f @@ -91,7 +91,7 @@ EXTERNAL CHKXER, SGBSV, SGBSVX, SGESV, SGESVX, SGTSV, $ SGTSVX, SPBSV, SPBSVX, SPOSV, SPOSVX, SPPSV, $ SPPSVX, SPTSV, SPTSVX, SSPSV, SSPSVX, SSYSV, - $ SSYSV_ROOK, SSYSVX + $ SSYSV_AASEN, SSYSV_ROOK, SSYSVX * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -626,6 +626,24 @@ CALL SSYSVX( 'N', 'U', 2, 0, A, 2, AF, 2, IP, B, 2, X, 2, $ RCOND, R1, R2, W, 3, IW, INFO ) CALL CHKXER( 'SSYSVX', INFOT, NOUT, LERR, OK ) +* + ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN +* +* SSYSV_AASEN +* + SRNAMT = 'SSYSV_AASEN' + INFOT = 1 + CALL SSYSV_AASEN( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'SSYSV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL SSYSV_AASEN( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'SSYSV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL SSYSV_AASEN( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'SSYSV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL SSYSV_AASEN( 'U', 2, 0, A, 2, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'SSYSV_AASEN', INFOT, NOUT, LERR, OK ) * ELSE IF( LSAMEN( 2, C2, 'SR' ) ) THEN * diff --git a/TESTING/LIN/ssyt01_aasen.f b/TESTING/LIN/ssyt01_aasen.f new file mode 100644 index 00000000..01163938 --- /dev/null +++ b/TESTING/LIN/ssyt01_aasen.f @@ -0,0 +1,262 @@ +*> \brief \b SSYT01_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE SSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, +* C, LDC, RWORK, RESID ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER LDA, LDAFAC, LDC, N +* REAL RESID +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* REAL A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), +* $ RWORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SSYT01_AASEN reconstructs a symmetric indefinite matrix A from its +*> block L*D*L' or U*D*U' factorization and computes the residual +*> norm( C - A ) / ( N * norm(A) * EPS ), +*> where C is the reconstructed matrix and EPS is the machine epsilon. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the upper or lower triangular part of the +*> symmetric matrix A is stored: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of rows and columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is REAL array, dimension (LDA,N) +*> The original symmetric matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N) +*> \endverbatim +*> +*> \param[in] AFAC +*> \verbatim +*> AFAC is REAL array, dimension (LDAFAC,N) +*> The factored form of the matrix A. AFAC contains the block +*> diagonal matrix D and the multipliers used to obtain the +*> factor L or U from the block L*D*L' or U*D*U' factorization +*> as computed by SSYTRF. +*> \endverbatim +*> +*> \param[in] LDAFAC +*> \verbatim +*> LDAFAC is INTEGER +*> The leading dimension of the array AFAC. LDAFAC >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> The pivot indices from SSYTRF. +*> \endverbatim +*> +*> \param[out] C +*> \verbatim +*> C is REAL array, dimension (LDC,N) +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (N) +*> \endverbatim +*> +*> \param[out] RESID +*> \verbatim +*> RESID is REAL +*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) +*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* +*> \ingroup real_lin +* +* ===================================================================== + SUBROUTINE SSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, + $ LDC, RWORK, RESID ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDAFAC, LDC, N + REAL RESID +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + REAL A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), + $ RWORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J + REAL ANORM, EPS +* .. +* .. External Functions .. + LOGICAL LSAME + REAL SLAMCH, SLANSY + EXTERNAL LSAME, SLAMCH, SLANSY +* .. +* .. External Subroutines .. + EXTERNAL SLASET, SLAVSY +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE +* .. +* .. Executable Statements .. +* +* Quick exit if N = 0. +* + IF( N.LE.0 ) THEN + RESID = ZERO + RETURN + END IF +* +* Determine EPS and the norm of A. +* + EPS = SLAMCH( 'Epsilon' ) + ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK ) +* +* Initialize C to the tridiagonal matrix T. +* + CALL SLASET( 'Full', N, N, ZERO, ZERO, C, LDC ) + CALL SLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 ) + IF( N.GT.1 ) THEN + IF( LSAME( UPLO, 'U' ) ) THEN + CALL SLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL SLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + ELSE + CALL SLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL SLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + ENDIF + ENDIF +* +* Call STRMM to form the product U' * D (or L * D ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL STRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N, + $ ONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) + ELSE + CALL STRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ ONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF +* +* Call STRMM again to multiply by U (or L ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL STRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ ONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL STRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1, + $ ONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) + END IF +* +* Apply symmetric pivots +* + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL SSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC ) + END DO + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL SSWAP( N, C( 1, J ), 1, C( 1, I ), 1 ) + END DO +* +* +* Compute the difference C - A . +* + IF( LSAME( UPLO, 'U' ) ) THEN + DO J = 1, N + DO I = 1, J + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + ELSE + DO J = 1, N + DO I = J, N + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + END IF +* +* Compute norm( C - A ) / ( N * norm(A) * EPS ) +* + RESID = SLANSY( '1', UPLO, N, C, LDC, RWORK ) +* + IF( ANORM.LE.ZERO ) THEN + IF( RESID.NE.ZERO ) + $ RESID = ONE / EPS + ELSE + RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS + END IF +* + RETURN +* +* End of SSYT01_AASEN +* + END diff --git a/TESTING/LIN/zchkaa.f b/TESTING/LIN/zchkaa.f index 9212afc9..90b98a2e 100644 --- a/TESTING/LIN/zchkaa.f +++ b/TESTING/LIN/zchkaa.f @@ -50,6 +50,7 @@ *> ZPB 8 List types on next line if 0 < NTYPES < 8 *> ZPT 12 List types on next line if 0 < NTYPES < 12 *> ZHE 10 List types on next line if 0 < NTYPES < 10 +*> ZHA 10 List types on next line if 0 < NTYPES < 10 *> ZHR 10 List types on next line if 0 < NTYPES < 10 *> ZHP 10 List types on next line if 0 < NTYPES < 10 *> ZSY 11 List types on next line if 0 < NTYPES < 11 @@ -102,17 +103,17 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date November 2015 +*> \date November 2016 * *> \ingroup complex16_lin * * ===================================================================== PROGRAM ZCHKAA * -* -- LAPACK test routine (version 3.6.0) -- +* -- LAPACK test routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2015 +* November 2016 * * ===================================================================== * @@ -638,6 +639,33 @@ ELSE WRITE( NOUT, FMT = 9988 )PATH END IF +* + ELSE IF( LSAMEN( 2, C2, 'HA' ) ) THEN +* +* HA: Hermitian indefinite matrices, +* with partial (Aasen's) pivoting algorithm +* + NTYPES = 10 + CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) +* + IF( TSTCHK ) THEN + CALL ZCHKHE_AASEN( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS, + $ NSVAL, THRESH, TSTERR, LDA, + $ A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), + $ WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9989 )PATH + END IF +* + IF( TSTDRV ) THEN + CALL ZDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), + $ WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9988 )PATH + END IF * ELSE IF( LSAMEN( 2, C2, 'HR' ) ) THEN * diff --git a/TESTING/LIN/zchkhe_aasen.f b/TESTING/LIN/zchkhe_aasen.f new file mode 100644 index 00000000..1ebe5ef4 --- /dev/null +++ b/TESTING/LIN/zchkhe_aasen.f @@ -0,0 +1,579 @@ +*> \brief \b ZCHKHE_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE ZCHKHE_AASEN( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, +* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, +* XACT, WORK, RWORK, IWORK, NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NN, NNB, NNS, NOUT +* DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) +* DOUBLE PRECISION RWORK( * ) +* COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), +* $ WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZCHKHE_AASEN tests ZHETRF_AASEN, -TRS_AASEN. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NNB +*> \verbatim +*> NNB is INTEGER +*> The number of values of NB contained in the vector NBVAL. +*> \endverbatim +*> +*> \param[in] NBVAL +*> \verbatim +*> NBVAL is INTEGER array, dimension (NBVAL) +*> The values of the blocksize NB. +*> \endverbatim +*> +*> \param[in] NNS +*> \verbatim +*> NNS is INTEGER +*> The number of values of NRHS contained in the vector NSVAL. +*> \endverbatim +*> +*> \param[in] NSVAL +*> \verbatim +*> NSVAL is INTEGER array, dimension (NNS) +*> The values of the number of right hand sides NRHS. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is DOUBLE PRECISION +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (NMAX*NSMAX) +*> where NSMAX is the largest entry in NSVAL. +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is COMPLEX*16 array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is COMPLEX*16 array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (NMAX*max(3,NSMAX)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (max(NMAX,2*NSMAX)) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* +*> \ingroup complex16_lin +* +* ===================================================================== + SUBROUTINE ZCHKHE_AASEN( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, + $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, + $ X, XACT, WORK, RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NMAX, NN, NNB, NNS, NOUT + DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), + $ WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) + COMPLEX*16 CZERO + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) + INTEGER NTYPES + PARAMETER ( NTYPES = 10 ) + INTEGER NTESTS + PARAMETER ( NTESTS = 9 ) +* .. +* .. Local Scalars .. + LOGICAL TRFCON, ZEROT + CHARACTER DIST, TYPE, UPLO, XTYPE + CHARACTER*3 PATH, MATPATH + INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, + $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, + $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT + DOUBLE PRECISION ANORM, CNDNUM, RCOND, RCONDC +* .. +* .. Local Arrays .. + CHARACTER UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + DOUBLE PRECISION RESULT( NTESTS ) +* .. +* .. External Functions .. + DOUBLE PRECISION DGET06, ZLANHE + EXTERNAL DGET06, ZLANHE +* .. +* .. External Subroutines .. + EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRHE, ZGET04, + $ ZHECON, ZHERFS, ZHET01, ZHETRF_AASEN, ZHETRI2, + $ ZHETRS_AASEN, ZLACPY, ZLAIPD, ZLARHS, ZLATB4, + $ ZLATMS, ZPOT02, ZPOT03, ZPOT05 +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Zomplex precision' + PATH( 2: 3 ) = 'HA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Zomplex precision' + MATPATH( 2: 3 ) = 'HE' + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE +* +* Test the error exits +* + IF( TSTERR ) + $ CALL ZERRHE( PATH, NOUT ) + INFOT = 0 +* +* Set the minimum block size for which the block routine should +* be used, which will be later returned by ILAENV +* + CALL XLAENV( 2, 2 ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + IF( N .GT. NMAX ) THEN + NFAIL = NFAIL + 1 + WRITE(NOUT, 9995) 'M ', N, NMAX + GO TO 180 + END IF + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + IZERO = 0 + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Set up parameters with ZLATB4 for the matrix generator +* based on the type of matrix to be generated. +* + CALL ZLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, + $ ANORM, MODE, CNDNUM, DIST ) +* +* Generate a matrix with ZLATMS. +* + SRNAMT = 'ZLATMS' + CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from ZLATMS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) +* +* Skip all tests for this generated matrix +* + GO TO 160 + END IF +* +* For types 3-6, zero one or more rows and columns of +* the matrix to test that INFO is returned correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = CZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = CZERO + 50 CONTINUE + END IF + ELSE + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + IOFF = 0 + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = CZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + IOFF = 0 + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = CZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* +* End generate test matrix A. +* +* +* Set the imaginary part of the diagonals. +* + CALL ZLAIPD( N, A, LDA+1, 0 ) +* +* Do for each value of NB in NBVAL +* + DO 150 INB = 1, NNB +* +* Set the optimal blocksize, which will be later +* returned by ILAENV. +* + NB = NBVAL( INB ) + CALL XLAENV( 1, NB ) +* +* Copy the test matrix A into matrix AFAC which +* will be factorized in place. This is needed to +* preserve the test matrix A for subsequent tests. +* + CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +* +* Compute the L*D*L**T or U*D*U**T factorization of the +* matrix. IWORK stores details of the interchanges and +* the block structure of D. AINV is a work array for +* block factorization, LWORK is the length of AINV. +* + LWORK = ( NB+1 )*LDA + SRNAMT = 'ZHETRF_AASEN' + CALL ZHETRF_AASEN( UPLO, N, AFAC, LDA, IWORK, AINV, + $ LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from ZHETRF and handle error. +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'ZHETRF_AASEN', INFO, K, UPLO, + $ N, N, -1, -1, NB, IMAT, NFAIL, NERRS, + $ NOUT ) + END IF +* +* Set the condition estimate flag if the INFO is not 0. +* + IF( INFO.NE.0 ) THEN + TRFCON = .TRUE. + ELSE + TRFCON = .FALSE. + END IF +* +*+ TEST 1 +* Reconstruct matrix from factors and compute residual. +* + CALL ZHET01_AASEN( UPLO, N, A, LDA, AFAC, LDA, IWORK, + $ AINV, LDA, RWORK, RESULT( 1 ) ) + NT = 1 +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, + $ RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT +* +* Do only the condition estimate if INFO is not 0. +* + IF( TRFCON ) THEN + RCONDC = ZERO + GO TO 140 + END IF +* +* Do for each value of NRHS in NSVAL. +* + DO 130 IRHS = 1, NNS + NRHS = NSVAL( IRHS ) +* +*+ TEST 3 (Using TRS) +* Solve and compute residual for A * X = B. +* +* Choose a set of NRHS random solution vectors +* stored in XACT and set up the right hand side B +* + SRNAMT = 'ZLARHS' + CALL ZLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, + $ KL, KU, NRHS, A, LDA, XACT, LDA, + $ B, LDA, ISEED, INFO ) + CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* + SRNAMT = 'ZHETRS_AASEN' + LWORK = 3*N-2 + CALL ZHETRS_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK, + $ X, LDA, WORK, LWORK, INFO ) +* +* Check error code from ZHETRS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'ZHETRS', INFO, 0, UPLO, N, + $ N, -1, -1, NRHS, IMAT, NFAIL, + $ NERRS, NOUT ) + END IF +* + CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) +* +* Compute the residual for the solution +* + CALL ZPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* Print information about the tests that did not pass +* the threshold. +* + DO 120 K = 2, 2 + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, + $ IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 120 CONTINUE + NRUN = NRUN + 1 +* +* End do for each value of NRHS in NSVAL. +* + 130 CONTINUE + 140 CONTINUE + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ', + $ I2, ', test ', I2, ', ratio =', G12.5 ) + 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', + $ I2, ', test(', I2, ') =', G12.5 ) +c 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2, +c $ ', test(', I2, ') =', G12.5 ) + 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=', + $ I6 ) + RETURN +* +* End of ZCHKHE_AASEN +* + END diff --git a/TESTING/LIN/zdrvhe_aasen.f b/TESTING/LIN/zdrvhe_aasen.f new file mode 100644 index 00000000..0170873e --- /dev/null +++ b/TESTING/LIN/zdrvhe_aasen.f @@ -0,0 +1,525 @@ +*> \brief \b ZDRVHE_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE ZDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NOUT, NRHS +* DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NVAL( * ) +* DOUBLE PRECISION RWORK( * ) +* COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), +* $ WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZDRVHE_AASEN tests the driver routine ZHESV_AASEN. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand side vectors to be generated for +*> each linear system. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is DOUBLE PRECISION +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is COMPLEX*16 array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is COMPLEX*16 array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (NMAX*max(2,NRHS)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (NMAX+2*NRHS) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup complex16_lin +* +* ===================================================================== + SUBROUTINE ZDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, + $ RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NMAX, NN, NOUT, NRHS + DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NVAL( * ) + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), + $ WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + INTEGER NTYPES, NTESTS + PARAMETER ( NTYPES = 10, NTESTS = 3 ) + INTEGER NFACT + PARAMETER ( NFACT = 2 ) +* .. +* .. Local Scalars .. + LOGICAL ZEROT + CHARACTER DIST, FACT, TYPE, UPLO, XTYPE + CHARACTER*3 MATPATH, PATH + INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, + $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N, + $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT + DOUBLE PRECISION AINVNM, ANORM, CNDNUM, RCOND, RCONDC +* .. +* .. Local Arrays .. + CHARACTER FACTS( NFACT ), UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + DOUBLE PRECISION RESULT( NTESTS ) +* .. +* .. External Functions .. + DOUBLE PRECISION DGET06, ZLANHE + EXTERNAL DGET06, ZLANHE +* .. +* .. External Subroutines .. + EXTERNAL ALADHD, ALAERH, ALASVM, XLAENV, ZERRVX, ZGET04, + $ ZHESV_AASEN, ZHET01_AASEN, ZHETRF_AASEN, + $ ZHETRI2, ZLACPY, ZLAIPD, ZLARHS, ZLATB4, ZLATMS, + $ ZPOT02 +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Intrinsic Functions .. + INTRINSIC DCMPLX, MAX, MIN +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Zomplex precision' + PATH( 2: 3 ) = 'HA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Zomplex precision' + MATPATH( 2: 3 ) = 'HE' +* + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE + LWORK = MAX( 2*NMAX, NMAX*NRHS ) +* +* Test the error exits +* + IF( TSTERR ) + $ CALL ZERRVX( PATH, NOUT ) + INFOT = 0 +* +* Set the block size and minimum block size for testing. +* + NB = 1 + NBMIN = 2 + CALL XLAENV( 1, NB ) + CALL XLAENV( 2, NBMIN ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Begin generate the test matrix A. +* +* Set up parameters with ZLATB4 and generate a test matrix +* with ZLATMS. +* + CALL ZLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM, + $ MODE, CNDNUM, DIST ) +* + SRNAMT = 'ZLATMS' + CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from ZLATMS. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) + GO TO 160 + END IF +* +* For types 3-6, zero one or more rows and columns of the +* matrix to test that INFO is returned correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = ZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = ZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = ZERO + 50 CONTINUE + END IF + ELSE + IOFF = 0 + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = ZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = ZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* +* Set the imaginary part of the diagonals. +* + CALL ZLAIPD( N, A, LDA+1, 0 ) +* + DO 150 IFACT = 1, NFACT +* +* Do first for FACT = 'F', then for other values. +* + FACT = FACTS( IFACT ) +* +* Compute the condition number for comparison with +* the value returned by ZHESVX. +* + IF( ZEROT ) THEN + IF( IFACT.EQ.1 ) + $ GO TO 150 + RCONDC = ZERO +* + ELSE IF( IFACT.EQ.1 ) THEN +* +* Compute the 1-norm of A. +* + ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK ) +* +* Factor the matrix A. +* +c CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +c CALL ZHETRF( UPLO, N, AFAC, LDA, IWORK, WORK, +c $ LWORK, INFO ) +* +* Compute inv(A) and take its norm. +* +c CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) +c LWORK = (N+NB+1)*(NB+3) +c CALL ZHETRI2( UPLO, N, AINV, LDA, IWORK, WORK, +c $ LWORK, INFO ) +c AINVNM = ZLANHE( '1', UPLO, N, AINV, LDA, RWORK ) +* +* Compute the 1-norm condition number of A. +* +c IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN +c RCONDC = ONE +c ELSE +c RCONDC = ( ONE / ANORM ) / AINVNM +c END IF + END IF +* +* Form an exact solution and set the right hand side. +* + SRNAMT = 'ZLARHS' + CALL ZLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU, + $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED, + $ INFO ) + XTYPE = 'C' +* +* --- Test ZHESV_AASEN --- +* + IF( IFACT.EQ.2 ) THEN + CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) + CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* +* Factor the matrix and solve the system using ZHESV. +* + SRNAMT = 'ZHESV_AASEN ' + CALL ZHESV_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK, + $ X, LDA, WORK, LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from ZHESV . +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'ZHESV_AASEN', INFO, K, UPLO, N, + $ N, -1, -1, NRHS, IMAT, NFAIL, + $ NERRS, NOUT ) + GO TO 120 + ELSE IF( INFO.NE.0 ) THEN + GO TO 120 + END IF +* +* Reconstruct matrix from factors and compute +* residual. +* + CALL ZHET01_AASEN( UPLO, N, A, LDA, AFAC, LDA, + $ IWORK, AINV, LDA, RWORK, + $ RESULT( 1 ) ) +* +* Compute residual of the computed solution. +* + CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL ZPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* Check solution from generated exact solution. +* + CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, + $ RESULT( 3 ) ) + NT = 3 +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALADHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )'ZHESV_AASEN', UPLO, N, + $ IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT + 120 CONTINUE + END IF +* + 150 CONTINUE +* + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2, + $ ', test ', I2, ', ratio =', G12.5 ) + RETURN +* +* End of ZDRVHE_AASEN +* + END diff --git a/TESTING/LIN/zerrvx.f b/TESTING/LIN/zerrvx.f index 11a35aa4..1bb2d222 100644 --- a/TESTING/LIN/zerrvx.f +++ b/TESTING/LIN/zerrvx.f @@ -93,7 +93,7 @@ $ ZGTSVX, ZHESV, ZHESV_ROOK, ZHESVX, ZHPSV, $ ZHPSVX, ZPBSV, ZPBSVX, ZPOSV, ZPOSVX, ZPPSV, $ ZPPSVX, ZPTSV, ZPTSVX, ZSPSV, ZSPSVX, ZSYSV, - $ ZSYSV_ROOK, ZSYSVX + $ ZSYSV_AASEN, ZSYSV_ROOK, ZSYSVX * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -634,6 +634,25 @@ $ RCOND, R1, R2, W, 3, RW, INFO ) CALL CHKXER( 'ZHESVX', INFOT, NOUT, LERR, OK ) * + ELSE IF( LSAMEN( 2, C2, 'HA' ) ) THEN +* +* ZHESV_AASEN +* + SRNAMT = 'ZHESV_AASEN' + INFOT = 1 + CALL ZHESV_AASEN( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'ZHESV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL ZHESV_AASEN( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'ZHESV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL ZHESV_AASEN( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'ZHESV_AASEN', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL ZHESV_AASEN( 'U', 2, 0, A, 2, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'ZHESV_AASEN', INFOT, NOUT, LERR, OK ) +* + ELSE IF( LSAMEN( 2, C2, 'HR' ) ) THEN * * ZHESV_ROOK diff --git a/TESTING/LIN/zhet01_aasen.f b/TESTING/LIN/zhet01_aasen.f new file mode 100644 index 00000000..89b87990 --- /dev/null +++ b/TESTING/LIN/zhet01_aasen.f @@ -0,0 +1,267 @@ +*> \brief \b ZHET01_AASEN +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE ZHET01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, +* C, LDC, RWORK, RESID ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER LDA, LDAFAC, LDC, N +* COMPLEX*16 RESID +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), +* $ RWORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZHET01_AASEN reconstructs a hermitian indefinite matrix A from its +*> block L*D*L' or U*D*U' factorization and computes the residual +*> norm( C - A ) / ( N * norm(A) * EPS ), +*> where C is the reconstructed matrix and EPS is the machine epsilon. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the upper or lower triangular part of the +*> hermitian matrix A is stored: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of rows and columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> The original hermitian matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N) +*> \endverbatim +*> +*> \param[in] AFAC +*> \verbatim +*> AFAC is COMPLEX*16 array, dimension (LDAFAC,N) +*> The factored form of the matrix A. AFAC contains the block +*> diagonal matrix D and the multipliers used to obtain the +*> factor L or U from the block L*D*L' or U*D*U' factorization +*> as computed by ZHETRF. +*> \endverbatim +*> +*> \param[in] LDAFAC +*> \verbatim +*> LDAFAC is INTEGER +*> The leading dimension of the array AFAC. LDAFAC >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> The pivot indices from ZHETRF. +*> \endverbatim +*> +*> \param[out] C +*> \verbatim +*> C is COMPLEX*16 array, dimension (LDC,N) +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is COMPLEX*16 array, dimension (N) +*> \endverbatim +*> +*> \param[out] RESID +*> \verbatim +*> RESID is COMPLEX*16 +*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) +*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* +*> \ingroup complex16_lin +* +* ===================================================================== + SUBROUTINE ZHET01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, + $ LDC, RWORK, RESID ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDAFAC, LDC, N + DOUBLE PRECISION RESID +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), + $ RWORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), + $ CONE = ( 1.0D+0, 0.0D+0 ) ) + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J + DOUBLE PRECISION ANORM, EPS +* .. +* .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, ZLANHE + EXTERNAL LSAME, DLAMCH, ZLANHE +* .. +* .. External Subroutines .. + EXTERNAL ZLASET, ZLAVHE +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE +* .. +* .. Executable Statements .. +* +* Quick exit if N = 0. +* + IF( N.LE.0 ) THEN + RESID = ZERO + RETURN + END IF +* +* Determine EPS and the norm of A. +* + EPS = DLAMCH( 'Epsilon' ) + ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK ) +* +* Initialize C to the tridiagonal matrix T. +* + CALL ZLASET( 'Full', N, N, CZERO, CZERO, C, LDC ) + CALL ZLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 ) + IF( N.GT.1 ) THEN + IF( LSAME( UPLO, 'U' ) ) THEN + CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + CALL ZLACGV( N-1, C( 2, 1 ), LDC+1 ) + ELSE + CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + CALL ZLACGV( N-1, C( 1, 2 ), LDC+1 ) + ENDIF + ENDIF +* +* Call ZTRMM to form the product U' * D (or L * D ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose', 'Unit', N-1, + $ N, CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) + ELSE + CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF +* +* Call ZTRMM again to multiply by U (or L ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose', 'Unit', N, + $ N-1, CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) + END IF +* +* Apply hermitian pivots +* + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL ZSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC ) + END DO + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL ZSWAP( N, C( 1, J ), 1, C( 1, I ), 1 ) + END DO +* +* +* Compute the difference C - A . +* + IF( LSAME( UPLO, 'U' ) ) THEN + DO J = 1, N + DO I = 1, J + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + ELSE + DO J = 1, N + DO I = J, N + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + END IF +* +* Compute norm( C - A ) / ( N * norm(A) * EPS ) +* + RESID = ZLANHE( '1', UPLO, N, C, LDC, RWORK ) +* + IF( ANORM.LE.ZERO ) THEN + IF( RESID.NE.ZERO ) + $ RESID = ONE / EPS + ELSE + RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS + END IF +* + RETURN +* +* End of ZHET01_AASEN +* + END -- cgit v1.2.3