OSchip
/
OpenBLAS
Not watched
1
0
Fork 0
Code
Releases 66 Wiki evaluate Activity
Issues 0 Pull Requests 0 Datasets Model Cloudbrain HPC
Browse Source

Merge pull request #4356 from martin-frbg/lapack736-2 

Add LAPACK tests for the Dynamic Mode Decomposition functions from Reference-LAPACK PR 736
tags/v0.3.26
Martin Kroeker GitHub 1 year ago
parent
e3039fa7f6 226a14c549
commit
5915a69734
No known key found for this signature in database GPG Key ID: 4AEE18F83AFDEB23
11 changed files with 3238 additions and 55 deletions
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. +41
    -9
      lapack-netlib/TESTING/EIG/Makefile
  2. +721
    -0
      lapack-netlib/TESTING/EIG/cchkdmd.f90
  3. +813
    -0
      lapack-netlib/TESTING/EIG/dchkdmd.f90
  4. +792
    -0
      lapack-netlib/TESTING/EIG/schkdmd.f90
  5. +745
    -0
      lapack-netlib/TESTING/EIG/zchkdmd.f90
  6. +39
    -4
      lapack-netlib/TESTING/Makefile
  7. +11
    -0
      lapack-netlib/TESTING/cdmd.in
  8. +11
    -0
      lapack-netlib/TESTING/ddmd.in
  9. +11
    -0
      lapack-netlib/TESTING/sdmd.in
  10. +11
    -0
      lapack-netlib/TESTING/zdmd.in
  11. +43
    -42
      lapack-netlib/lapack_testing.py

+ 41
- 9
lapack-netlib/TESTING/EIG/Makefile View File

@@ -64,6 +64,8 @@ SEIGTST = schkee.o \
sort03.o ssbt21.o ssgt01.o sslect.o sspt21.o sstt21.o \ sort03.o ssbt21.o ssgt01.o sslect.o sspt21.o sstt21.o \
sstt22.o ssyl01.o ssyt21.o ssyt22.o sstt22.o ssyl01.o ssyt21.o ssyt22.o


SDMDEIGTST = schkdmd.o

CEIGTST = cchkee.o \ CEIGTST = cchkee.o \
cbdt01.o cbdt02.o cbdt03.o cbdt05.o \ cbdt01.o cbdt02.o cbdt03.o cbdt05.o \
cchkbb.o cchkbd.o cchkbk.o cchkbl.o cchkec.o \ cchkbb.o cchkbd.o cchkbk.o cchkbl.o cchkec.o \
@@ -81,6 +83,8 @@ CEIGTST = cchkee.o \
csgt01.o cslect.o csyl01.o\ csgt01.o cslect.o csyl01.o\
cstt21.o cstt22.o cunt01.o cunt03.o cstt21.o cstt22.o cunt01.o cunt03.o


CDMDEIGTST = cchkdmd.o

DZIGTST = dlafts.o dlahd2.o dlasum.o dlatb9.o dstech.o dstect.o \ DZIGTST = dlafts.o dlahd2.o dlasum.o dlatb9.o dstech.o dstect.o \
dsvdch.o dsvdct.o dsxt1.o dsvdch.o dsvdct.o dsxt1.o


@@ -101,6 +105,8 @@ DEIGTST = dchkee.o \
dort03.o dsbt21.o dsgt01.o dslect.o dspt21.o dstt21.o \ dort03.o dsbt21.o dsgt01.o dslect.o dspt21.o dstt21.o \
dstt22.o dsyl01.o dsyt21.o dsyt22.o dstt22.o dsyl01.o dsyt21.o dsyt22.o


DDMDEIGTST = dchkdmd.o

ZEIGTST = zchkee.o \ ZEIGTST = zchkee.o \
zbdt01.o zbdt02.o zbdt03.o zbdt05.o \ zbdt01.o zbdt02.o zbdt03.o zbdt05.o \
zchkbb.o zchkbd.o zchkbk.o zchkbl.o zchkec.o \ zchkbb.o zchkbd.o zchkbk.o zchkbl.o zchkec.o \
@@ -118,27 +124,45 @@ ZEIGTST = zchkee.o \
zsgt01.o zslect.o zsyl01.o\ zsgt01.o zslect.o zsyl01.o\
zstt21.o zstt22.o zunt01.o zunt03.o zstt21.o zstt22.o zunt01.o zunt03.o


ZDMDEIGTST = zchkdmd.o

.PHONY: all .PHONY: all
all: single complex double complex16 all: single complex double complex16


.PHONY: single complex double complex16 .PHONY: single complex double complex16
single: xeigtsts
complex: xeigtstc
double: xeigtstd
complex16: xeigtstz
single: xeigtsts xdmdeigtsts
complex: xeigtstc xdmdeigtstc
double: xeigtstd xdmdeigtstd
complex16: xeigtstz xdmdeigtstz

xdmdeigtsts: $(SDMDEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB)
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^

xdmdeigtstc: $(CDMDEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB)
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^

xdmdeigtstd: $(DDMDEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB)
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^

xdmdeigtstz: $(ZDMDEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB)
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^


xeigtsts: $(SEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) xeigtsts: $(SEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^


xeigtstc: $(CEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) xeigtstc: $(CEIGTST) $(SCIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^


xeigtstd: $(DEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) xeigtstd: $(DEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^


xeigtstz: $(ZEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB) xeigtstz: $(ZEIGTST) $(DZIGTST) $(AEIGTST) $(TMGLIB) ../$(LAPACKLIB) $(BLASLIB)
$(LOADER) $(FFLAGS) $(LDFLAGS) -o $@ $^
$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^


$(SDMDEIGTST): $(FRC)
$(CDMDEIGTST): $(FRC)
$(DDMDEIGTST): $(FRC)
$(ZDMDEIGTST): $(FRC)
$(AEIGTST): $(FRC) $(AEIGTST): $(FRC)
$(SCIGTST): $(FRC) $(SCIGTST): $(FRC)
$(DZIGTST): $(FRC) $(DZIGTST): $(FRC)
@@ -155,7 +179,7 @@ clean: cleanobj cleanexe
cleanobj: cleanobj:
rm -f *.o rm -f *.o
cleanexe: cleanexe:
rm -f xeigtst*
rm -f xeigtst* xdmdeigtst*


schkee.o: schkee.F schkee.o: schkee.F
$(FC) $(FFLAGS_DRV) -c -o $@ $< $(FC) $(FFLAGS_DRV) -c -o $@ $<
@@ -165,3 +189,11 @@ cchkee.o: cchkee.F
$(FC) $(FFLAGS_DRV) -c -o $@ $< $(FC) $(FFLAGS_DRV) -c -o $@ $<
zchkee.o: zchkee.F zchkee.o: zchkee.F
$(FC) $(FFLAGS_DRV) -c -o $@ $< $(FC) $(FFLAGS_DRV) -c -o $@ $<
schkdmd.o: schkdmd.f90
$(FC) $(FFLAGS_DRV) -c -o $@ $<
cchkdmd.o: cchkdmd.f90
$(FC) $(FFLAGS_DRV) -c -o $@ $<
dchkdmd.o: dchkdmd.f90
$(FC) $(FFLAGS_DRV) -c -o $@ $<
zchkdmd.o: zchkdmd.f90
$(FC) $(FFLAGS_DRV) -c -o $@ $<

+ 721
- 0
lapack-netlib/TESTING/EIG/cchkdmd.f90 View File

@@ -0,0 +1,721 @@
! This is a test program for checking the implementations of
! the implementations of the following subroutines
!
! CGEDMD, for computation of the
! Dynamic Mode Decomposition (DMD)
! CGEDMDQ, for computation of a
! QR factorization based compressed DMD
!
! Developed and supported by:
! ===========================
! Developed and coded by Zlatko Drmac, Faculty of Science,
! University of Zagreb; drmac@math.hr
! In cooperation with
! AIMdyn Inc., Santa Barbara, CA.
! ========================================================
! How to run the code (compiler, link info)
! ========================================================
! Compile as FORTRAN 90 (or later) and link with BLAS and
! LAPACK libraries.
! NOTE: The code is developed and tested on top of the
! Intel MKL library (versions 2022.0.3 and 2022.2.0),
! using the Intel Fortran compiler.
!
! For developers of the C++ implementation
! ========================================================
! See the LAPACK++ and Template Numerical Toolkit (TNT)
!
! Note on a development of the GPU HP implementation
! ========================================================
! Work in progress. See CUDA, MAGMA, SLATE.
! NOTE: The four SVD subroutines used in this code are
! included as a part of R&D and for the completeness.
! This was also an opportunity to test those SVD codes.
! If the scaling option is used all four are essentially
! equally good. For implementations on HP platforms,
! one can use whichever SVD is available.
!............................................................
!............................................................
!............................................................
!
PROGRAM DMD_TEST
use iso_fortran_env
IMPLICIT NONE
integer, parameter :: WP = real32
!............................................................
REAL(KIND=WP), PARAMETER :: ONE = 1.0_WP
REAL(KIND=WP), PARAMETER :: ZERO = 0.0_WP
COMPLEX(KIND=WP), PARAMETER :: CONE = ( 1.0_WP, 0.0_WP )
COMPLEX(KIND=WP), PARAMETER :: CZERO = ( 0.0_WP, 0.0_WP )
!............................................................
REAL(KIND=WP), ALLOCATABLE, DIMENSION(:) :: RES, &
RES1, RESEX, SINGVX, SINGVQX, WORK
INTEGER , ALLOCATABLE, DIMENSION(:) :: IWORK
REAL(KIND=WP) :: WDUMMY(2)
INTEGER :: IDUMMY(4), ISEED(4)
REAL(KIND=WP) :: ANORM, COND, CONDL, CONDR, EPS, &
TOL, TOL2, SVDIFF, TMP, TMP_AU, &
TMP_FQR, TMP_REZ, TMP_REZQ, TMP_XW, &
TMP_EX
!............................................................
COMPLEX(KIND=WP) :: CMAX
INTEGER :: LCWORK
COMPLEX(KIND=WP), ALLOCATABLE, DIMENSION(:,:) :: A, AC, &
AU, F, F0, F1, S, W, &
X, X0, Y, Y0, Y1, Z, Z1
COMPLEX(KIND=WP), ALLOCATABLE, DIMENSION(:) :: CDA, CDR, &
CDL, CEIGS, CEIGSA, CWORK
COMPLEX(KIND=WP) :: CDUMMY(22), CDUM2X2(2,2)
!............................................................
INTEGER :: K, KQ, LDF, LDS, LDA, LDAU, LDW, LDX, LDY, &
LDZ, LIWORK, LWORK, M, N, LLOOP, NRNK
INTEGER :: i, iJOBREF, iJOBZ, iSCALE, INFO, j, &
NFAIL, NFAIL_AU, NFAIL_F_QR, NFAIL_REZ, &
NFAIL_REZQ, NFAIL_SVDIFF, NFAIL_TOTAL, NFAILQ_TOTAL, &
NFAIL_Z_XV, MODE, MODEL, MODER, WHTSVD
INTEGER :: iNRNK, iWHTSVD, K_traj, LWMINOPT
CHARACTER :: GRADE, JOBREF, JOBZ, PIVTNG, RSIGN, &
SCALE, RESIDS, WANTQ, WANTR
LOGICAL :: TEST_QRDMD
!..... external subroutines (BLAS and LAPACK)
EXTERNAL CAXPY, CGEEV, CGEMM, CGEMV, CLASCL
!.....external subroutines DMD package
! subroutines under test
EXTERNAL CGEDMD, CGEDMDQ
!..... external functions (BLAS and LAPACK)
EXTERNAL SCNRM2, SLAMCH
REAL(KIND=WP) :: SCNRM2, SLAMCH
EXTERNAL CLANGE
REAL(KIND=WP) :: CLANGE
EXTERNAL ICAMAX
INTEGER ICAMAX
EXTERNAL LSAME
LOGICAL LSAME
INTRINSIC ABS, INT, MIN, MAX, SIGN
!............................................................
WRITE(*,*) 'COMPLEX CODE TESTING'
! The test is always in pairs : ( CGEDMD and CGEDMDQ)
! because the test includes comparing the results (in pairs).
!.....................................................................................
! This code by default performs tests on CGEDMDQ
! Since the QR factorizations based algorithm is designed for
! single trajectory data, only single trajectory tests will
! be performed with xGEDMDQ.
WANTQ = 'Q'
WANTR = 'R'
!.................................................................................
EPS = SLAMCH( 'P' ) ! machine precision WP
! Global counters of failures of some particular tests
NFAIL = 0
NFAIL_REZ = 0
NFAIL_REZQ = 0
NFAIL_Z_XV = 0
NFAIL_F_QR = 0
NFAIL_AU = 0
NFAIL_SVDIFF = 0
NFAIL_TOTAL = 0
NFAILQ_TOTAL = 0
DO LLOOP = 1, 4
WRITE(*,*) 'L Loop Index = ', LLOOP
! Set the dimensions of the problem ...
READ(*,*) M
WRITE(*,*) 'M = ', M
! ... and the number of snapshots.
READ(*,*) N
WRITE(*,*) 'N = ', N
! Test the dimensions
IF ( ( MIN(M,N) == 0 ) .OR. ( M < N ) ) THEN
WRITE(*,*) 'Bad dimensions. Required: M >= N > 0.'
STOP
END IF
!.............
! The seed inside the LLOOP so that each pass can be reproduced easily.
ISEED(1) = 4
ISEED(2) = 3
ISEED(3) = 2
ISEED(4) = 1
LDA = M
LDF = M
LDX = M
LDY = M
LDW = N
LDZ = M
LDAU = M
LDS = N
TMP_XW = ZERO
TMP_AU = ZERO
TMP_REZ = ZERO
TMP_REZQ = ZERO
SVDIFF = ZERO
TMP_EX = ZERO
ALLOCATE( A(LDA,M) )
ALLOCATE( AC(LDA,M) )
ALLOCATE( F(LDF,N+1) )
ALLOCATE( F0(LDF,N+1) )
ALLOCATE( F1(LDF,N+1) )
ALLOCATE( X(LDX,N) )
ALLOCATE( X0(LDX,N) )
ALLOCATE( Y(LDY,N+1) )
ALLOCATE( Y0(LDY,N+1) )
ALLOCATE( Y1(LDY,N+1) )
ALLOCATE( AU(LDAU,N) )
ALLOCATE( W(LDW,N) )
ALLOCATE( S(LDS,N) )
ALLOCATE( Z(LDZ,N) )
ALLOCATE( Z1(LDZ,N) )
ALLOCATE( RES(N) )
ALLOCATE( RES1(N) )
ALLOCATE( RESEX(N) )
ALLOCATE( CEIGS(N) )
ALLOCATE( SINGVX(N) )
ALLOCATE( SINGVQX(N) )
TOL = 10*M*EPS
TOL2 = 10*M*N*EPS
!.............
DO K_traj = 1, 2
! Number of intial conditions in the simulation/trajectories (1 or 2)
COND = 1.0D4
CMAX = (1.0D1,1.0D1)
RSIGN = 'F'
GRADE = 'N'
MODEL = 6
CONDL = 1.0D1
MODER = 6
CONDR = 1.0D1
PIVTNG = 'N'
! Loop over all parameter MODE values for CLATMR (+-1,..,+-6)
DO MODE = 1, 6
ALLOCATE( IWORK(2*M) )
ALLOCATE( CDA(M) )
ALLOCATE( CDL(M) )
ALLOCATE( CDR(M) )
CALL CLATMR( M, M, 'N', ISEED, 'N', CDA, MODE, COND, &
CMAX, RSIGN, GRADE, CDL, MODEL, CONDL, &
CDR, MODER, CONDR, PIVTNG, IWORK, M, M, &
ZERO, -ONE, 'N', A, LDA, IWORK(M+1), INFO )
DEALLOCATE( CDR )
DEALLOCATE( CDL )
DEALLOCATE( CDA )
DEALLOCATE( IWORK )
LCWORK = MAX(1,2*M)
ALLOCATE( CEIGSA(M) )
ALLOCATE( CWORK(LCWORK) )
ALLOCATE( WORK(2*M) )
AC(1:M,1:M) = A(1:M,1:M)
CALL CGEEV( 'N','N', M, AC, LDA, CEIGSA, CDUM2X2, 2, &
CDUM2X2, 2, CWORK, LCWORK, WORK, INFO ) ! LAPACK CALL
DEALLOCATE(WORK)
DEALLOCATE(CWORK)
TMP = ABS(CEIGSA(ICAMAX(M, CEIGSA, 1))) ! The spectral radius of A
! Scale the matrix A to have unit spectral radius.
CALL CLASCL( 'G',0, 0, TMP, ONE, M, M, &
A, LDA, INFO )
CALL CLASCL( 'G',0, 0, TMP, ONE, M, 1, &
CEIGSA, M, INFO )
ANORM = CLANGE( 'F', M, M, A, LDA, WDUMMY )
IF ( K_traj == 2 ) THEN
! generate data as two trajectories
! with two inital conditions
CALL CLARNV(2, ISEED, M, F(1,1) )
DO i = 1, N/2
CALL CGEMV( 'N', M, M, CONE, A, LDA, F(1,i), 1, &
CZERO, F(1,i+1), 1 )
END DO
X0(1:M,1:N/2) = F(1:M,1:N/2)
Y0(1:M,1:N/2) = F(1:M,2:N/2+1)
CALL CLARNV(2, ISEED, M, F(1,1) )
DO i = 1, N-N/2
CALL CGEMV( 'N', M, M, CONE, A, LDA, F(1,i), 1, &
CZERO, F(1,i+1), 1 )
END DO
X0(1:M,N/2+1:N) = F(1:M,1:N-N/2)
Y0(1:M,N/2+1:N) = F(1:M,2:N-N/2+1)
ELSE
CALL CLARNV(2, ISEED, M, F(1,1) )
DO i = 1, N
CALL CGEMV( 'N', M, M, CONE, A, M, F(1,i), 1, &
CZERO, F(1,i+1), 1 )
END DO
F0(1:M,1:N+1) = F(1:M,1:N+1)
X0(1:M,1:N) = F0(1:M,1:N)
Y0(1:M,1:N) = F0(1:M,2:N+1)
END IF
DEALLOCATE( CEIGSA )
!........................................................................
DO iJOBZ = 1, 4
SELECT CASE ( iJOBZ )
CASE(1)
JOBZ = 'V'
RESIDS = 'R'
CASE(2)
JOBZ = 'V'
RESIDS = 'N'
CASE(3)
JOBZ = 'F'
RESIDS = 'N'
CASE(4)
JOBZ = 'N'
RESIDS = 'N'
END SELECT
DO iJOBREF = 1, 3
SELECT CASE ( iJOBREF )
CASE(1)
JOBREF = 'R'
CASE(2)
JOBREF = 'E'
CASE(3)
JOBREF = 'N'
END SELECT
DO iSCALE = 1, 4
SELECT CASE ( iSCALE )
CASE(1)
SCALE = 'S'
CASE(2)
SCALE = 'C'
CASE(3)
SCALE = 'Y'
CASE(4)
SCALE = 'N'
END SELECT
DO iNRNK = -1, -2, -1
NRNK = iNRNK
DO iWHTSVD = 1, 3
! Check all four options to compute the POD basis
! via the SVD.
WHTSVD = iWHTSVD
DO LWMINOPT = 1, 2
! Workspace query for the minimal (1) and for the optimal
! (2) workspace lengths determined by workspace query.
! CGEDMD is always tested and its results are also used for
! comparisons with CGEDMDQ.
X(1:M,1:N) = X0(1:M,1:N)
Y(1:M,1:N) = Y0(1:M,1:N)
CALL CGEDMD( SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, X, LDX, Y, LDY, NRNK, TOL, &
K, CEIGS, Z, LDZ, RES, &
AU, LDAU, W, LDW, S, LDS, &
CDUMMY, -1, WDUMMY, -1, IDUMMY, -1, INFO )
IF ( (INFO .EQ. 2) .OR. ( INFO .EQ. 3 ) &
.OR. ( INFO < 0 ) ) THEN
WRITE(*,*) 'Call to CGEDMD workspace query failed. &
&Check the calling sequence and the code.'
WRITE(*,*) 'The error code is ', INFO
WRITE(*,*) 'The input parameters were ', &
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL, LDZ, LDAU, LDW, LDS
STOP
ELSE
!WRITE(*,*) '... done. Workspace length computed.'
END IF
LCWORK = INT(CDUMMY(LWMINOPT))
ALLOCATE(CWORK(LCWORK))
LIWORK = IDUMMY(1)
ALLOCATE(IWORK(LIWORK))
LWORK = INT(WDUMMY(1))
ALLOCATE(WORK(LWORK))
CALL CGEDMD( SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, X, LDX, Y, LDY, NRNK, TOL, &
K, CEIGS, Z, LDZ, RES, &
AU, LDAU, W, LDW, S, LDS, &
CWORK, LCWORK, WORK, LWORK, IWORK, LIWORK, INFO )
IF ( INFO /= 0 ) THEN
WRITE(*,*) 'Call to CGEDMD failed. &
&Check the calling sequence and the code.'
WRITE(*,*) 'The error code is ', INFO
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
STOP
END IF
SINGVX(1:N) = WORK(1:N)
!...... CGEDMD check point
IF ( LSAME(JOBZ,'V') ) THEN
! Check that Z = X*W, on return from CGEDMD
! This checks that the returned eigenvectors in Z are
! the product of the SVD'POD basis returned in X
! and the eigenvectors of the Rayleigh quotient
! returned in W
CALL CGEMM( 'N', 'N', M, K, K, CONE, X, LDX, W, LDW, &
CZERO, Z1, LDZ )
TMP = ZERO
DO i = 1, K
CALL CAXPY( M, -CONE, Z(1,i), 1, Z1(1,i), 1)
TMP = MAX(TMP, SCNRM2( M, Z1(1,i), 1 ) )
END DO
TMP_XW = MAX(TMP_XW, TMP )
IF ( TMP_XW <= TOL ) THEN
!WRITE(*,*) ' :) .... OK .........CGEDMD PASSED.'
ELSE
NFAIL_Z_XV = NFAIL_Z_XV + 1
WRITE(*,*) ':( .................CGEDMD FAILED!', &
'Check the code for implementation errors.'
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
END IF
END IF
!...... CGEDMD check point
IF ( LSAME(JOBREF,'R') ) THEN
! The matrix A*U is returned for computing refined Ritz vectors.
! Check that A*U is computed correctly using the formula
! A*U = Y * V * inv(SIGMA). This depends on the
! accuracy in the computed singular values and vectors of X.
! See the paper for an error analysis.
! Note that the left singular vectors of the input matrix X
! are returned in the array X.
CALL CGEMM( 'N', 'N', M, K, M, CONE, A, LDA, X, LDX, &
CZERO, Z1, LDZ )
TMP = ZERO
DO i = 1, K
CALL CAXPY( M, -CONE, AU(1,i), 1, Z1(1,i), 1)
TMP = MAX( TMP, SCNRM2( M, Z1(1,i),1 ) * &
SINGVX(K)/(ANORM*SINGVX(1)) )
END DO
TMP_AU = MAX( TMP_AU, TMP )
IF ( TMP <= TOL2 ) THEN
!WRITE(*,*) ':) .... OK .........CGEDMD PASSED.'
ELSE
NFAIL_AU = NFAIL_AU + 1
WRITE(*,*) ':( .................CGEDMD FAILED!', &
'Check the code for implementation errors.'
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL2
END IF
ELSEIF ( LSAME(JOBREF,'E') ) THEN
! The unscaled vectors of the Exact DMD are computed.
! This option is included for the sake of completeness,
! for users who prefer the Exact DMD vectors. The
! returned vectors are in the real form, in the same way
! as the Ritz vectors. Here we just save the vectors
! and test them separately using a Matlab script.
CALL CGEMM( 'N', 'N', M, K, M, CONE, A, LDA, AU, LDAU, CZERO, Y1, LDY )
DO i=1, K
CALL CAXPY( M, -CEIGS(i), AU(1,i), 1, Y1(1,i), 1 )
RESEX(i) = SCNRM2( M, Y1(1,i), 1) / SCNRM2(M,AU(1,i),1)
END DO
END IF
!...... CGEDMD check point
IF ( LSAME(RESIDS, 'R') ) THEN
! Compare the residuals returned by CGEDMD with the
! explicitly computed residuals using the matrix A.
! Compute explicitly Y1 = A*Z
CALL CGEMM( 'N', 'N', M, K, M, CONE, A, LDA, Z, LDZ, CZERO, Y1, LDY )
! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
! of the invariant subspaces that correspond to complex conjugate
! pairs of eigencalues. (See the description of Z in CGEDMD,)
DO i=1, K
! have a real eigenvalue with real eigenvector
CALL CAXPY( M, -CEIGS(i), Z(1,i), 1, Y1(1,i), 1 )
RES1(i) = SCNRM2( M, Y1(1,i), 1)
END DO
TMP = ZERO
DO i = 1, K
TMP = MAX( TMP, ABS(RES(i) - RES1(i)) * &
SINGVX(K)/(ANORM*SINGVX(1)) )
END DO
TMP_REZ = MAX( TMP_REZ, TMP )
IF ( TMP <= TOL2 ) THEN
!WRITE(*,*) ':) .... OK ..........CGEDMD PASSED.'
ELSE
NFAIL_REZ = NFAIL_REZ + 1
WRITE(*,*) ':( ..................CGEDMD FAILED!', &
'Check the code for implementation errors.'
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
END IF
IF ( LSAME(JOBREF,'E') ) THEN
TMP = ZERO
DO i = 1, K
TMP = MAX( TMP, ABS(RES1(i) - RESEX(i))/(RES1(i)+RESEX(i)) )
END DO
TMP_EX = MAX(TMP_EX,TMP)
END IF
END IF
DEALLOCATE(CWORK)
DEALLOCATE(WORK)
DEALLOCATE(IWORK)
!.......................................................................................................
IF ( K_traj == 1 ) THEN
F(1:M,1:N+1) = F0(1:M,1:N+1)
CALL CGEDMDQ( SCALE, JOBZ, RESIDS, WANTQ, WANTR, JOBREF, &
WHTSVD, M, N+1, F, LDF, X, LDX, Y, LDY, &
NRNK, TOL, K, CEIGS, Z, LDZ, RES, AU, &
LDAU, W, LDW, S, LDS, CDUMMY, -1, &
WDUMMY, -1, IDUMMY, -1, INFO )
LCWORK = INT(CDUMMY(LWMINOPT))
ALLOCATE(CWORK(LCWORK))
LIWORK = IDUMMY(1)
ALLOCATE(IWORK(LIWORK))
LWORK = INT(WDUMMY(1))
ALLOCATE(WORK(LWORK))
CALL CGEDMDQ( SCALE, JOBZ, RESIDS, WANTQ, WANTR, JOBREF, &
WHTSVD, M, N+1, F, LDF, X, LDX, Y, LDY, &
NRNK, TOL, KQ, CEIGS, Z, LDZ, RES, AU, &
LDAU, W, LDW, S, LDS, CWORK, LCWORK, &
WORK, LWORK, IWORK, LIWORK, INFO )
IF ( INFO /= 0 ) THEN
WRITE(*,*) 'Call to CGEDMDQ failed. &
&Check the calling sequence and the code.'
WRITE(*,*) 'The error code is ', INFO
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, WANTQ, WANTR, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
STOP
END IF
SINGVQX(1:N) =WORK(1:N)
!..... ZGEDMDQ check point
TMP = ZERO
DO i = 1, MIN(K, KQ)
TMP = MAX(TMP, ABS(SINGVX(i)-SINGVQX(i)) / &
SINGVX(1) )
END DO
SVDIFF = MAX( SVDIFF, TMP )
IF ( TMP > TOL2 ) THEN
WRITE(*,*) 'FAILED! Something was wrong with the run.'
NFAIL_SVDIFF = NFAIL_SVDIFF + 1
END IF
!..... CGEDMDQ check point
!..... CGEDMDQ check point
IF ( LSAME(WANTQ,'Q') .AND. LSAME(WANTR,'R') ) THEN
! Check that the QR factors are computed and returned
! as requested. The residual ||F-Q*R||_F / ||F||_F
! is compared to M*N*EPS.
F1(1:M,1:N+1) = F0(1:M,1:N+1)
CALL CGEMM( 'N', 'N', M, N+1, MIN(M,N+1), -CONE, F, &
LDF, Y, LDY, CONE, F1, LDF )
TMP_FQR = CLANGE( 'F', M, N+1, F1, LDF, WORK ) / &
CLANGE( 'F', M, N+1, F0, LDF, WORK )
IF ( TMP_FQR <= TOL2 ) THEN
!WRITE(*,*) ':) CGEDMDQ ........ PASSED.'
ELSE
WRITE(*,*) ':( CGEDMDQ ........ FAILED.'
NFAIL_F_QR = NFAIL_F_QR + 1
END IF
END IF
!..... ZGEDMDQ checkpoint
!..... ZGEDMDQ checkpoint
IF ( LSAME(RESIDS, 'R') ) THEN
! Compare the residuals returned by ZGEDMDQ with the
! explicitly computed residuals using the matrix A.
! Compute explicitly Y1 = A*Z
CALL CGEMM( 'N', 'N', M, KQ, M, CONE, A, LDA, Z, LDZ, CZERO, Y1, LDY )
! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
! of the invariant subspaces that correspond to complex conjugate
! pairs of eigencalues. (See the description of Z in ZGEDMDQ)
DO i = 1, KQ
! have a real eigenvalue with real eigenvector
CALL CAXPY( M, -CEIGS(i), Z(1,i), 1, Y1(1,i), 1 )
! Y(1:M,i) = Y(1:M,i) - REIG(i)*Z(1:M,i)
RES1(i) = SCNRM2( M, Y1(1,i), 1)
END DO
TMP = ZERO
DO i = 1, KQ
TMP = MAX( TMP, ABS(RES(i) - RES1(i)) * &
SINGVQX(KQ)/(ANORM*SINGVQX(1)) )
END DO
TMP_REZQ = MAX( TMP_REZQ, TMP )
IF ( TMP <= TOL2 ) THEN
!WRITE(*,*) '.... OK ........ CGEDMDQ PASSED.'
ELSE
NFAIL_REZQ = NFAIL_REZQ + 1
WRITE(*,*) '................ CGEDMDQ FAILED!', &
'Check the code for implementation errors.'
END IF
END IF
DEALLOCATE(CWORK)
DEALLOCATE(WORK)
DEALLOCATE(IWORK)
END IF
END DO ! LWMINOPT
!write(*,*) 'LWMINOPT loop completed'
END DO ! iWHTSVD
!write(*,*) 'WHTSVD loop completed'
END DO ! iNRNK -2:-1
!write(*,*) 'NRNK loop completed'
END DO ! iSCALE 1:4
!write(*,*) 'SCALE loop completed'
END DO
!write(*,*) 'JOBREF loop completed'
END DO ! iJOBZ
!write(*,*) 'JOBZ loop completed'
END DO ! MODE -6:6
!write(*,*) 'MODE loop completed'
END DO ! 1 or 2 trajectories
!write(*,*) 'trajectories loop completed'
DEALLOCATE( A )
DEALLOCATE( AC )
DEALLOCATE( Z )
DEALLOCATE( F )
DEALLOCATE( F0 )
DEALLOCATE( F1 )
DEALLOCATE( X )
DEALLOCATE( X0 )
DEALLOCATE( Y )
DEALLOCATE( Y0 )
DEALLOCATE( Y1 )
DEALLOCATE( AU )
DEALLOCATE( W )
DEALLOCATE( S )
DEALLOCATE( Z1 )
DEALLOCATE( RES )
DEALLOCATE( RES1 )
DEALLOCATE( RESEX )
DEALLOCATE( CEIGS )
DEALLOCATE( SINGVX )
DEALLOCATE( SINGVQX )
END DO ! LLOOP
WRITE(*,*)
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*) ' Test summary for CGEDMD :'
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*)
IF ( NFAIL_Z_XV == 0 ) THEN
WRITE(*,*) '>>>> Z - U*V test PASSED.'
ELSE
WRITE(*,*) 'Z - U*V test FAILED ', NFAIL_Z_XV, ' time(s)'
WRITE(*,*) 'Max error ||Z-U*V||_F was ', TMP_XW
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_z_XV
END IF
IF ( NFAIL_AU == 0 ) THEN
WRITE(*,*) '>>>> A*U test PASSED. '
ELSE
WRITE(*,*) 'A*U test FAILED ', NFAIL_AU, ' time(s)'
WRITE(*,*) 'Max A*U test adjusted error measure was ', TMP_AU
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_AU
END IF
IF ( NFAIL_REZ == 0 ) THEN
WRITE(*,*) '>>>> Rezidual computation test PASSED.'
ELSE
WRITE(*,*) 'Rezidual computation test FAILED ', NFAIL_REZ, 'time(s)'
WRITE(*,*) 'Max residual computing test adjusted error measure was ', TMP_REZ
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_REZ
END IF
IF ( NFAIL_TOTAL == 0 ) THEN
WRITE(*,*) '>>>> CGEDMD :: ALL TESTS PASSED.'
ELSE
WRITE(*,*) NFAIL_TOTAL, 'FAILURES!'
WRITE(*,*) '>>>>>>>>>>>>>> CGEDMD :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
END IF
WRITE(*,*)
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*) ' Test summary for CGEDMDQ :'
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*)
IF ( NFAIL_SVDIFF == 0 ) THEN
WRITE(*,*) '>>>> CGEDMD and CGEDMDQ computed singular &
&values test PASSED.'
ELSE
WRITE(*,*) 'ZGEDMD and ZGEDMDQ discrepancies in &
&the singular values unacceptable ', &
NFAIL_SVDIFF, ' times. Test FAILED.'
WRITE(*,*) 'The maximal discrepancy in the singular values (relative to the norm) was ', SVDIFF
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_SVDIFF
END IF
IF ( NFAIL_F_QR == 0 ) THEN
WRITE(*,*) '>>>> F - Q*R test PASSED.'
ELSE
WRITE(*,*) 'F - Q*R test FAILED ', NFAIL_F_QR, ' time(s)'
WRITE(*,*) 'The largest relative residual was ', TMP_FQR
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_F_QR
END IF
IF ( NFAIL_REZQ == 0 ) THEN
WRITE(*,*) '>>>> Rezidual computation test PASSED.'
ELSE
WRITE(*,*) 'Rezidual computation test FAILED ', NFAIL_REZQ, 'time(s)'
WRITE(*,*) 'Max residual computing test adjusted error measure was ', TMP_REZQ
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_REZQ
END IF
IF ( NFAILQ_TOTAL == 0 ) THEN
WRITE(*,*) '>>>>>>> CGEDMDQ :: ALL TESTS PASSED.'
ELSE
WRITE(*,*) NFAILQ_TOTAL, 'FAILURES!'
WRITE(*,*) '>>>>>>> CGEDMDQ :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
END IF
WRITE(*,*)
WRITE(*,*) 'Test completed.'
STOP
END

+ 813
- 0
lapack-netlib/TESTING/EIG/dchkdmd.f90 View File

@@ -0,0 +1,813 @@
! This is a test program for checking the implementations of
! the implementations of the following subroutines
!
! DGEDMD for computation of the
! Dynamic Mode Decomposition (DMD)
! DGEDMDQ for computation of a
! QR factorization based compressed DMD
!
! Developed and supported by:
! ===========================
! Developed and coded by Zlatko Drmac, Faculty of Science,
! University of Zagreb; drmac@math.hr
! In cooperation with
! AIMdyn Inc., Santa Barbara, CA.
! ========================================================
! How to run the code (compiler, link info)
! ========================================================
! Compile as FORTRAN 90 (or later) and link with BLAS and
! LAPACK libraries.
! NOTE: The code is developed and tested on top of the
! Intel MKL library (versions 2022.0.3 and 2022.2.0),
! using the Intel Fortran compiler.
!
! For developers of the C++ implementation
! ========================================================
! See the LAPACK++ and Template Numerical Toolkit (TNT)
!
! Note on a development of the GPU HP implementation
! ========================================================
! Work in progress. See CUDA, MAGMA, SLATE.
! NOTE: The four SVD subroutines used in this code are
! included as a part of R&D and for the completeness.
! This was also an opportunity to test those SVD codes.
! If the scaling option is used all four are essentially
! equally good. For implementations on HP platforms,
! one can use whichever SVD is available.
!... .........................................................
! NOTE:
! When using the Intel MKL 2022.0.3 the subroutine xGESVDQ
! (optionally used in xGEDMD) may cause access violation
! error for x = S, D, C, Z, but only if called with the
! work space query. (At least in our Windows 10 MSVS 2019.)
! The problem can be mitigated by downloading the source
! code of xGESVDQ from the LAPACK repository and use it
! localy instead of the one in the MKL. This seems to
! indicate that the problem is indeed in the MKL.
! This problem did not appear whith Intel MKL 2022.2.0.
!
! NOTE:
! xGESDD seems to have a problem with workspace. In some
! cases the length of the optimal workspace is returned
! smaller than the minimal workspace, as specified in the
! code. As a precaution, all optimal workspaces are
! set as MAX(minimal, optimal).
! Latest implementations of complex xGESDD have different
! length of the real worksapce. We use max value over
! two versions.
!............................................................
!............................................................
!
PROGRAM DMD_TEST
use iso_fortran_env, only: real64
IMPLICIT NONE
integer, parameter :: WP = real64
!............................................................
REAL(KIND=WP), PARAMETER :: ONE = 1.0_WP
REAL(KIND=WP), PARAMETER :: ZERO = 0.0_WP
!............................................................
REAL(KIND=WP), ALLOCATABLE, DIMENSION(:,:) :: &
A, AC, EIGA, LAMBDA, LAMBDAQ, F, F1, F2,&
Z, Z1, S, AU, W, VA, X, X0, Y, Y0, Y1
REAL(KIND=WP), ALLOCATABLE, DIMENSION(:) :: &
DA, DL, DR, REIG, REIGA, REIGQ, IEIG, &
IEIGA, IEIGQ, RES, RES1, RESEX, SINGVX,&
SINGVQX, WORK
INTEGER , ALLOCATABLE, DIMENSION(:) :: IWORK
REAL(KIND=WP) :: AB(2,2), WDUMMY(2)
INTEGER :: IDUMMY(2), ISEED(4), RJOBDATA(8)
REAL(KIND=WP) :: ANORM, COND, CONDL, CONDR, DMAX, EPS, &
TOL, TOL2, SVDIFF, TMP, TMP_AU, &
TMP_FQR, TMP_REZ, TMP_REZQ, TMP_ZXW, &
TMP_EX, XNORM, YNORM
!............................................................
INTEGER :: K, KQ, LDF, LDS, LDA, LDAU, LDW, LDX, LDY, &
LDZ, LIWORK, LWORK, M, N, L, LLOOP, NRNK
INTEGER :: i, iJOBREF, iJOBZ, iSCALE, INFO, j, KDIFF, &
NFAIL, NFAIL_AU, NFAIL_F_QR, NFAIL_REZ, &
NFAIL_REZQ, NFAIL_SVDIFF, NFAIL_TOTAL, NFAILQ_TOTAL, &
NFAIL_Z_XV, MODE, MODEL, MODER, WHTSVD
INTEGER iNRNK, iWHTSVD, K_TRAJ, LWMINOPT
CHARACTER(LEN=1) GRADE, JOBREF, JOBZ, PIVTNG, RSIGN, &
SCALE, RESIDS, WANTQ, WANTR
LOGICAL TEST_QRDMD
!..... external subroutines (BLAS and LAPACK)
EXTERNAL DAXPY, DGEEV, DGEMM, DGEMV, DLACPY, DLASCL
EXTERNAL DLARNV, DLATMR
!.....external subroutines DMD package, part 1
! subroutines under test
EXTERNAL DGEDMD, DGEDMDQ
!..... external functions (BLAS and LAPACK)
EXTERNAL DLAMCH, DLANGE, DNRM2
REAL(KIND=WP) :: DLAMCH, DLANGE, DNRM2
EXTERNAL LSAME
LOGICAL LSAME
INTRINSIC ABS, INT, MIN, MAX
!............................................................
! The test is always in pairs : ( DGEDMD and DGEDMDQ )
! because the test includes comparing the results (in pairs).
!.....................................................................................
TEST_QRDMD = .TRUE. ! This code by default performs tests on DGEDMDQ
! Since the QR factorizations based algorithm is designed for
! single trajectory data, only single trajectory tests will
! be performed with xGEDMDQ.
WANTQ = 'Q'
WANTR = 'R'
!.................................................................................
EPS = DLAMCH( 'P' ) ! machine precision DP
! Global counters of failures of some particular tests
NFAIL = 0
NFAIL_REZ = 0
NFAIL_REZQ = 0
NFAIL_Z_XV = 0
NFAIL_F_QR = 0
NFAIL_AU = 0
KDIFF = 0
NFAIL_SVDIFF = 0
NFAIL_TOTAL = 0
NFAILQ_TOTAL = 0
DO LLOOP = 1, 4
WRITE(*,*) 'L Loop Index = ', LLOOP
! Set the dimensions of the problem ...
WRITE(*,*) 'M = '
READ(*,*) M
WRITE(*,*) M
! ... and the number of snapshots.
WRITE(*,*) 'N = '
READ(*,*) N
WRITE(*,*) N
! ... Test the dimensions
IF ( ( MIN(M,N) == 0 ) .OR. ( M < N ) ) THEN
WRITE(*,*) 'Bad dimensions. Required: M >= N > 0.'
STOP
END IF
!.............
! The seed inside the LLOOP so that each pass can be reproduced easily.
ISEED(1) = 4
ISEED(2) = 3
ISEED(3) = 2
ISEED(4) = 1
LDA = M
LDF = M
LDX = MAX(M,N+1)
LDY = MAX(M,N+1)
LDW = N
LDZ = M
LDAU = MAX(M,N+1)
LDS = N
TMP_ZXW = ZERO
TMP_AU = ZERO
TMP_REZ = ZERO
TMP_REZQ = ZERO
SVDIFF = ZERO
TMP_EX = ZERO
!
! Test the subroutines on real data snapshots. All
! computation is done in real arithmetic, even when
! Koopman eigenvalues and modes are real.
!
! Allocate memory space
ALLOCATE( A(LDA,M) )
ALLOCATE( AC(LDA,M) )
ALLOCATE( DA(M) )
ALLOCATE( DL(M) )
ALLOCATE( F(LDF,N+1) )
ALLOCATE( F1(LDF,N+1) )
ALLOCATE( F2(LDF,N+1) )
ALLOCATE( X(LDX,N) )
ALLOCATE( X0(LDX,N) )
ALLOCATE( SINGVX(N) )
ALLOCATE( SINGVQX(N) )
ALLOCATE( Y(LDY,N+1) )
ALLOCATE( Y0(LDY,N+1) )
ALLOCATE( Y1(M,N+1) )
ALLOCATE( Z(LDZ,N) )
ALLOCATE( Z1(LDZ,N) )
ALLOCATE( RES(N) )
ALLOCATE( RES1(N) )
ALLOCATE( RESEX(N) )
ALLOCATE( REIG(N) )
ALLOCATE( IEIG(N) )
ALLOCATE( REIGQ(N) )
ALLOCATE( IEIGQ(N) )
ALLOCATE( REIGA(M) )
ALLOCATE( IEIGA(M) )
ALLOCATE( VA(LDA,M) )
ALLOCATE( LAMBDA(N,2) )
ALLOCATE( LAMBDAQ(N,2) )
ALLOCATE( EIGA(M,2) )
ALLOCATE( W(LDW,N) )
ALLOCATE( AU(LDAU,N) )
ALLOCATE( S(N,N) )
TOL = M*EPS
! This mimics O(M*N)*EPS bound for accumulated roundoff error.
! The factor 10 is somewhat arbitrary.
TOL2 = 10*M*N*EPS
!.............
DO K_TRAJ = 1, 2
! Number of intial conditions in the simulation/trajectories (1 or 2)
COND = 1.0D8
DMAX = 1.0D2
RSIGN = 'F'
GRADE = 'N'
MODEL = 6
CONDL = 1.0D2
MODER = 6
CONDR = 1.0D2
PIVTNG = 'N'
! Loop over all parameter MODE values for ZLATMR (+1,..,+6)
DO MODE = 1, 6
ALLOCATE( IWORK(2*M) )
ALLOCATE(DR(N))
CALL DLATMR( M, M, 'S', ISEED, 'N', DA, MODE, COND, &
DMAX, RSIGN, GRADE, DL, MODEL, CONDL, &
DR, MODER, CONDR, PIVTNG, IWORK, M, M, &
ZERO, -ONE, 'N', A, LDA, IWORK(M+1), INFO )
DEALLOCATE(IWORK)
DEALLOCATE(DR)
LWORK = 4*M+1
ALLOCATE(WORK(LWORK))
AC = A
CALL DGEEV( 'N','V', M, AC, M, REIGA, IEIGA, VA, M, &
VA, M, WORK, LWORK, INFO ) ! LAPACK CALL
DEALLOCATE(WORK)
TMP = ZERO
DO i = 1, M
EIGA(i,1) = REIGA(i)
EIGA(i,2) = IEIGA(i)
TMP = MAX( TMP, SQRT(REIGA(i)**2+IEIGA(i)**2))
END DO
! Scale A to have the desirable spectral radius.
CALL DLASCL( 'G', 0, 0, TMP, ONE, M, M, A, M, INFO )
CALL DLASCL( 'G', 0, 0, TMP, ONE, M, 2, EIGA, M, INFO )
! Compute the norm of A
ANORM = DLANGE( 'F', N, N, A, M, WDUMMY )
IF ( K_TRAJ == 2 ) THEN
! generate data with two inital conditions
CALL DLARNV(2, ISEED, M, F1(1,1) )
F1(1:M,1) = 1.0E-10*F1(1:M,1)
DO i = 1, N/2
CALL DGEMV( 'N', M, M, ONE, A, M, F1(1,i), 1, ZERO, &
F1(1,i+1), 1 )
END DO
X0(1:M,1:N/2) = F1(1:M,1:N/2)
Y0(1:M,1:N/2) = F1(1:M,2:N/2+1)
CALL DLARNV(2, ISEED, M, F1(1,1) )
DO i = 1, N-N/2
CALL DGEMV( 'N', M, M, ONE, A, M, F1(1,i), 1, ZERO, &
F1(1,i+1), 1 )
END DO
X0(1:M,N/2+1:N) = F1(1:M,1:N-N/2)
Y0(1:M,N/2+1:N) = F1(1:M,2:N-N/2+1)
ELSE
CALL DLARNV(2, ISEED, M, F(1,1) )
DO i = 1, N
CALL DGEMV( 'N', M, M, ONE, A, M, F(1,i), 1, ZERO, &
F(1,i+1), 1 )
END DO
X0(1:M,1:N) = F(1:M,1:N)
Y0(1:M,1:N) = F(1:M,2:N+1)
END IF
XNORM = DLANGE( 'F', M, N, X0, LDX, WDUMMY )
YNORM = DLANGE( 'F', M, N, Y0, LDX, WDUMMY )
!............................................................
DO iJOBZ = 1, 4
SELECT CASE ( iJOBZ )
CASE(1)
JOBZ = 'V' ! Ritz vectors will be computed
RESIDS = 'R' ! Residuals will be computed
CASE(2)
JOBZ = 'V'
RESIDS = 'N'
CASE(3)
JOBZ = 'F' ! Ritz vectors in factored form
RESIDS = 'N'
CASE(4)
JOBZ = 'N'
RESIDS = 'N'
END SELECT
DO iJOBREF = 1, 3
SELECT CASE ( iJOBREF )
CASE(1)
JOBREF = 'R' ! Data for refined Ritz vectors
CASE(2)
JOBREF = 'E' ! Exact DMD vectors
CASE(3)
JOBREF = 'N'
END SELECT
DO iSCALE = 1, 4
SELECT CASE ( iSCALE )
CASE(1)
SCALE = 'S' ! X data normalized
CASE(2)
SCALE = 'C' ! X normalized, consist. check
CASE(3)
SCALE = 'Y' ! Y data normalized
CASE(4)
SCALE = 'N'
END SELECT
DO iNRNK = -1, -2, -1
! Two truncation strategies. The "-2" case for R&D
! purposes only - it uses possibly low accuracy small
! singular values, in which case the formulas used in
! the DMD are highly sensitive.
NRNK = iNRNK
DO iWHTSVD = 1, 4
! Check all four options to compute the POD basis
! via the SVD.
WHTSVD = iWHTSVD
DO LWMINOPT = 1, 2
! Workspace query for the minimal (1) and for the optimal
! (2) workspace lengths determined by workspace query.
X(1:M,1:N) = X0(1:M,1:N)
Y(1:M,1:N) = Y0(1:M,1:N)
! DGEDMD: Workspace query and workspace allocation
CALL DGEDMD( SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, M, &
N, X, LDX, Y, LDY, NRNK, TOL, K, REIG, IEIG, Z, &
LDZ, RES, AU, LDAU, W, LDW, S, LDS, WDUMMY, -1, &
IDUMMY, -1, INFO )
LIWORK = IDUMMY(1)
ALLOCATE( IWORK(LIWORK) )
LWORK = INT(WDUMMY(LWMINOPT))
ALLOCATE( WORK(LWORK) )
! DGEDMD test: CALL DGEDMD
CALL DGEDMD( SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, M, &
N, X, LDX, Y, LDY, NRNK, TOL, K, REIG, IEIG, Z, &
LDZ, RES, AU, LDAU, W, LDW, S, LDS, WORK, LWORK,&
IWORK, LIWORK, INFO )
SINGVX(1:N) = WORK(1:N)
!...... DGEDMD check point
IF ( LSAME(JOBZ,'V') ) THEN
! Check that Z = X*W, on return from DGEDMD
! This checks that the returned aigenvectors in Z are
! the product of the SVD'POD basis returned in X
! and the eigenvectors of the rayleigh quotient
! returned in W
CALL DGEMM( 'N', 'N', M, K, K, ONE, X, LDX, W, LDW, &
ZERO, Z1, LDZ )
TMP = ZERO
DO i = 1, K
CALL DAXPY( M, -ONE, Z(1,i), 1, Z1(1,i), 1)
TMP = MAX(TMP, DNRM2( M, Z1(1,i), 1 ) )
END DO
TMP_ZXW = MAX(TMP_ZXW, TMP )
IF ( TMP_ZXW > 10*M*EPS ) THEN
NFAIL_Z_XV = NFAIL_Z_XV + 1
WRITE(*,*) ':( .................DGEDMD FAILED!', &
'Check the code for implementation errors.'
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
END IF
END IF
!...... DGEDMD check point
IF ( LSAME(JOBREF,'R') ) THEN
! The matrix A*U is returned for computing refined Ritz vectors.
! Check that A*U is computed correctly using the formula
! A*U = Y * V * inv(SIGMA). This depends on the
! accuracy in the computed singular values and vectors of X.
! See the paper for an error analysis.
! Note that the left singular vectors of the input matrix X
! are returned in the array X.
CALL DGEMM( 'N', 'N', M, K, M, ONE, A, LDA, X, LDX, &
ZERO, Z1, LDZ )
TMP = ZERO
DO i = 1, K
CALL DAXPY( M, -ONE, AU(1,i), 1, Z1(1,i), 1)
TMP = MAX( TMP, DNRM2( M, Z1(1,i),1 ) * &
SINGVX(K)/(ANORM*SINGVX(1)) )
END DO
TMP_AU = MAX( TMP_AU, TMP )
IF ( TMP > TOL2 ) THEN
NFAIL_AU = NFAIL_AU + 1
WRITE(*,*) ':( .................DGEDMD FAILED!', &
'Check the code for implementation errors.'
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
END IF
ELSEIF ( LSAME(JOBREF,'E') ) THEN
! The unscaled vectors of the Exact DMD are computed.
! This option is included for the sake of completeness,
! for users who prefer the Exact DMD vectors. The
! returned vectors are in the real form, in the same way
! as the Ritz vectors. Here we just save the vectors
! and test them separately using a Matlab script.
CALL DGEMM( 'N', 'N', M, K, M, ONE, A, LDA, AU, LDAU, ZERO, Y1, M )
i=1
DO WHILE ( i <= K )
IF ( IEIG(i) == ZERO ) THEN
! have a real eigenvalue with real eigenvector
CALL DAXPY( M, -REIG(i), AU(1,i), 1, Y1(1,i), 1 )
RESEX(i) = DNRM2( M, Y1(1,i), 1) / DNRM2(M,AU(1,i),1)
i = i + 1
ELSE
! Have a complex conjugate pair
! REIG(i) +- sqrt(-1)*IMEIG(i).
! Since all computation is done in real
! arithmetic, the formula for the residual
! is recast for real representation of the
! complex conjugate eigenpair. See the
! description of RES.
AB(1,1) = REIG(i)
AB(2,1) = -IEIG(i)
AB(1,2) = IEIG(i)
AB(2,2) = REIG(i)
CALL DGEMM( 'N', 'N', M, 2, 2, -ONE, AU(1,i), &
M, AB, 2, ONE, Y1(1,i), M )
RESEX(i) = DLANGE( 'F', M, 2, Y1(1,i), M, &
WORK )/ DLANGE( 'F', M, 2, AU(1,i), M, &
WORK )
RESEX(i+1) = RESEX(i)
i = i + 2
END IF
END DO
END IF
!...... DGEDMD check point
IF ( LSAME(RESIDS, 'R') ) THEN
! Compare the residuals returned by DGEDMD with the
! explicitly computed residuals using the matrix A.
! Compute explicitly Y1 = A*Z
CALL DGEMM( 'N', 'N', M, K, M, ONE, A, LDA, Z, LDZ, ZERO, Y1, M )
! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
! of the invariant subspaces that correspond to complex conjugate
! pairs of eigencalues. (See the description of Z in DGEDMD,)
i = 1
DO WHILE ( i <= K )
IF ( IEIG(i) == ZERO ) THEN
! have a real eigenvalue with real eigenvector
CALL DAXPY( M, -REIG(i), Z(1,i), 1, Y1(1,i), 1 )
RES1(i) = DNRM2( M, Y1(1,i), 1)
i = i + 1
ELSE
! Have a complex conjugate pair
! REIG(i) +- sqrt(-1)*IMEIG(i).
! Since all computation is done in real
! arithmetic, the formula for the residual
! is recast for real representation of the
! complex conjugate eigenpair. See the
! description of RES.
AB(1,1) = REIG(i)
AB(2,1) = -IEIG(i)
AB(1,2) = IEIG(i)
AB(2,2) = REIG(i)
CALL DGEMM( 'N', 'N', M, 2, 2, -ONE, Z(1,i), &
M, AB, 2, ONE, Y1(1,i), M )
RES1(i) = DLANGE( 'F', M, 2, Y1(1,i), M, &
WORK )
RES1(i+1) = RES1(i)
i = i + 2
END IF
END DO
TMP = ZERO
DO i = 1, K
TMP = MAX( TMP, ABS(RES(i) - RES1(i)) * &
SINGVX(K)/(ANORM*SINGVX(1)) )
END DO
TMP_REZ = MAX( TMP_REZ, TMP )
IF ( TMP > TOL2 ) THEN
NFAIL_REZ = NFAIL_REZ + 1
WRITE(*,*) ':( ..................DGEDMD FAILED!', &
'Check the code for implementation errors.'
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
END IF
IF ( LSAME(JOBREF,'E') ) THEN
TMP = ZERO
DO i = 1, K
TMP = MAX( TMP, ABS(RES1(i) - RESEX(i))/(RES1(i)+RESEX(i)) )
END DO
TMP_EX = MAX(TMP_EX,TMP)
END IF
END IF
!..... store the results for inspection
DO i = 1, K
LAMBDA(i,1) = REIG(i)
LAMBDA(i,2) = IEIG(i)
END DO
DEALLOCATE(IWORK)
DEALLOCATE(WORK)
!======================================================================
! Now test the DGEDMDQ
!======================================================================
IF ( TEST_QRDMD .AND. (K_TRAJ == 1) ) THEN
RJOBDATA(2) = 1
F1 = F
! DGEDMDQ test: Workspace query and workspace allocation
CALL DGEDMDQ( SCALE, JOBZ, RESIDS, WANTQ, WANTR, &
JOBREF, WHTSVD, M, N+1, F1, LDF, X, LDX, Y, &
LDY, NRNK, TOL, KQ, REIGQ, IEIGQ, Z, LDZ, &
RES, AU, LDAU, W, LDW, S, LDS, WDUMMY, &
-1, IDUMMY, -1, INFO )
LIWORK = IDUMMY(1)
ALLOCATE( IWORK(LIWORK) )
LWORK = INT(WDUMMY(LWMINOPT))
ALLOCATE(WORK(LWORK))
! DGEDMDQ test: CALL DGEDMDQ
CALL DGEDMDQ( SCALE, JOBZ, RESIDS, WANTQ, WANTR, &
JOBREF, WHTSVD, M, N+1, F1, LDF, X, LDX, Y, &
LDY, NRNK, TOL, KQ, REIGQ, IEIGQ, Z, LDZ, &
RES, AU, LDAU, W, LDW, S, LDS, &
WORK, LWORK, IWORK, LIWORK, INFO )
SINGVQX(1:KQ) = WORK(MIN(M,N+1)+1: MIN(M,N+1)+KQ)
!..... DGEDMDQ check point
IF ( KQ /= K ) THEN
KDIFF = KDIFF+1
END IF
TMP = ZERO
DO i = 1, MIN(K, KQ)
TMP = MAX(TMP, ABS(SINGVX(i)-SINGVQX(i)) / &
SINGVX(1) )
END DO
SVDIFF = MAX( SVDIFF, TMP )
IF ( TMP > M*N*EPS ) THEN
WRITE(*,*) 'FAILED! Something was wrong with the run.'
NFAIL_SVDIFF = NFAIL_SVDIFF + 1
DO j =1, 3
write(*,*) j, SINGVX(j), SINGVQX(j)
read(*,*)
END DO
END IF
!..... DGEDMDQ check point
IF ( LSAME(WANTQ,'Q') .AND. LSAME(WANTR,'R') ) THEN
! Check that the QR factors are computed and returned
! as requested. The residual ||F-Q*R||_F / ||F||_F
! is compared to M*N*EPS.
F2 = F
CALL DGEMM( 'N', 'N', M, N+1, MIN(M,N+1), -ONE, F1, &
LDF, Y, LDY, ONE, F2, LDF )
TMP_FQR = DLANGE( 'F', M, N+1, F2, LDF, WORK ) / &
DLANGE( 'F', M, N+1, F, LDF, WORK )
IF ( TMP_FQR > TOL2 ) THEN
WRITE(*,*) 'FAILED! Something was wrong with the run.'
NFAIL_F_QR = NFAIL_F_QR + 1
END IF
END IF
!..... DGEDMDQ check point
IF ( LSAME(RESIDS, 'R') ) THEN
! Compare the residuals returned by DGEDMDQ with the
! explicitly computed residuals using the matrix A.
! Compute explicitly Y1 = A*Z
CALL DGEMM( 'N', 'N', M, KQ, M, ONE, A, M, Z, M, ZERO, Y1, M )
! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
! of the invariant subspaces that correspond to complex conjugate
! pairs of eigencalues. (See the description of Z in DGEDMDQ)
i = 1
DO WHILE ( i <= KQ )
IF ( IEIGQ(i) == ZERO ) THEN
! have a real eigenvalue with real eigenvector
CALL DAXPY( M, -REIGQ(i), Z(1,i), 1, Y1(1,i), 1 )
! Y(1:M,i) = Y(1:M,i) - REIG(i)*Z(1:M,i)
RES1(i) = DNRM2( M, Y1(1,i), 1)
i = i + 1
ELSE
! Have a complex conjugate pair
! REIG(i) +- sqrt(-1)*IMEIG(i).
! Since all computation is done in real
! arithmetic, the formula for the residual
! is recast for real representation of the
! complex conjugate eigenpair. See the
! description of RES.
AB(1,1) = REIGQ(i)
AB(2,1) = -IEIGQ(i)
AB(1,2) = IEIGQ(i)
AB(2,2) = REIGQ(i)
CALL DGEMM( 'N', 'N', M, 2, 2, -ONE, Z(1,i), &
M, AB, 2, ONE, Y1(1,i), M ) ! BLAS CALL
! Y(1:M,i:i+1) = Y(1:M,i:i+1) - Z(1:M,i:i+1) * AB ! INTRINSIC
RES1(i) = DLANGE( 'F', M, 2, Y1(1,i), M, &
WORK ) ! LAPACK CALL
RES1(i+1) = RES1(i)
i = i + 2
END IF
END DO
TMP = ZERO
DO i = 1, KQ
TMP = MAX( TMP, ABS(RES(i) - RES1(i)) * &
SINGVQX(K)/(ANORM*SINGVQX(1)) )
END DO
TMP_REZQ = MAX( TMP_REZQ, TMP )
IF ( TMP > TOL2 ) THEN
NFAIL_REZQ = NFAIL_REZQ + 1
WRITE(*,*) '................ DGEDMDQ FAILED!', &
'Check the code for implementation errors.'
STOP
END IF
END IF
DO i = 1, KQ
LAMBDAQ(i,1) = REIGQ(i)
LAMBDAQ(i,2) = IEIGQ(i)
END DO
DEALLOCATE(WORK)
DEALLOCATE(IWORK)
END IF ! TEST_QRDMD
!======================================================================
END DO ! LWMINOPT
!write(*,*) 'LWMINOPT loop completed'
END DO ! WHTSVD LOOP
!write(*,*) 'WHTSVD loop completed'
END DO ! NRNK LOOP
!write(*,*) 'NRNK loop completed'
END DO ! SCALE LOOP
!write(*,*) 'SCALE loop completed'
END DO ! JOBF LOOP
!write(*,*) 'JOBREF loop completed'
END DO ! JOBZ LOOP
!write(*,*) 'JOBZ loop completed'
END DO ! MODE -6:6
!write(*,*) 'MODE loop completed'
END DO ! 1 or 2 trajectories
!write(*,*) 'trajectories loop completed'
DEALLOCATE(A)
DEALLOCATE(AC)
DEALLOCATE(DA)
DEALLOCATE(DL)
DEALLOCATE(F)
DEALLOCATE(F1)
DEALLOCATE(F2)
DEALLOCATE(X)
DEALLOCATE(X0)
DEALLOCATE(SINGVX)
DEALLOCATE(SINGVQX)
DEALLOCATE(Y)
DEALLOCATE(Y0)
DEALLOCATE(Y1)
DEALLOCATE(Z)
DEALLOCATE(Z1)
DEALLOCATE(RES)
DEALLOCATE(RES1)
DEALLOCATE(RESEX)
DEALLOCATE(REIG)
DEALLOCATE(IEIG)
DEALLOCATE(REIGQ)
DEALLOCATE(IEIGQ)
DEALLOCATE(REIGA)
DEALLOCATE(IEIGA)
DEALLOCATE(VA)
DEALLOCATE(LAMBDA)
DEALLOCATE(LAMBDAQ)
DEALLOCATE(EIGA)
DEALLOCATE(W)
DEALLOCATE(AU)
DEALLOCATE(S)
!............................................................
! Generate random M-by-M matrix A. Use DLATMR from
END DO ! LLOOP
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*) ' Test summary for DGEDMD :'
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*)
IF ( NFAIL_Z_XV == 0 ) THEN
WRITE(*,*) '>>>> Z - U*V test PASSED.'
ELSE
WRITE(*,*) 'Z - U*V test FAILED ', NFAIL_Z_XV, ' time(s)'
WRITE(*,*) 'Max error ||Z-U*V||_F was ', TMP_ZXW
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_Z_XV
END IF
IF ( NFAIL_AU == 0 ) THEN
WRITE(*,*) '>>>> A*U test PASSED. '
ELSE
WRITE(*,*) 'A*U test FAILED ', NFAIL_AU, ' time(s)'
WRITE(*,*) 'Max A*U test adjusted error measure was ', TMP_AU
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_AU
END IF
IF ( NFAIL_REZ == 0 ) THEN
WRITE(*,*) '>>>> Rezidual computation test PASSED.'
ELSE
WRITE(*,*) 'Rezidual computation test FAILED ', NFAIL_REZ, 'time(s)'
WRITE(*,*) 'Max residual computing test adjusted error measure was ', TMP_REZ
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_REZ
END IF
IF ( NFAIL_TOTAL == 0 ) THEN
WRITE(*,*) '>>>> DGEDMD :: ALL TESTS PASSED.'
ELSE
WRITE(*,*) NFAIL_TOTAL, 'FAILURES!'
WRITE(*,*) '>>>>>>>>>>>>>> DGEDMD :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
END IF
IF ( TEST_QRDMD ) THEN
WRITE(*,*)
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*) ' Test summary for DGEDMDQ :'
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*)
IF ( NFAIL_SVDIFF == 0 ) THEN
WRITE(*,*) '>>>> DGEDMD and DGEDMDQ computed singular &
&values test PASSED.'
ELSE
WRITE(*,*) 'DGEDMD and DGEDMDQ discrepancies in &
&the singular values unacceptable ', &
NFAIL_SVDIFF, ' times. Test FAILED.'
WRITE(*,*) 'The maximal discrepancy in the singular values (relative to the norm) was ', SVDIFF
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_SVDIFF
END IF
IF ( NFAIL_F_QR == 0 ) THEN
WRITE(*,*) '>>>> F - Q*R test PASSED.'
ELSE
WRITE(*,*) 'F - Q*R test FAILED ', NFAIL_F_QR, ' time(s)'
WRITE(*,*) 'The largest relative residual was ', TMP_FQR
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_F_QR
END IF
IF ( NFAIL_REZQ == 0 ) THEN
WRITE(*,*) '>>>> Rezidual computation test PASSED.'
ELSE
WRITE(*,*) 'Rezidual computation test FAILED ', NFAIL_REZQ, 'time(s)'
WRITE(*,*) 'Max residual computing test adjusted error measure was ', TMP_REZQ
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_REZQ
END IF
IF ( NFAILQ_TOTAL == 0 ) THEN
WRITE(*,*) '>>>>>>> DGEDMDQ :: ALL TESTS PASSED.'
ELSE
WRITE(*,*) NFAILQ_TOTAL, 'FAILURES!'
WRITE(*,*) '>>>>>>> DGEDMDQ :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
END IF
END IF
WRITE(*,*)
WRITE(*,*) 'Test completed.'
STOP
END

+ 792
- 0
lapack-netlib/TESTING/EIG/schkdmd.f90 View File

@@ -0,0 +1,792 @@
! This is a test program for checking the implementations of
! the implementations of the following subroutines
!
! SGEDMD for computation of the
! Dynamic Mode Decomposition (DMD)
! SGEDMDQ for computation of a
! QR factorization based compressed DMD
!
! Developed and supported by:
! ===========================
! Developed and coded by Zlatko Drmac, Faculty of Science,
! University of Zagreb; drmac@math.hr
! In cooperation with
! AIMdyn Inc., Santa Barbara, CA.
! ========================================================
! How to run the code (compiler, link info)
! ========================================================
! Compile as FORTRAN 90 (or later) and link with BLAS and
! LAPACK libraries.
! NOTE: The code is developed and tested on top of the
! Intel MKL library (versions 2022.0.3 and 2022.2.0),
! using the Intel Fortran compiler.
!
! For developers of the C++ implementation
! ========================================================
! See the LAPACK++ and Template Numerical Toolkit (TNT)
!
! Note on a development of the GPU HP implementation
! ========================================================
! Work in progress. See CUDA, MAGMA, SLATE.
! NOTE: The four SVD subroutines used in this code are
! included as a part of R&D and for the completeness.
! This was also an opportunity to test those SVD codes.
! If the scaling option is used all four are essentially
! equally good. For implementations on HP platforms,
! one can use whichever SVD is available.
!... .........................................................
! NOTE:
! When using the Intel MKL 2022.0.3 the subroutine xGESVDQ
! (optionally used in xGEDMD) may cause access violation
! error for x = S, D, C, Z, but only if called with the
! work space query. (At least in our Windows 10 MSVS 2019.)
! The problem can be mitigated by downloading the source
! code of xGESVDQ from the LAPACK repository and use it
! localy instead of the one in the MKL. This seems to
! indicate that the problem is indeed in the MKL.
! This problem did not appear whith Intel MKL 2022.2.0.
!
! NOTE:
! xGESDD seems to have a problem with workspace. In some
! cases the length of the optimal workspace is returned
! smaller than the minimal workspace, as specified in the
! code. As a precaution, all optimal workspaces are
! set as MAX(minimal, optimal).
! Latest implementations of complex xGESDD have different
! length of the real worksapce. We use max value over
! two versions.
!............................................................
!............................................................
!
PROGRAM DMD_TEST
use iso_fortran_env, only: real32
IMPLICIT NONE
integer, parameter :: WP = real32
!............................................................
REAL(KIND=WP), PARAMETER :: ONE = 1.0_WP
REAL(KIND=WP), PARAMETER :: ZERO = 0.0_WP
!............................................................
REAL(KIND=WP), ALLOCATABLE, DIMENSION(:,:) :: &
A, AC, EIGA, LAMBDA, LAMBDAQ, F, F1, F2,&
Z, Z1, S, AU, W, VA, X, X0, Y, Y0, Y1
REAL(KIND=WP), ALLOCATABLE, DIMENSION(:) :: &
DA, DL, DR, REIG, REIGA, REIGQ, IEIG, &
IEIGA, IEIGQ, RES, RES1, RESEX, SINGVX,&
SINGVQX, WORK
INTEGER , ALLOCATABLE, DIMENSION(:) :: IWORK
REAL(KIND=WP) :: AB(2,2), WDUMMY(2)
INTEGER :: IDUMMY(2), ISEED(4), RJOBDATA(8)
REAL(KIND=WP) :: ANORM, COND, CONDL, CONDR, DMAX, EPS, &
TOL, TOL2, SVDIFF, TMP, TMP_AU, &
TMP_FQR, TMP_REZ, TMP_REZQ, TMP_ZXW, &
TMP_EX, XNORM, YNORM
!............................................................
INTEGER :: K, KQ, LDF, LDS, LDA, LDAU, LDW, LDX, LDY, &
LDZ, LIWORK, LWORK, M, N, L, LLOOP, NRNK
INTEGER :: i, iJOBREF, iJOBZ, iSCALE, INFO, KDIFF, &
NFAIL, NFAIL_AU, NFAIL_F_QR, NFAIL_REZ, &
NFAIL_REZQ, NFAIL_SVDIFF, NFAIL_TOTAL, NFAILQ_TOTAL, &
NFAIL_Z_XV, MODE, MODEL, MODER, WHTSVD
INTEGER iNRNK, iWHTSVD, K_TRAJ, LWMINOPT
CHARACTER(LEN=1) GRADE, JOBREF, JOBZ, PIVTNG, RSIGN, &
SCALE, RESIDS, WANTQ, WANTR
LOGICAL TEST_QRDMD
!..... external subroutines (BLAS and LAPACK)
EXTERNAL SAXPY, SGEEV, SGEMM, SGEMV, SLACPY, SLASCL
EXTERNAL SLARNV, SLATMR
!.....external subroutines DMD package, part 1
! subroutines under test
EXTERNAL SGEDMD, SGEDMDQ
!..... external functions (BLAS and LAPACK)
EXTERNAL SLAMCH, SLANGE, SNRM2
REAL(KIND=WP) :: SLAMCH, SLANGE, SNRM2
EXTERNAL LSAME
LOGICAL LSAME
INTRINSIC ABS, INT, MIN, MAX
!............................................................
! The test is always in pairs : ( SGEDMD and SGEDMDQ )
! because the test includes comparing the results (in pairs).
!.....................................................................................
TEST_QRDMD = .TRUE. ! This code by default performs tests on SGEDMDQ
! Since the QR factorizations based algorithm is designed for
! single trajectory data, only single trajectory tests will
! be performed with xGEDMDQ.
WANTQ = 'Q'
WANTR = 'R'
!.................................................................................
EPS = SLAMCH( 'P' ) ! machine precision SP
! Global counters of failures of some particular tests
NFAIL = 0
NFAIL_REZ = 0
NFAIL_REZQ = 0
NFAIL_Z_XV = 0
NFAIL_F_QR = 0
NFAIL_AU = 0
KDIFF = 0
NFAIL_SVDIFF = 0
NFAIL_TOTAL = 0
NFAILQ_TOTAL = 0
DO LLOOP = 1, 4
WRITE(*,*) 'L Loop Index = ', LLOOP
! Set the dimensions of the problem ...
WRITE(*,*) 'M = '
READ(*,*) M
WRITE(*,*) M
! ... and the number of snapshots.
WRITE(*,*) 'N = '
READ(*,*) N
WRITE(*,*) N
! ... Test the dimensions
IF ( ( MIN(M,N) == 0 ) .OR. ( M < N ) ) THEN
WRITE(*,*) 'Bad dimensions. Required: M >= N > 0.'
STOP
END IF
!.............
! The seed inside the LLOOP so that each pass can be reproduced easily.
ISEED(1) = 4
ISEED(2) = 3
ISEED(3) = 2
ISEED(4) = 1
LDA = M
LDF = M
LDX = MAX(M,N+1)
LDY = MAX(M,N+1)
LDW = N
LDZ = M
LDAU = MAX(M,N+1)
LDS = N
TMP_ZXW = ZERO
TMP_AU = ZERO
TMP_REZ = ZERO
TMP_REZQ = ZERO
SVDIFF = ZERO
TMP_EX = ZERO
!
! Test the subroutines on real data snapshots. All
! computation is done in real arithmetic, even when
! Koopman eigenvalues and modes are real.
!
! Allocate memory space
ALLOCATE( A(LDA,M) )
ALLOCATE( AC(LDA,M) )
ALLOCATE( DA(M) )
ALLOCATE( DL(M) )
ALLOCATE( F(LDF,N+1) )
ALLOCATE( F1(LDF,N+1) )
ALLOCATE( F2(LDF,N+1) )
ALLOCATE( X(LDX,N) )
ALLOCATE( X0(LDX,N) )
ALLOCATE( SINGVX(N) )
ALLOCATE( SINGVQX(N) )
ALLOCATE( Y(LDY,N+1) )
ALLOCATE( Y0(LDY,N+1) )
ALLOCATE( Y1(M,N+1) )
ALLOCATE( Z(LDZ,N) )
ALLOCATE( Z1(LDZ,N) )
ALLOCATE( RES(N) )
ALLOCATE( RES1(N) )
ALLOCATE( RESEX(N) )
ALLOCATE( REIG(N) )
ALLOCATE( IEIG(N) )
ALLOCATE( REIGQ(N) )
ALLOCATE( IEIGQ(N) )
ALLOCATE( REIGA(M) )
ALLOCATE( IEIGA(M) )
ALLOCATE( VA(LDA,M) )
ALLOCATE( LAMBDA(N,2) )
ALLOCATE( LAMBDAQ(N,2) )
ALLOCATE( EIGA(M,2) )
ALLOCATE( W(LDW,N) )
ALLOCATE( AU(LDAU,N) )
ALLOCATE( S(N,N) )
TOL = M*EPS
! This mimics O(M*N)*EPS bound for accumulated roundoff error.
! The factor 10 is somewhat arbitrary.
TOL2 = 10*M*N*EPS
!.............
DO K_TRAJ = 1, 2
! Number of intial conditions in the simulation/trajectories (1 or 2)
COND = 1.0D8
DMAX = 1.0D2
RSIGN = 'F'
GRADE = 'N'
MODEL = 6
CONDL = 1.0D2
MODER = 6
CONDR = 1.0D2
PIVTNG = 'N'
! Loop over all parameter MODE values for ZLATMR (+1,..,+6)
DO MODE = 1, 6
ALLOCATE( IWORK(2*M) )
ALLOCATE(DR(N))
CALL SLATMR( M, M, 'S', ISEED, 'N', DA, MODE, COND, &
DMAX, RSIGN, GRADE, DL, MODEL, CONDL, &
DR, MODER, CONDR, PIVTNG, IWORK, M, M, &
ZERO, -ONE, 'N', A, LDA, IWORK(M+1), INFO )
DEALLOCATE(IWORK)
DEALLOCATE(DR)
LWORK = 4*M+1
ALLOCATE(WORK(LWORK))
AC = A
CALL SGEEV( 'N','V', M, AC, M, REIGA, IEIGA, VA, M, &
VA, M, WORK, LWORK, INFO ) ! LAPACK CALL
DEALLOCATE(WORK)
TMP = ZERO
DO i = 1, M
EIGA(i,1) = REIGA(i)
EIGA(i,2) = IEIGA(i)
TMP = MAX( TMP, SQRT(REIGA(i)**2+IEIGA(i)**2))
END DO
! Scale A to have the desirable spectral radius.
CALL SLASCL( 'G', 0, 0, TMP, ONE, M, M, A, M, INFO )
CALL SLASCL( 'G', 0, 0, TMP, ONE, M, 2, EIGA, M, INFO )
! Compute the norm of A
ANORM = SLANGE( 'F', N, N, A, M, WDUMMY )
IF ( K_TRAJ == 2 ) THEN
! generate data with two inital conditions
CALL SLARNV(2, ISEED, M, F1(1,1) )
F1(1:M,1) = 1.0E-10*F1(1:M,1)
DO i = 1, N/2
CALL SGEMV( 'N', M, M, ONE, A, M, F1(1,i), 1, ZERO, &
F1(1,i+1), 1 )
END DO
X0(1:M,1:N/2) = F1(1:M,1:N/2)
Y0(1:M,1:N/2) = F1(1:M,2:N/2+1)
CALL SLARNV(2, ISEED, M, F1(1,1) )
DO i = 1, N-N/2
CALL SGEMV( 'N', M, M, ONE, A, M, F1(1,i), 1, ZERO, &
F1(1,i+1), 1 )
END DO
X0(1:M,N/2+1:N) = F1(1:M,1:N-N/2)
Y0(1:M,N/2+1:N) = F1(1:M,2:N-N/2+1)
ELSE
! single trajectory
CALL SLARNV(2, ISEED, M, F(1,1) )
DO i = 1, N
CALL SGEMV( 'N', M, M, ONE, A, M, F(1,i), 1, ZERO, &
F(1,i+1), 1 )
END DO
X0(1:M,1:N) = F(1:M,1:N)
Y0(1:M,1:N) = F(1:M,2:N+1)
END IF
XNORM = SLANGE( 'F', M, N, X0, LDX, WDUMMY )
YNORM = SLANGE( 'F', M, N, Y0, LDX, WDUMMY )
!............................................................
DO iJOBZ = 1, 4
SELECT CASE ( iJOBZ )
CASE(1)
JOBZ = 'V' ! Ritz vectors will be computed
RESIDS = 'R' ! Residuals will be computed
CASE(2)
JOBZ = 'V'
RESIDS = 'N'
CASE(3)
JOBZ = 'F' ! Ritz vectors in factored form
RESIDS = 'N'
CASE(4)
JOBZ = 'N'
RESIDS = 'N'
END SELECT
DO iJOBREF = 1, 3
SELECT CASE ( iJOBREF )
CASE(1)
JOBREF = 'R' ! Data for refined Ritz vectors
CASE(2)
JOBREF = 'E' ! Exact DMD vectors
CASE(3)
JOBREF = 'N'
END SELECT
DO iSCALE = 1, 4
SELECT CASE ( iSCALE )
CASE(1)
SCALE = 'S' ! X data normalized
CASE(2)
SCALE = 'C' ! X normalized, consist. check
CASE(3)
SCALE = 'Y' ! Y data normalized
CASE(4)
SCALE = 'N'
END SELECT
DO iNRNK = -1, -2, -1
! Two truncation strategies. The "-2" case for R&D
! purposes only - it uses possibly low accuracy small
! singular values, in which case the formulas used in
! the DMD are highly sensitive.
NRNK = iNRNK
DO iWHTSVD = 1, 4
! Check all four options to compute the POD basis
! via the SVD.
WHTSVD = iWHTSVD
DO LWMINOPT = 1, 2
! Workspace query for the minimal (1) and for the optimal
! (2) workspace lengths determined by workspace query.
X(1:M,1:N) = X0(1:M,1:N)
Y(1:M,1:N) = Y0(1:M,1:N)
! SGEDMD: Workspace query and workspace allocation
CALL SGEDMD( SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, M, &
N, X, LDX, Y, LDY, NRNK, TOL, K, REIG, IEIG, Z, &
LDZ, RES, AU, LDAU, W, LDW, S, LDS, WDUMMY, -1, &
IDUMMY, -1, INFO )
LIWORK = IDUMMY(1)
ALLOCATE( IWORK(LIWORK) )
LWORK = INT(WDUMMY(LWMINOPT))
ALLOCATE( WORK(LWORK) )
! SGEDMD test: CALL SGEDMD
CALL SGEDMD( SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, M, &
N, X, LDX, Y, LDY, NRNK, TOL, K, REIG, IEIG, Z, &
LDZ, RES, AU, LDAU, W, LDW, S, LDS, WORK, LWORK,&
IWORK, LIWORK, INFO )
SINGVX(1:N) = WORK(1:N)
!...... SGEDMD check point
IF ( LSAME(JOBZ,'V') ) THEN
! Check that Z = X*W, on return from SGEDMD
! This checks that the returned aigenvectors in Z are
! the product of the SVD'POD basis returned in X
! and the eigenvectors of the rayleigh quotient
! returned in W
CALL SGEMM( 'N', 'N', M, K, K, ONE, X, LDX, W, LDW, &
ZERO, Z1, LDZ )
TMP = ZERO
DO i = 1, K
CALL SAXPY( M, -ONE, Z(1,i), 1, Z1(1,i), 1)
TMP = MAX(TMP, SNRM2( M, Z1(1,i), 1 ) )
END DO
TMP_ZXW = MAX(TMP_ZXW, TMP )
IF ( TMP_ZXW > 10*M*EPS ) THEN
NFAIL_Z_XV = NFAIL_Z_XV + 1
END IF
END IF
!...... SGEDMD check point
IF ( LSAME(JOBREF,'R') ) THEN
! The matrix A*U is returned for computing refined Ritz vectors.
! Check that A*U is computed correctly using the formula
! A*U = Y * V * inv(SIGMA). This depends on the
! accuracy in the computed singular values and vectors of X.
! See the paper for an error analysis.
! Note that the left singular vectors of the input matrix X
! are returned in the array X.
CALL SGEMM( 'N', 'N', M, K, M, ONE, A, LDA, X, LDX, &
ZERO, Z1, LDZ )
TMP = ZERO
DO i = 1, K
CALL SAXPY( M, -ONE, AU(1,i), 1, Z1(1,i), 1)
TMP = MAX( TMP, SNRM2( M, Z1(1,i),1 ) * &
SINGVX(K)/(ANORM*SINGVX(1)) )
END DO
TMP_AU = MAX( TMP_AU, TMP )
IF ( TMP > TOL2 ) THEN
NFAIL_AU = NFAIL_AU + 1
END IF
ELSEIF ( LSAME(JOBREF,'E') ) THEN
! The unscaled vectors of the Exact DMD are computed.
! This option is included for the sake of completeness,
! for users who prefer the Exact DMD vectors. The
! returned vectors are in the real form, in the same way
! as the Ritz vectors. Here we just save the vectors
! and test them separately using a Matlab script.
CALL SGEMM( 'N', 'N', M, K, M, ONE, A, LDA, AU, LDAU, ZERO, Y1, M )
i=1
DO WHILE ( i <= K )
IF ( IEIG(i) == ZERO ) THEN
! have a real eigenvalue with real eigenvector
CALL SAXPY( M, -REIG(i), AU(1,i), 1, Y1(1,i), 1 )
RESEX(i) = SNRM2( M, Y1(1,i), 1) / SNRM2(M,AU(1,i),1)
i = i + 1
ELSE
! Have a complex conjugate pair
! REIG(i) +- sqrt(-1)*IMEIG(i).
! Since all computation is done in real
! arithmetic, the formula for the residual
! is recast for real representation of the
! complex conjugate eigenpair. See the
! description of RES.
AB(1,1) = REIG(i)
AB(2,1) = -IEIG(i)
AB(1,2) = IEIG(i)
AB(2,2) = REIG(i)
CALL SGEMM( 'N', 'N', M, 2, 2, -ONE, AU(1,i), &
M, AB, 2, ONE, Y1(1,i), M )
RESEX(i) = SLANGE( 'F', M, 2, Y1(1,i), M, &
WORK )/ SLANGE( 'F', M, 2, AU(1,i), M, &
WORK )
RESEX(i+1) = RESEX(i)
i = i + 2
END IF
END DO
END IF
!...... SGEDMD check point
IF ( LSAME(RESIDS, 'R') ) THEN
! Compare the residuals returned by SGEDMD with the
! explicitly computed residuals using the matrix A.
! Compute explicitly Y1 = A*Z
CALL SGEMM( 'N', 'N', M, K, M, ONE, A, LDA, Z, LDZ, ZERO, Y1, M )
! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
! of the invariant subspaces that correspond to complex conjugate
! pairs of eigencalues. (See the description of Z in SGEDMD,)
i = 1
DO WHILE ( i <= K )
IF ( IEIG(i) == ZERO ) THEN
! have a real eigenvalue with real eigenvector
CALL SAXPY( M, -REIG(i), Z(1,i), 1, Y1(1,i), 1 )
RES1(i) = SNRM2( M, Y1(1,i), 1)
i = i + 1
ELSE
! Have a complex conjugate pair
! REIG(i) +- sqrt(-1)*IMEIG(i).
! Since all computation is done in real
! arithmetic, the formula for the residual
! is recast for real representation of the
! complex conjugate eigenpair. See the
! description of RES.
AB(1,1) = REIG(i)
AB(2,1) = -IEIG(i)
AB(1,2) = IEIG(i)
AB(2,2) = REIG(i)
CALL SGEMM( 'N', 'N', M, 2, 2, -ONE, Z(1,i), &
M, AB, 2, ONE, Y1(1,i), M )
RES1(i) = SLANGE( 'F', M, 2, Y1(1,i), M, &
WORK )
RES1(i+1) = RES1(i)
i = i + 2
END IF
END DO
TMP = ZERO
DO i = 1, K
TMP = MAX( TMP, ABS(RES(i) - RES1(i)) * &
SINGVX(K)/(ANORM*SINGVX(1)) )
END DO
TMP_REZ = MAX( TMP_REZ, TMP )
IF ( TMP > TOL2 ) THEN
NFAIL_REZ = NFAIL_REZ + 1
END IF
IF ( LSAME(JOBREF,'E') ) THEN
TMP = ZERO
DO i = 1, K
TMP = MAX( TMP, ABS(RES1(i) - RESEX(i))/(RES1(i)+RESEX(i)) )
END DO
TMP_EX = MAX(TMP_EX,TMP)
END IF
END IF
! ... store the results for inspection
DO i = 1, K
LAMBDA(i,1) = REIG(i)
LAMBDA(i,2) = IEIG(i)
END DO
DEALLOCATE(IWORK)
DEALLOCATE(WORK)
!======================================================================
! Now test the SGEDMDQ, if requested.
!======================================================================
IF ( TEST_QRDMD .AND. (K_TRAJ == 1) ) THEN
RJOBDATA(2) = 1
F1 = F
! SGEDMDQ test: Workspace query and workspace allocation
CALL SGEDMDQ( SCALE, JOBZ, RESIDS, WANTQ, WANTR, &
JOBREF, WHTSVD, M, N+1, F1, LDF, X, LDX, Y, &
LDY, NRNK, TOL, KQ, REIGQ, IEIGQ, Z, LDZ, &
RES, AU, LDAU, W, LDW, S, LDS, WDUMMY, &
-1, IDUMMY, -1, INFO )
LIWORK = IDUMMY(1)
ALLOCATE( IWORK(LIWORK) )
LWORK = INT(WDUMMY(LWMINOPT))
ALLOCATE(WORK(LWORK))
! SGEDMDQ test: CALL SGEDMDQ
CALL SGEDMDQ( SCALE, JOBZ, RESIDS, WANTQ, WANTR, &
JOBREF, WHTSVD, M, N+1, F1, LDF, X, LDX, Y, &
LDY, NRNK, TOL, KQ, REIGQ, IEIGQ, Z, LDZ, &
RES, AU, LDAU, W, LDW, S, LDS, &
WORK, LWORK, IWORK, LIWORK, INFO )
SINGVQX(1:KQ) = WORK(MIN(M,N+1)+1: MIN(M,N+1)+KQ)
!..... SGEDMDQ check point
IF ( KQ /= K ) THEN
KDIFF = KDIFF+1
END IF
TMP = ZERO
DO i = 1, MIN(K, KQ)
TMP = MAX(TMP, ABS(SINGVX(i)-SINGVQX(i)) / &
SINGVX(1) )
END DO
SVDIFF = MAX( SVDIFF, TMP )
IF ( TMP > M*N*EPS ) THEN
NFAIL_SVDIFF = NFAIL_SVDIFF + 1
END IF
!..... SGEDMDQ check point
IF ( LSAME(WANTQ,'Q') .AND. LSAME(WANTR,'R') ) THEN
! Check that the QR factors are computed and returned
! as requested. The residual ||F-Q*R||_F / ||F||_F
! is compared to M*N*EPS.
F2 = F
CALL SGEMM( 'N', 'N', M, N+1, MIN(M,N+1), -ONE, F1, &
LDF, Y, LDY, ONE, F2, LDF )
TMP_FQR = SLANGE( 'F', M, N+1, F2, LDF, WORK ) / &
SLANGE( 'F', M, N+1, F, LDF, WORK )
IF ( TMP_FQR > TOL2 ) THEN
NFAIL_F_QR = NFAIL_F_QR + 1
END IF
END IF
!..... SGEDMDQ checkpoint
IF ( LSAME(RESIDS, 'R') ) THEN
! Compare the residuals returned by SGEDMDQ with the
! explicitly computed residuals using the matrix A.
! Compute explicitly Y1 = A*Z
CALL SGEMM( 'N', 'N', M, KQ, M, ONE, A, M, Z, M, ZERO, Y1, M )
! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
! of the invariant subspaces that correspond to complex conjugate
! pairs of eigencalues. (See the description of Z in SGEDMDQ)
i = 1
DO WHILE ( i <= KQ )
IF ( IEIGQ(i) == ZERO ) THEN
! have a real eigenvalue with real eigenvector
CALL SAXPY( M, -REIGQ(i), Z(1,i), 1, Y1(1,i), 1 )
! Y(1:M,i) = Y(1:M,i) - REIG(i)*Z(1:M,i)
RES1(i) = SNRM2( M, Y1(1,i), 1)
i = i + 1
ELSE
! Have a complex conjugate pair
! REIG(i) +- sqrt(-1)*IMEIG(i).
! Since all computation is done in real
! arithmetic, the formula for the residual
! is recast for real representation of the
! complex conjugate eigenpair. See the
! description of RES.
AB(1,1) = REIGQ(i)
AB(2,1) = -IEIGQ(i)
AB(1,2) = IEIGQ(i)
AB(2,2) = REIGQ(i)
CALL SGEMM( 'N', 'N', M, 2, 2, -ONE, Z(1,i), &
M, AB, 2, ONE, Y1(1,i), M ) ! BLAS CALL
! Y(1:M,i:i+1) = Y(1:M,i:i+1) - Z(1:M,i:i+1) * AB ! INTRINSIC
RES1(i) = SLANGE( 'F', M, 2, Y1(1,i), M, &
WORK ) ! LAPACK CALL
RES1(i+1) = RES1(i)
i = i + 2
END IF
END DO
TMP = ZERO
DO i = 1, KQ
TMP = MAX( TMP, ABS(RES(i) - RES1(i)) * &
SINGVQX(K)/(ANORM*SINGVQX(1)) )
END DO
TMP_REZQ = MAX( TMP_REZQ, TMP )
IF ( TMP > TOL2 ) THEN
NFAIL_REZQ = NFAIL_REZQ + 1
END IF
END IF
DO i = 1, KQ
LAMBDAQ(i,1) = REIGQ(i)
LAMBDAQ(i,2) = IEIGQ(i)
END DO
DEALLOCATE(WORK)
DEALLOCATE(IWORK)
END IF ! TEST_QRDMD
!======================================================================
END DO ! LWMINOPT
!write(*,*) 'LWMINOPT loop completed'
END DO ! WHTSVD LOOP
!write(*,*) 'WHTSVD loop completed'
END DO ! NRNK LOOP
!write(*,*) 'NRNK loop completed'
END DO ! SCALE LOOP
!write(*,*) 'SCALE loop completed'
END DO ! JOBF LOOP
!write(*,*) 'JOBREF loop completed'
END DO ! JOBZ LOOP
!write(*,*) 'JOBZ loop completed'
END DO ! MODE -6:6
!write(*,*) 'MODE loop completed'
END DO ! 1 or 2 trajectories
!write(*,*) 'trajectories loop completed'
DEALLOCATE(A)
DEALLOCATE(AC)
DEALLOCATE(DA)
DEALLOCATE(DL)
DEALLOCATE(F)
DEALLOCATE(F1)
DEALLOCATE(F2)
DEALLOCATE(X)
DEALLOCATE(X0)
DEALLOCATE(SINGVX)
DEALLOCATE(SINGVQX)
DEALLOCATE(Y)
DEALLOCATE(Y0)
DEALLOCATE(Y1)
DEALLOCATE(Z)
DEALLOCATE(Z1)
DEALLOCATE(RES)
DEALLOCATE(RES1)
DEALLOCATE(RESEX)
DEALLOCATE(REIG)
DEALLOCATE(IEIG)
DEALLOCATE(REIGQ)
DEALLOCATE(IEIGQ)
DEALLOCATE(REIGA)
DEALLOCATE(IEIGA)
DEALLOCATE(VA)
DEALLOCATE(LAMBDA)
DEALLOCATE(LAMBDAQ)
DEALLOCATE(EIGA)
DEALLOCATE(W)
DEALLOCATE(AU)
DEALLOCATE(S)
!............................................................
! Generate random M-by-M matrix A. Use DLATMR from
END DO ! LLOOP
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*) ' Test summary for SGEDMD :'
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*)
IF ( NFAIL_Z_XV == 0 ) THEN
WRITE(*,*) '>>>> Z - U*V test PASSED.'
ELSE
WRITE(*,*) 'Z - U*V test FAILED ', NFAIL_Z_XV, ' time(s)'
WRITE(*,*) 'Max error ||Z-U*V||_F was ', TMP_ZXW
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_Z_XV
END IF
IF ( NFAIL_AU == 0 ) THEN
WRITE(*,*) '>>>> A*U test PASSED. '
ELSE
WRITE(*,*) 'A*U test FAILED ', NFAIL_AU, ' time(s)'
WRITE(*,*) 'Max A*U test adjusted error measure was ', TMP_AU
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_AU
END IF
IF ( NFAIL_REZ == 0 ) THEN
WRITE(*,*) '>>>> Rezidual computation test PASSED.'
ELSE
WRITE(*,*) 'Rezidual computation test FAILED ', NFAIL_REZ, 'time(s)'
WRITE(*,*) 'Max residual computing test adjusted error measure was ', TMP_REZ
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_REZ
END IF
IF ( NFAIL_TOTAL == 0 ) THEN
WRITE(*,*) '>>>> SGEDMD :: ALL TESTS PASSED.'
ELSE
WRITE(*,*) NFAIL_TOTAL, 'FAILURES!'
WRITE(*,*) '>>>>>>>>>>>>>> SGEDMD :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
END IF
IF ( TEST_QRDMD ) THEN
WRITE(*,*)
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*) ' Test summary for SGEDMDQ :'
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*)
IF ( NFAIL_SVDIFF == 0 ) THEN
WRITE(*,*) '>>>> SGEDMD and SGEDMDQ computed singular &
&values test PASSED.'
ELSE
WRITE(*,*) 'SGEDMD and SGEDMDQ discrepancies in &
&the singular values unacceptable ', &
NFAIL_SVDIFF, ' times. Test FAILED.'
WRITE(*,*) 'The maximal discrepancy in the singular values (relative to the norm) was ', SVDIFF
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_SVDIFF
END IF
IF ( NFAIL_F_QR == 0 ) THEN
WRITE(*,*) '>>>> F - Q*R test PASSED.'
ELSE
WRITE(*,*) 'F - Q*R test FAILED ', NFAIL_F_QR, ' time(s)'
WRITE(*,*) 'The largest relative residual was ', TMP_FQR
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_F_QR
END IF
IF ( NFAIL_REZQ == 0 ) THEN
WRITE(*,*) '>>>> Rezidual computation test PASSED.'
ELSE
WRITE(*,*) 'Rezidual computation test FAILED ', NFAIL_REZQ, 'time(s)'
WRITE(*,*) 'Max residual computing test adjusted error measure was ', TMP_REZQ
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_REZQ
END IF
IF ( NFAILQ_TOTAL == 0 ) THEN
WRITE(*,*) '>>>>>>> SGEDMDQ :: ALL TESTS PASSED.'
ELSE
WRITE(*,*) NFAILQ_TOTAL, 'FAILURES!'
WRITE(*,*) '>>>>>>> SGEDMDQ :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
END IF
END IF
WRITE(*,*)
WRITE(*,*) 'Test completed.'
STOP
END

+ 745
- 0
lapack-netlib/TESTING/EIG/zchkdmd.f90 View File

@@ -0,0 +1,745 @@
! This is a test program for checking the implementations of
! the implementations of the following subroutines
!
! ZGEDMD, for computation of the
! Dynamic Mode Decomposition (DMD)
! ZGEDMDQ, for computation of a
! QR factorization based compressed DMD
!
! Developed and supported by:
! ===========================
! Developed and coded by Zlatko Drmac, Faculty of Science,
! University of Zagreb; drmac@math.hr
! In cooperation with
! AIMdyn Inc., Santa Barbara, CA.
! ========================================================
! How to run the code (compiler, link info)
! ========================================================
! Compile as FORTRAN 90 (or later) and link with BLAS and
! LAPACK libraries.
! NOTE: The code is developed and tested on top of the
! Intel MKL library (versions 2022.0.3 and 2022.2.0),
! using the Intel Fortran compiler.
!
! For developers of the C++ implementation
! ========================================================
! See the LAPACK++ and Template Numerical Toolkit (TNT)
!
! Note on a development of the GPU HP implementation
! ========================================================
! Work in progress. See CUDA, MAGMA, SLATE.
! NOTE: The four SVD subroutines used in this code are
! included as a part of R&D and for the completeness.
! This was also an opportunity to test those SVD codes.
! If the scaling option is used all four are essentially
! equally good. For implementations on HP platforms,
! one can use whichever SVD is available.
!............................................................
!............................................................
!............................................................
!
PROGRAM DMD_TEST
use iso_fortran_env, only: real64
IMPLICIT NONE
integer, parameter :: WP = real64
!............................................................
REAL(KIND=WP), PARAMETER :: ONE = 1.0_WP
REAL(KIND=WP), PARAMETER :: ZERO = 0.0_WP
COMPLEX(KIND=WP), PARAMETER :: ZONE = ( 1.0_WP, 0.0_WP )
COMPLEX(KIND=WP), PARAMETER :: ZZERO = ( 0.0_WP, 0.0_WP )
!............................................................
REAL(KIND=WP), ALLOCATABLE, DIMENSION(:) :: RES, &
RES1, RESEX, SINGVX, SINGVQX, WORK
INTEGER , ALLOCATABLE, DIMENSION(:) :: IWORK
REAL(KIND=WP) :: WDUMMY(2)
INTEGER :: IDUMMY(4), ISEED(4)
REAL(KIND=WP) :: ANORM, COND, CONDL, CONDR, EPS, &
TOL, TOL2, SVDIFF, TMP, TMP_AU, &
TMP_FQR, TMP_REZ, TMP_REZQ, TMP_ZXW, &
TMP_EX
!............................................................
COMPLEX(KIND=WP) :: ZMAX
INTEGER :: LZWORK
COMPLEX(KIND=WP), ALLOCATABLE, DIMENSION(:,:) :: ZA, ZAC, &
ZAU, ZF, ZF0, ZF1, ZS, ZW, &
ZX, ZX0, ZY, ZY0, ZY1, ZZ, ZZ1
COMPLEX(KIND=WP), ALLOCATABLE, DIMENSION(:) :: ZDA, ZDR, &
ZDL, ZEIGS, ZEIGSA, ZWORK
COMPLEX(KIND=WP) :: ZDUMMY(22), ZDUM2X2(2,2)
!............................................................
INTEGER :: K, KQ, LDF, LDS, LDA, LDAU, LDW, LDX, LDY, &
LDZ, LIWORK, LWORK, M, N, LLOOP, NRNK, NRNKsp
INTEGER :: i, iJOBREF, iJOBZ, iSCALE, INFO, j, &
NFAIL, NFAIL_AU, NFAIL_F_QR, NFAIL_REZ, &
NFAIL_REZQ, NFAIL_SVDIFF, NFAIL_TOTAL, NFAILQ_TOTAL, &
NFAIL_Z_XV, MODE, MODEL, MODER, WHTSVD, &
WHTSVDsp
INTEGER :: iNRNK, iWHTSVD, K_TRAJ, LWMINOPT
CHARACTER :: GRADE, JOBREF, JOBZ, PIVTNG, RSIGN, &
SCALE, RESIDS, WANTQ, WANTR
LOGICAL :: TEST_QRDMD
!.....external subroutines (BLAS and LAPACK)
EXTERNAL DAXPY, DGEEV, DGEMM, DGEMV, DLACPY, DLASCL
EXTERNAL ZGEEV, ZGEMV, ZLASCL
EXTERNAL ZLARNV, ZLATMR
EXTERNAL ZAXPY, ZGEMM
!.....external subroutines DMD package, part 1
! subroutines under test
EXTERNAL ZGEDMD, ZGEDMDQ
!.....external functions (BLAS and LAPACK)
EXTERNAL DLAMCH, DZNRM2
REAL(KIND=WP) :: DLAMCH, DZNRM2
REAL(KIND=WP) :: ZLANGE
EXTERNAL IZAMAX
INTEGER IZAMAX
EXTERNAL LSAME
LOGICAL LSAME
INTRINSIC ABS, INT, MIN, MAX, SIGN
!............................................................
! The test is always in pairs : ( ZGEDMD and ZGEDMDQ )
! because the test includes comparing the results (in pairs).
!.....................................................................................
TEST_QRDMD = .TRUE. ! This code by default performs tests on ZGEDMDQ
! Since the QR factorizations based algorithm is designed for
! single trajectory data, only single trajectory tests will
! be performed with xGEDMDQ.
WANTQ = 'Q'
WANTR = 'R'
!.................................................................................
EPS = DLAMCH( 'P' ) ! machine precision DP
! Global counters of failures of some particular tests
NFAIL = 0
NFAIL_REZ = 0
NFAIL_REZQ = 0
NFAIL_Z_XV = 0
NFAIL_F_QR = 0
NFAIL_AU = 0
NFAIL_SVDIFF = 0
NFAIL_TOTAL = 0
NFAILQ_TOTAL = 0
DO LLOOP = 1, 4
WRITE(*,*) 'L Loop Index = ', LLOOP
! Set the dimensions of the problem ...
WRITE(*,*) 'M = '
READ(*,*) M
WRITE(*,*) M
! ... and the number of snapshots.
WRITE(*,*) 'N = '
READ(*,*) N
WRITE(*,*) N
! ... Test the dimensions
IF ( ( MIN(M,N) == 0 ) .OR. ( M < N ) ) THEN
WRITE(*,*) 'Bad dimensions. Required: M >= N > 0.'
STOP
END IF
!.............
! The seed inside the LLOOP so that each pass can be reproduced easily.
ISEED(1) = 4
ISEED(2) = 3
ISEED(3) = 2
ISEED(4) = 1
LDA = M
LDF = M
LDX = M
LDY = M
LDW = N
LDZ = M
LDAU = M
LDS = N
TMP_ZXW = ZERO
TMP_AU = ZERO
TMP_REZ = ZERO
TMP_REZQ = ZERO
SVDIFF = ZERO
TMP_EX = ZERO
ALLOCATE( ZA(LDA,M) )
ALLOCATE( ZAC(LDA,M) )
ALLOCATE( ZF(LDF,N+1) )
ALLOCATE( ZF0(LDF,N+1) )
ALLOCATE( ZF1(LDF,N+1) )
ALLOCATE( ZX(LDX,N) )
ALLOCATE( ZX0(LDX,N) )
ALLOCATE( ZY(LDY,N+1) )
ALLOCATE( ZY0(LDY,N+1) )
ALLOCATE( ZY1(LDY,N+1) )
ALLOCATE( ZAU(LDAU,N) )
ALLOCATE( ZW(LDW,N) )
ALLOCATE( ZS(LDS,N) )
ALLOCATE( ZZ(LDZ,N) )
ALLOCATE( ZZ1(LDZ,N) )
ALLOCATE( RES(N) )
ALLOCATE( RES1(N) )
ALLOCATE( RESEX(N) )
ALLOCATE( ZEIGS(N) )
ALLOCATE( SINGVX(N) )
ALLOCATE( SINGVQX(N) )
TOL = M*EPS
! This mimics O(M*N)*EPS bound for accumulated roundoff error.
! The factor 10 is somewhat arbitrary.
TOL2 = 10*M*N*EPS
!.............
DO K_TRAJ = 1, 2
! Number of intial conditions in the simulation/trajectories (1 or 2)
COND = 1.0D4
ZMAX = (1.0D1,1.0D1)
RSIGN = 'F'
GRADE = 'N'
MODEL = 6
CONDL = 1.0D1
MODER = 6
CONDR = 1.0D1
PIVTNG = 'N'
! Loop over all parameter MODE values for ZLATMR (+1,..,+6)
DO MODE = 1, 6
ALLOCATE( IWORK(2*M) )
ALLOCATE( ZDA(M) )
ALLOCATE( ZDL(M) )
ALLOCATE( ZDR(M) )
CALL ZLATMR( M, M, 'N', ISEED, 'N', ZDA, MODE, COND, &
ZMAX, RSIGN, GRADE, ZDL, MODEL, CONDL, &
ZDR, MODER, CONDR, PIVTNG, IWORK, M, M, &
ZERO, -ONE, 'N', ZA, LDA, IWORK(M+1), INFO )
DEALLOCATE( ZDR )
DEALLOCATE( ZDL )
DEALLOCATE( ZDA )
DEALLOCATE( IWORK )
LZWORK = MAX(1,2*M)
ALLOCATE( ZEIGSA(M) )
ALLOCATE( ZWORK(LZWORK) )
ALLOCATE( WORK(2*M) )
ZAC(1:M,1:M) = ZA(1:M,1:M)
CALL ZGEEV( 'N','N', M, ZAC, LDA, ZEIGSA, ZDUM2X2, 2, &
ZDUM2X2, 2, ZWORK, LZWORK, WORK, INFO ) ! LAPACK CALL
DEALLOCATE(WORK)
DEALLOCATE(ZWORK)
TMP = ABS(ZEIGSA(IZAMAX(M, ZEIGSA, 1))) ! The spectral radius of ZA
! Scale the matrix ZA to have unit spectral radius.
CALL ZLASCL( 'G',0, 0, TMP, ONE, M, M, &
ZA, LDA, INFO )
CALL ZLASCL( 'G',0, 0, TMP, ONE, M, 1, &
ZEIGSA, M, INFO )
ANORM = ZLANGE( 'F', M, M, ZA, LDA, WDUMMY )
IF ( K_TRAJ == 2 ) THEN
! generate data as two trajectories
! with two inital conditions
CALL ZLARNV(2, ISEED, M, ZF(1,1) )
DO i = 1, N/2
CALL ZGEMV( 'N', M, M, ZONE, ZA, LDA, ZF(1,i), 1, &
ZZERO, ZF(1,i+1), 1 )
END DO
ZX0(1:M,1:N/2) = ZF(1:M,1:N/2)
ZY0(1:M,1:N/2) = ZF(1:M,2:N/2+1)
CALL ZLARNV(2, ISEED, M, ZF(1,1) )
DO i = 1, N-N/2
CALL ZGEMV( 'N', M, M, ZONE, ZA, LDA, ZF(1,i), 1, &
ZZERO, ZF(1,i+1), 1 )
END DO
ZX0(1:M,N/2+1:N) = ZF(1:M,1:N-N/2)
ZY0(1:M,N/2+1:N) = ZF(1:M,2:N-N/2+1)
ELSE
CALL ZLARNV(2, ISEED, M, ZF(1,1) )
DO i = 1, N
CALL ZGEMV( 'N', M, M, ZONE, ZA, M, ZF(1,i), 1, &
ZZERO, ZF(1,i+1), 1 )
END DO
ZF0(1:M,1:N+1) = ZF(1:M,1:N+1)
ZX0(1:M,1:N) = ZF0(1:M,1:N)
ZY0(1:M,1:N) = ZF0(1:M,2:N+1)
END IF
DEALLOCATE( ZEIGSA )
!........................................................................
DO iJOBZ = 1, 4
SELECT CASE ( iJOBZ )
CASE(1)
JOBZ = 'V'
RESIDS = 'R'
CASE(2)
JOBZ = 'V'
RESIDS = 'N'
CASE(3)
JOBZ = 'F'
RESIDS = 'N'
CASE(4)
JOBZ = 'N'
RESIDS = 'N'
END SELECT
DO iJOBREF = 1, 3
SELECT CASE ( iJOBREF )
CASE(1)
JOBREF = 'R'
CASE(2)
JOBREF = 'E'
CASE(3)
JOBREF = 'N'
END SELECT
DO iSCALE = 1, 4
SELECT CASE ( iSCALE )
CASE(1)
SCALE = 'S'
CASE(2)
SCALE = 'C'
CASE(3)
SCALE = 'Y'
CASE(4)
SCALE = 'N'
END SELECT
DO iNRNK = -1, -2, -1
NRNK = iNRNK
NRNKsp = iNRNK
DO iWHTSVD = 1, 3
! Check all four options to compute the POD basis
! via the SVD.
WHTSVD = iWHTSVD
WHTSVDsp = iWHTSVD
DO LWMINOPT = 1, 2
! Workspace query for the minimal (1) and for the optimal
! (2) workspace lengths determined by workspace query.
! ZGEDMD is always tested and its results are also used for
! comparisons with ZGEDMDQ.
ZX(1:M,1:N) = ZX0(1:M,1:N)
ZY(1:M,1:N) = ZY0(1:M,1:N)
CALL ZGEDMD( SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, ZX, LDX, ZY, LDY, NRNK, TOL, &
K, ZEIGS, ZZ, LDZ, RES, ZAU, LDAU, &
ZW, LDW, ZS, LDS, ZDUMMY, -1, &
WDUMMY, -1, IDUMMY, -1, INFO )
IF ( (INFO .EQ. 2) .OR. ( INFO .EQ. 3 ) &
.OR. ( INFO < 0 ) ) THEN
WRITE(*,*) 'Call to ZGEDMD workspace query failed. &
&Check the calling sequence and the code.'
WRITE(*,*) 'The error code is ', INFO
WRITE(*,*) 'The input parameters were ', &
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL, LDZ, LDAU, LDW, LDS
STOP
END IF
LZWORK = INT(ZDUMMY(LWMINOPT))
LWORK = INT(WDUMMY(1))
LIWORK = IDUMMY(1)
ALLOCATE(ZWORK(LZWORK))
ALLOCATE( WORK(LWORK))
ALLOCATE(IWORK(LIWORK))
CALL ZGEDMD( SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, ZX, LDX, ZY, LDY, NRNK, TOL, &
K, ZEIGS, ZZ, LDZ, RES, ZAU, LDAU, &
ZW, LDW, ZS, LDS, ZWORK, LZWORK, &
WORK, LWORK, IWORK, LIWORK, INFO )
IF ( INFO /= 0 ) THEN
WRITE(*,*) 'Call to ZGEDMD failed. &
&Check the calling sequence and the code.'
WRITE(*,*) 'The error code is ', INFO
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
STOP
END IF
SINGVX(1:N) = WORK(1:N)
!...... ZGEDMD check point
IF ( LSAME(JOBZ,'V') ) THEN
! Check that Z = X*W, on return from ZGEDMD
! This checks that the returned eigenvectors in Z are
! the product of the SVD'POD basis returned in X
! and the eigenvectors of the rayleigh quotient
! returned in W
CALL ZGEMM( 'N', 'N', M, K, K, ZONE, ZX, LDX, ZW, LDW, &
ZZERO, ZZ1, LDZ )
TMP = ZERO
DO i = 1, K
CALL ZAXPY( M, -ZONE, ZZ(1,i), 1, ZZ1(1,i), 1)
TMP = MAX(TMP, DZNRM2( M, ZZ1(1,i), 1 ) )
END DO
TMP_ZXW = MAX(TMP_ZXW, TMP )
IF ( TMP_ZXW <= 10*M*EPS ) THEN
!WRITE(*,*) ' :) .... OK .........ZGEDMD PASSED.'
ELSE
NFAIL_Z_XV = NFAIL_Z_XV + 1
WRITE(*,*) ':( .................ZGEDMD FAILED!', &
'Check the code for implementation errors.'
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
END IF
END IF
!...... ZGEDMD check point
IF ( LSAME(JOBREF,'R') ) THEN
! The matrix A*U is returned for computing refined Ritz vectors.
! Check that A*U is computed correctly using the formula
! A*U = Y * V * inv(SIGMA). This depends on the
! accuracy in the computed singular values and vectors of X.
! See the paper for an error analysis.
! Note that the left singular vectors of the input matrix X
! are returned in the array X.
CALL ZGEMM( 'N', 'N', M, K, M, ZONE, ZA, LDA, ZX, LDX, &
ZZERO, ZZ1, LDZ )
TMP = ZERO
DO i = 1, K
CALL ZAXPY( M, -ZONE, ZAU(1,i), 1, ZZ1(1,i), 1)
TMP = MAX( TMP, DZNRM2( M, ZZ1(1,i),1 ) * &
SINGVX(K)/(ANORM*SINGVX(1)) )
END DO
TMP_AU = MAX( TMP_AU, TMP )
IF ( TMP <= TOL2 ) THEN
!WRITE(*,*) ':) .... OK .........ZGEDMD PASSED.'
ELSE
NFAIL_AU = NFAIL_AU + 1
WRITE(*,*) ':( .................ZGEDMD FAILED!', &
'Check the code for implementation errors.'
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
END IF
ELSEIF ( LSAME(JOBREF,'E') ) THEN
! The unscaled vectors of the Exact DMD are computed.
! This option is included for the sake of completeness,
! for users who prefer the Exact DMD vectors. The
! returned vectors are in the real form, in the same way
! as the Ritz vectors. Here we just save the vectors
! and test them separately using a Matlab script.
CALL ZGEMM( 'N', 'N', M, K, M, ZONE, ZA, LDA, ZAU, LDAU, ZZERO, ZY1, LDY )
DO i=1, K
! have a real eigenvalue with real eigenvector
CALL ZAXPY( M, -ZEIGS(i), ZAU(1,i), 1, ZY1(1,i), 1 )
RESEX(i) = DZNRM2( M, ZY1(1,i), 1) / DZNRM2(M,ZAU(1,i),1)
END DO
END IF
!...... ZGEDMD check point
IF ( LSAME(RESIDS, 'R') ) THEN
! Compare the residuals returned by ZGEDMD with the
! explicitly computed residuals using the matrix A.
! Compute explicitly Y1 = A*Z
CALL ZGEMM( 'N', 'N', M, K, M, ZONE, ZA, LDA, ZZ, LDZ, ZZERO, ZY1, LDY )
! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
! of the invariant subspaces that correspond to complex conjugate
! pairs of eigencalues. (See the description of Z in ZGEDMD,)
DO i=1, K
! have a real eigenvalue with real eigenvector
CALL ZAXPY( M, -ZEIGS(i), ZZ(1,i), 1, ZY1(1,i), 1 )
RES1(i) = DZNRM2( M, ZY1(1,i), 1)
END DO
TMP = ZERO
DO i = 1, K
TMP = MAX( TMP, ABS(RES(i) - RES1(i)) * &
SINGVX(K)/(ANORM*SINGVX(1)) )
END DO
TMP_REZ = MAX( TMP_REZ, TMP )
IF ( TMP <= TOL2 ) THEN
!WRITE(*,*) ':) .... OK ..........ZGEDMD PASSED.'
ELSE
NFAIL_REZ = NFAIL_REZ + 1
WRITE(*,*) ':( ..................ZGEDMD FAILED!', &
'Check the code for implementation errors.'
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, JOBREF, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
END IF
IF ( LSAME(JOBREF,'E') ) THEN
TMP = ZERO
DO i = 1, K
TMP = MAX( TMP, ABS(RES1(i) - RESEX(i))/(RES1(i)+RESEX(i)) )
END DO
TMP_EX = MAX(TMP_EX,TMP)
END IF
END IF
DEALLOCATE(ZWORK)
DEALLOCATE(WORK)
DEALLOCATE(IWORK)
IF ( TEST_QRDMD .AND. (K_TRAJ == 1) ) THEN
ZF(1:M,1:N+1) = ZF0(1:M,1:N+1)
CALL ZGEDMDQ( SCALE, JOBZ, RESIDS, WANTQ, WANTR, JOBREF, &
WHTSVD, M, N+1, ZF, LDF, ZX, LDX, ZY, LDY, &
NRNK, TOL, K, ZEIGS, ZZ, LDZ, RES, ZAU, &
LDAU, ZW, LDW, ZS, LDS, ZDUMMY, -1, &
WDUMMY, -1, IDUMMY, -1, INFO )
LZWORK = INT(ZDUMMY(LWMINOPT))
ALLOCATE( ZWORK(LZWORK) )
LIWORK = IDUMMY(1)
ALLOCATE(IWORK(LIWORK))
LWORK = INT(WDUMMY(1))
ALLOCATE(WORK(LWORK))
CALL ZGEDMDQ( SCALE, JOBZ, RESIDS, WANTQ, WANTR, JOBREF, &
WHTSVD, M, N+1, ZF, LDF, ZX, LDX, ZY, LDY, &
NRNK, TOL, KQ, ZEIGS, ZZ, LDZ, RES, ZAU, &
LDAU, ZW, LDW, ZS, LDS, ZWORK, LZWORK, &
WORK, LWORK, IWORK, LIWORK, INFO )
IF ( INFO /= 0 ) THEN
WRITE(*,*) 'Call to ZGEDMDQ failed. &
&Check the calling sequence and the code.'
WRITE(*,*) 'The error code is ', INFO
WRITE(*,*) 'The input parameters were ',&
SCALE, JOBZ, RESIDS, WANTQ, WANTR, WHTSVD, &
M, N, LDX, LDY, NRNK, TOL
STOP
END IF
SINGVQX(1:N) = WORK(1:N)
!..... ZGEDMDQ check point
IF ( 1 == 0 ) THEN
! Comparison of ZGEDMD and ZGEDMDQ singular values disabled
TMP = ZERO
DO i = 1, MIN(K, KQ)
TMP = MAX(TMP, ABS(SINGVX(i)-SINGVQX(i)) / &
SINGVX(1) )
END DO
SVDIFF = MAX( SVDIFF, TMP )
IF ( TMP > M*N*EPS ) THEN
WRITE(*,*) 'FAILED! Something was wrong with the run.'
NFAIL_SVDIFF = NFAIL_SVDIFF + 1
DO j =1, 3
write(*,*) j, SINGVX(j), SINGVQX(j)
read(*,*)
END DO
END IF
END IF
!..... ZGEDMDQ check point
IF ( LSAME(WANTQ,'Q') .AND. LSAME(WANTR,'R') ) THEN
! Check that the QR factors are computed and returned
! as requested. The residual ||F-Q*R||_F / ||F||_F
! is compared to M*N*EPS.
ZF1(1:M,1:N+1) = ZF0(1:M,1:N+1)
CALL ZGEMM( 'N', 'N', M, N+1, MIN(M,N+1), -ZONE, ZF, &
LDF, ZY, LDY, ZONE, ZF1, LDF )
TMP_FQR = ZLANGE( 'F', M, N+1, ZF1, LDF, WORK ) / &
ZLANGE( 'F', M, N+1, ZF0, LDF, WORK )
IF ( TMP_FQR > TOL2 ) THEN
WRITE(*,*) 'FAILED! Something was wrong with the run.'
NFAIL_F_QR = NFAIL_F_QR + 1
ELSE
!WRITE(*,*) '........ PASSED.'
END IF
END IF
!..... ZGEDMDQ check point
IF ( LSAME(RESIDS, 'R') ) THEN
! Compare the residuals returned by ZGEDMDQ with the
! explicitly computed residuals using the matrix A.
! Compute explicitly Y1 = A*Z
CALL ZGEMM( 'N', 'N', M, KQ, M, ZONE, ZA, LDA, ZZ, LDZ, ZZERO, ZY1, LDY )
! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
! of the invariant subspaces that correspond to complex conjugate
! pairs of eigencalues. (See the description of Z in ZGEDMDQ)
DO i=1, KQ
! have a real eigenvalue with real eigenvector
CALL ZAXPY( M, -ZEIGS(i), ZZ(1,i), 1, ZY1(1,i), 1 )
! Y(1:M,i) = Y(1:M,i) - REIG(i)*Z(1:M,i)
RES1(i) = DZNRM2( M, ZY1(1,i), 1)
END DO
TMP = ZERO
DO i = 1, KQ
TMP = MAX( TMP, ABS(RES(i) - RES1(i)) * &
SINGVQX(KQ)/(ANORM*SINGVQX(1)) )
END DO
TMP_REZQ = MAX( TMP_REZQ, TMP )
IF ( TMP <= TOL2 ) THEN
!WRITE(*,*) '.... OK ........ ZGEDMDQ PASSED.'
ELSE
NFAIL_REZQ = NFAIL_REZQ + 1
WRITE(*,*) '................ ZGEDMDQ FAILED!', &
'Check the code for implementation errors.'
STOP
END IF
END IF
DEALLOCATE( ZWORK )
DEALLOCATE( WORK )
DEALLOCATE( IWORK )
END IF ! ZGEDMDQ
!.......................................................................................................
END DO ! LWMINOPT
!write(*,*) 'LWMINOPT loop completed'
END DO ! iWHTSVD
!write(*,*) 'WHTSVD loop completed'
END DO ! iNRNK -2:-1
!write(*,*) 'NRNK loop completed'
END DO ! iSCALE 1:4
!write(*,*) 'SCALE loop completed'
END DO
!write(*,*) 'JOBREF loop completed'
END DO ! iJOBZ
!write(*,*) 'JOBZ loop completed'
END DO ! MODE -6:6
!write(*,*) 'MODE loop completed'
END DO ! 1 or 2 trajectories
!write(*,*) 'trajectories loop completed'
DEALLOCATE( ZA )
DEALLOCATE( ZAC )
DEALLOCATE( ZZ )
DEALLOCATE( ZF )
DEALLOCATE( ZF0 )
DEALLOCATE( ZF1 )
DEALLOCATE( ZX )
DEALLOCATE( ZX0 )
DEALLOCATE( ZY )
DEALLOCATE( ZY0 )
DEALLOCATE( ZY1 )
DEALLOCATE( ZAU )
DEALLOCATE( ZW )
DEALLOCATE( ZS )
DEALLOCATE( ZZ1 )
DEALLOCATE( RES )
DEALLOCATE( RES1 )
DEALLOCATE( RESEX )
DEALLOCATE( ZEIGS )
DEALLOCATE( SINGVX )
DEALLOCATE( SINGVQX )
END DO ! LLOOP
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*) ' Test summary for ZGEDMD :'
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*)
IF ( NFAIL_Z_XV == 0 ) THEN
WRITE(*,*) '>>>> Z - U*V test PASSED.'
ELSE
WRITE(*,*) 'Z - U*V test FAILED ', NFAIL_Z_XV, ' time(s)'
WRITE(*,*) 'Max error ||Z-U*V||_F was ', TMP_ZXW
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_Z_XV
END IF
IF ( NFAIL_AU == 0 ) THEN
WRITE(*,*) '>>>> A*U test PASSED. '
ELSE
WRITE(*,*) 'A*U test FAILED ', NFAIL_AU, ' time(s)'
WRITE(*,*) 'Max A*U test adjusted error measure was ', TMP_AU
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_AU
END IF
IF ( NFAIL_REZ == 0 ) THEN
WRITE(*,*) '>>>> Rezidual computation test PASSED.'
ELSE
WRITE(*,*) 'Rezidual computation test FAILED ', NFAIL_REZ, 'time(s)'
WRITE(*,*) 'Max residual computing test adjusted error measure was ', TMP_REZ
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAIL_TOTAL = NFAIL_TOTAL + NFAIL_REZ
END IF
IF ( NFAIL_TOTAL == 0 ) THEN
WRITE(*,*) '>>>> ZGEDMD :: ALL TESTS PASSED.'
ELSE
WRITE(*,*) NFAIL_TOTAL, 'FAILURES!'
WRITE(*,*) '>>>>>>>>>>>>>> ZGEDMD :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
END IF
IF ( TEST_QRDMD ) THEN
WRITE(*,*)
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*) ' Test summary for ZGEDMDQ :'
WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
WRITE(*,*)
IF ( NFAIL_SVDIFF == 0 ) THEN
WRITE(*,*) '>>>> ZGEDMD and ZGEDMDQ computed singular &
&values test PASSED.'
ELSE
WRITE(*,*) 'ZGEDMD and ZGEDMDQ discrepancies in &
&the singular values unacceptable ', &
NFAIL_SVDIFF, ' times. Test FAILED.'
WRITE(*,*) 'The maximal discrepancy in the singular values (relative to the norm) was ', SVDIFF
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_SVDIFF
END IF
IF ( NFAIL_F_QR == 0 ) THEN
WRITE(*,*) '>>>> F - Q*R test PASSED.'
ELSE
WRITE(*,*) 'F - Q*R test FAILED ', NFAIL_F_QR, ' time(s)'
WRITE(*,*) 'The largest relative residual was ', TMP_FQR
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_F_QR
END IF
IF ( NFAIL_REZQ == 0 ) THEN
WRITE(*,*) '>>>> Rezidual computation test PASSED.'
ELSE
WRITE(*,*) 'Rezidual computation test FAILED ', NFAIL_REZQ, 'time(s)'
WRITE(*,*) 'Max residual computing test adjusted error measure was ', TMP_REZQ
WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', EPS
NFAILQ_TOTAL = NFAILQ_TOTAL + NFAIL_REZQ
END IF
IF ( NFAILQ_TOTAL == 0 ) THEN
WRITE(*,*) '>>>>>>> ZGEDMDQ :: ALL TESTS PASSED.'
ELSE
WRITE(*,*) NFAILQ_TOTAL, 'FAILURES!'
WRITE(*,*) '>>>>>>> ZGEDMDQ :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
END IF
END IF
WRITE(*,*)
WRITE(*,*) 'Test completed.'
STOP
END

+ 39
- 4
lapack-netlib/TESTING/Makefile View File

@@ -61,6 +61,8 @@ SEIGTST= snep.out \
scsd.out \ scsd.out \
slse.out slse.out


SDMDEIGTST= sdmd.out

CEIGTST= cnep.out \ CEIGTST= cnep.out \
csep.out \ csep.out \
cse2.out \ cse2.out \
@@ -82,6 +84,8 @@ CEIGTST= cnep.out \
ccsd.out \ ccsd.out \
clse.out clse.out


CDMDEIGTST= cdmd.out

DEIGTST= dnep.out \ DEIGTST= dnep.out \
dsep.out \ dsep.out \
dse2.out \ dse2.out \
@@ -103,6 +107,8 @@ DEIGTST= dnep.out \
dcsd.out \ dcsd.out \
dlse.out dlse.out


DDMDEIGTST= ddmd.out

ZEIGTST= znep.out \ ZEIGTST= znep.out \
zsep.out \ zsep.out \
zse2.out \ zse2.out \
@@ -124,6 +130,7 @@ ZEIGTST= znep.out \
zcsd.out \ zcsd.out \
zlse.out zlse.out


ZDMDEIGTST= zdmd.out


SLINTST= stest.out SLINTST= stest.out


@@ -142,10 +149,10 @@ ZLINTST= ztest.out
ZLINTSTPROTO= zctest.out ztest_rfp.out ZLINTSTPROTO= zctest.out ztest_rfp.out


.PHONY: single complex double complex16 .PHONY: single complex double complex16
single: $(SLINTST) $(SEIGTST)
complex: $(CLINTST) $(CEIGTST)
double: $(DLINTST) $(DEIGTST)
complex16: $(ZLINTST) $(ZEIGTST)
single: $(SLINTST) $(SEIGTST) $(SDMDEIGTST)
complex: $(CLINTST) $(CEIGTST) $(CDMDEIGTST)
double: $(DLINTST) $(DEIGTST) $(DDMDEIGTST)
complex16: $(ZLINTST) $(ZEIGTST) $(ZDMDEIGTST)


.PHONY: singleproto complexproto doubleproto complex16proto .PHONY: singleproto complexproto doubleproto complex16proto
singleproto: $(SLINTSTPROTO) singleproto: $(SLINTSTPROTO)
@@ -297,6 +304,10 @@ scsd.out: csd.in EIG/xeigtsts
slse.out: lse.in EIG/xeigtsts slse.out: lse.in EIG/xeigtsts
@echo LSE: Testing Constrained Linear Least Squares routines @echo LSE: Testing Constrained Linear Least Squares routines
./EIG/xeigtsts < lse.in > $@ 2>&1 ./EIG/xeigtsts < lse.in > $@ 2>&1

sdmd.out: sdmd.in EIG/xdmdeigtsts
@echo DMD: Testing Dynamic Mode Decomposition routines
./EIG/xdmdeigtsts < sdmd.in > $@ 2>&1
# #
# ======== COMPLEX EIG TESTS =========================== # ======== COMPLEX EIG TESTS ===========================


@@ -379,6 +390,10 @@ ccsd.out: csd.in EIG/xeigtstc
clse.out: lse.in EIG/xeigtstc clse.out: lse.in EIG/xeigtstc
@echo LSE: Testing Constrained Linear Least Squares routines @echo LSE: Testing Constrained Linear Least Squares routines
./EIG/xeigtstc < lse.in > $@ 2>&1 ./EIG/xeigtstc < lse.in > $@ 2>&1

cdmd.out: cdmd.in EIG/xdmdeigtstc
@echo DMD: Testing Dynamic Mode Decomposition routines
./EIG/xdmdeigtstc < cdmd.in > $@ 2>&1
# #
# ======== DOUBLE EIG TESTS =========================== # ======== DOUBLE EIG TESTS ===========================


@@ -461,6 +476,10 @@ dcsd.out: csd.in EIG/xeigtstd
dlse.out: lse.in EIG/xeigtstd dlse.out: lse.in EIG/xeigtstd
@echo LSE: Testing Constrained Linear Least Squares routines @echo LSE: Testing Constrained Linear Least Squares routines
./EIG/xeigtstd < lse.in > $@ 2>&1 ./EIG/xeigtstd < lse.in > $@ 2>&1

ddmd.out: ddmd.in EIG/xdmdeigtstd
@echo DMD: Testing Dynamic Mode Decomposition routines
./EIG/xdmdeigtstd < ddmd.in > $@ 2>&1
# #
# ======== COMPLEX16 EIG TESTS =========================== # ======== COMPLEX16 EIG TESTS ===========================


@@ -543,6 +562,10 @@ zcsd.out: csd.in EIG/xeigtstz
zlse.out: lse.in EIG/xeigtstz zlse.out: lse.in EIG/xeigtstz
@echo LSE: Testing Constrained Linear Least Squares routines @echo LSE: Testing Constrained Linear Least Squares routines
./EIG/xeigtstz < lse.in > $@ 2>&1 ./EIG/xeigtstz < lse.in > $@ 2>&1

zdmd.out: zdmd.in EIG/xdmdeigtstz
@echo DMD: Testing Dynamic Mode Decomposition routines
./EIG/xdmdeigtstz < zdmd.in > $@ 2>&1
# ============================================================================== # ==============================================================================


LIN/xlintsts: $(FRCLIN) $(FRC) LIN/xlintsts: $(FRCLIN) $(FRC)
@@ -578,15 +601,27 @@ LIN/xlintstzc: $(FRCLIN) $(FRC)
EIG/xeigtsts: $(FRCEIG) $(FRC) EIG/xeigtsts: $(FRCEIG) $(FRC)
$(MAKE) -C EIG xeigtsts $(MAKE) -C EIG xeigtsts


EIG/xdmdeigtsts: $(FRCEIG) $(FRC)
$(MAKE) -C EIG xdmdeigtsts

EIG/xeigtstc: $(FRCEIG) $(FRC) EIG/xeigtstc: $(FRCEIG) $(FRC)
$(MAKE) -C EIG xeigtstc $(MAKE) -C EIG xeigtstc


EIG/xdmdeigtstc: $(FRCEIG) $(FRC)
$(MAKE) -C EIG xdmdeigtstc

EIG/xeigtstd: $(FRCEIG) $(FRC) EIG/xeigtstd: $(FRCEIG) $(FRC)
$(MAKE) -C EIG xeigtstd $(MAKE) -C EIG xeigtstd


EIG/xdmdeigtstd: $(FRCEIG) $(FRC)
$(MAKE) -C EIG xdmdeigtstd

EIG/xeigtstz: $(FRCEIG) $(FRC) EIG/xeigtstz: $(FRCEIG) $(FRC)
$(MAKE) -C EIG xeigtstz $(MAKE) -C EIG xeigtstz


EIG/xdmdeigtstz: $(FRCEIG) $(FRC)
$(MAKE) -C EIG xdmdeigtstz

.PHONY: clean cleantest .PHONY: clean cleantest
clean: cleantest clean: cleantest
cleantest: cleantest:


+ 11
- 0
lapack-netlib/TESTING/cdmd.in View File

@@ -0,0 +1,11 @@
10
5

20
10

30
11

50
20

+ 11
- 0
lapack-netlib/TESTING/ddmd.in View File

@@ -0,0 +1,11 @@
10
5

20
10

30
11

50
20

+ 11
- 0
lapack-netlib/TESTING/sdmd.in View File

@@ -0,0 +1,11 @@
10
5

20
10

30
11

50
20

+ 11
- 0
lapack-netlib/TESTING/zdmd.in View File

@@ -0,0 +1,11 @@
10
5

20
10

30
11

50
20

+ 43
- 42
lapack-netlib/lapack_testing.py View File

@@ -1,31 +1,29 @@
#! /usr/bin/env python
# -*- coding: utf-8 -*-
#!/usr/bin/env python3




############################################################################### ###############################################################################
# lapack_testing.py # lapack_testing.py
############################################################################### ###############################################################################


from __future__ import print_function
from subprocess import Popen, STDOUT, PIPE from subprocess import Popen, STDOUT, PIPE
import os, sys, math import os, sys, math
import getopt import getopt
# Arguments # Arguments
try: try:
opts, args = getopt.getopt(sys.argv[1:], "hd:b:srep:t:n", opts, args = getopt.getopt(sys.argv[1:], "hd:b:srep:t:n",
["help", "dir", "bin", "short", "run", "error","prec=","test=","number"])
["help", "dir=", "bin=", "short", "run", "error","prec=","test=","number"])


except getopt.error as msg: except getopt.error as msg:
print(msg) print(msg)
print("for help use --help") print("for help use --help")
sys.exit(2) sys.exit(2)


short_summary=0
with_file=1
just_errors = 0
short_summary = False
with_file = True
just_errors = False
prec='x' prec='x'
test='all' test='all'
only_numbers=0
only_numbers = False
test_dir='TESTING' test_dir='TESTING'
bin_dir='bin/Release' bin_dir='bin/Release'


@@ -34,10 +32,9 @@ for o, a in opts:
print(sys.argv[0]+" [-h|--help] [-d dir |--dir dir] [-s |--short] [-r |--run] [-e |--error] [-p p |--prec p] [-t test |--test test] [-n | --number]") print(sys.argv[0]+" [-h|--help] [-d dir |--dir dir] [-s |--short] [-r |--run] [-e |--error] [-p p |--prec p] [-t test |--test test] [-n | --number]")
print(" - h is to print this message") print(" - h is to print this message")
print(" - r is to use to run the LAPACK tests then analyse the output (.out files). By default, the script will not run all the LAPACK tests") print(" - r is to use to run the LAPACK tests then analyse the output (.out files). By default, the script will not run all the LAPACK tests")
print(" - d [dir] is to indicate where is the LAPACK testing directory (.out files). By default, the script will use .")
print(" - b [bin] is to indicate where is the LAPACK binary files are located. By default, the script will use .")
print(" - d [dir] indicates the location of the LAPACK testing directory (.out files). By default, the script will use {:s}.".format(test_dir))
print(" - b [bin] indicates the location of the LAPACK binary files. By default, the script will use {:s}.".format(bin_dir))
print(" LEVEL OF OUTPUT") print(" LEVEL OF OUTPUT")
print(" - x is to print a detailed summary")
print(" - e is to print only the error summary") print(" - e is to print only the error summary")
print(" - s is to print a short summary") print(" - s is to print a short summary")
print(" - n is to print the numbers of failing tests (turn on summary mode)") print(" - n is to print the numbers of failing tests (turn on summary mode)")
@@ -63,15 +60,14 @@ for o, a in opts:
print(" Will return the numbers of failed tests in REAL precision by running the LAPACK Tests then analyzing the output") print(" Will return the numbers of failed tests in REAL precision by running the LAPACK Tests then analyzing the output")
print(" ./lapack_testing.py -n -p s -t eig ") print(" ./lapack_testing.py -n -p s -t eig ")
print(" Will return the numbers of failed tests in REAL precision by analyzing only the LAPACK output of EIGEN testings") print(" Will return the numbers of failed tests in REAL precision by analyzing only the LAPACK output of EIGEN testings")
print("Written by Julie Langou (June 2011) ")
sys.exit(0) sys.exit(0)
else: else:
if o in ("-s", "--short"): if o in ("-s", "--short"):
short_summary = 1
short_summary = True
if o in ("-r", "--run"): if o in ("-r", "--run"):
with_file = 0
with_file = False
if o in ("-e", "--error"): if o in ("-e", "--error"):
just_errors = 1
just_errors = True
if o in ( '-p', '--prec' ): if o in ( '-p', '--prec' ):
prec = a prec = a
if o in ( '-b', '--bin' ): if o in ( '-b', '--bin' ):
@@ -81,12 +77,12 @@ for o, a in opts:
if o in ( '-t', '--test' ): if o in ( '-t', '--test' ):
test = a test = a
if o in ( '-n', '--number' ): if o in ( '-n', '--number' ):
only_numbers = 1
short_summary = 1
only_numbers = True
short_summary = True


# process options # process options


abs_bin_dir=os.path.normpath(os.path.join(os.getcwd(),bin_dir))
abs_bin_dir=os.path.abspath(bin_dir)


os.chdir(test_dir) os.chdir(test_dir)


@@ -108,7 +104,7 @@ def run_summary_test( f, cmdline, short_summary):
nb_test_illegal=0 nb_test_illegal=0
nb_test_info=0 nb_test_info=0


if (with_file):
if with_file:
if not os.path.exists(cmdline): if not os.path.exists(cmdline):
error_message=cmdline+" file not found" error_message=cmdline+" file not found"
r=1 r=1
@@ -145,16 +141,16 @@ def run_summary_test( f, cmdline, short_summary):
whereisrun=words_in_line.index("run)") whereisrun=words_in_line.index("run)")
nb_test_run+=int(words_in_line[whereisrun-2]) nb_test_run+=int(words_in_line[whereisrun-2])
if (line.find("out of")!=-1): if (line.find("out of")!=-1):
if (short_summary==0): print(line, end=' ')
if not short_summary: print(line, end=' ')
whereisout= words_in_line.index("out") whereisout= words_in_line.index("out")
nb_test_fail+=int(words_in_line[whereisout-1]) nb_test_fail+=int(words_in_line[whereisout-1])
if ((line.find("illegal")!=-1) or (line.find("Illegal")!=-1)): if ((line.find("illegal")!=-1) or (line.find("Illegal")!=-1)):
if (short_summary==0):print(line, end=' ')
if not short_summary: print(line, end=' ')
nb_test_illegal+=1 nb_test_illegal+=1
if (line.find(" INFO")!=-1): if (line.find(" INFO")!=-1):
if (short_summary==0):print(line, end=' ')
if not short_summary: print(line, end=' ')
nb_test_info+=1 nb_test_info+=1
if (with_file==1):
if with_file:
pipe.close() pipe.close()


f.flush(); f.flush();
@@ -169,7 +165,7 @@ try:
except IOError: except IOError:
f = sys.stdout f = sys.stdout


if (short_summary==0):
if not short_summary:
print(" ") print(" ")
print("---------------- Testing LAPACK Routines ----------------") print("---------------- Testing LAPACK Routines ----------------")
print(" ") print(" ")
@@ -203,6 +199,8 @@ elif test=='mixed':
range_prec=[1,3] range_prec=[1,3]
elif test=='rfp': elif test=='rfp':
range_test=[18] range_test=[18]
elif test=='dmd':
range_test=[20]
elif test=='eig': elif test=='eig':
range_test=list(range(16)) range_test=list(range(16))
else: else:
@@ -219,7 +217,7 @@ for dtype in range_prec:
letter = dtypes[0][dtype] letter = dtypes[0][dtype]
name = dtypes[1][dtype] name = dtypes[1][dtype]


if (short_summary==0):
if not short_summary:
print(" ") print(" ")
print("------------------------- %s ------------------------" % name) print("------------------------- %s ------------------------" % name)
print(" ") print(" ")
@@ -231,19 +229,19 @@ for dtype in range_prec:
letter+"gd",letter+"sb",letter+"sg", letter+"gd",letter+"sb",letter+"sg",
letter+"bb","glm","gqr", letter+"bb","glm","gqr",
"gsv","csd","lse", "gsv","csd","lse",
letter+"test", letter+dtypes[0][dtype-1]+"test",letter+"test_rfp"),
letter+"test", letter+dtypes[0][dtype-1]+"test",letter+"test_rfp",letter+"dmd"),
("Nonsymmetric-Eigenvalue-Problem", "Symmetric-Eigenvalue-Problem", "Symmetric-Eigenvalue-Problem-2-stage", "Singular-Value-Decomposition", ("Nonsymmetric-Eigenvalue-Problem", "Symmetric-Eigenvalue-Problem", "Symmetric-Eigenvalue-Problem-2-stage", "Singular-Value-Decomposition",
"Eigen-Condition","Nonsymmetric-Eigenvalue","Nonsymmetric-Generalized-Eigenvalue-Problem", "Eigen-Condition","Nonsymmetric-Eigenvalue","Nonsymmetric-Generalized-Eigenvalue-Problem",
"Nonsymmetric-Generalized-Eigenvalue-Problem-driver", "Symmetric-Eigenvalue-Problem", "Symmetric-Eigenvalue-Generalized-Problem", "Nonsymmetric-Generalized-Eigenvalue-Problem-driver", "Symmetric-Eigenvalue-Problem", "Symmetric-Eigenvalue-Generalized-Problem",
"Banded-Singular-Value-Decomposition-routines", "Generalized-Linear-Regression-Model-routines", "Generalized-QR-and-RQ-factorization-routines", "Banded-Singular-Value-Decomposition-routines", "Generalized-Linear-Regression-Model-routines", "Generalized-QR-and-RQ-factorization-routines",
"Generalized-Singular-Value-Decomposition-routines", "CS-Decomposition-routines", "Constrained-Linear-Least-Squares-routines", "Generalized-Singular-Value-Decomposition-routines", "CS-Decomposition-routines", "Constrained-Linear-Least-Squares-routines",
"Linear-Equation-routines", "Mixed-Precision-linear-equation-routines","RFP-linear-equation-routines"),
"Linear-Equation-routines", "Mixed-Precision-linear-equation-routines","RFP-linear-equation-routines","Dynamic-Mode-Decomposition"),
(letter+"nep", letter+"sep", letter+"se2", letter+"svd", (letter+"nep", letter+"sep", letter+"se2", letter+"svd",
letter+"ec",letter+"ed",letter+"gg", letter+"ec",letter+"ed",letter+"gg",
letter+"gd",letter+"sb",letter+"sg", letter+"gd",letter+"sb",letter+"sg",
letter+"bb",letter+"glm",letter+"gqr", letter+"bb",letter+"glm",letter+"gqr",
letter+"gsv",letter+"csd",letter+"lse", letter+"gsv",letter+"csd",letter+"lse",
letter+"test", letter+dtypes[0][dtype-1]+"test",letter+"test_rfp"),
letter+"test", letter+dtypes[0][dtype-1]+"test",letter+"test_rfp",letter+"dmd"),
) )




@@ -252,22 +250,25 @@ for dtype in range_prec:
# NEED TO SKIP SOME PRECISION (namely s and c) FOR PROTO MIXED PRECISION TESTING # NEED TO SKIP SOME PRECISION (namely s and c) FOR PROTO MIXED PRECISION TESTING
if dtest==17 and (letter=="s" or letter=="c"): if dtest==17 and (letter=="s" or letter=="c"):
continue continue
if (with_file==1):
if with_file:
cmdbase=dtests[2][dtest]+".out" cmdbase=dtests[2][dtest]+".out"
else: else:
if dtest==16: if dtest==16:
# LIN TESTS # LIN TESTS
cmdbase="LIN/xlintst"+letter+" < "+dtests[0][dtest]+".in > "+dtests[2][dtest]+".out"
cmdbase="xlintst"+letter+" < "+dtests[0][dtest]+".in > "+dtests[2][dtest]+".out"
elif dtest==17: elif dtest==17:
# PROTO LIN TESTS # PROTO LIN TESTS
cmdbase="LIN/xlintst"+letter+dtypes[0][dtype-1]+" < "+dtests[0][dtest]+".in > "+dtests[2][dtest]+".out"
cmdbase="xlintst"+letter+dtypes[0][dtype-1]+" < "+dtests[0][dtest]+".in > "+dtests[2][dtest]+".out"
elif dtest==18: elif dtest==18:
# PROTO LIN TESTS # PROTO LIN TESTS
cmdbase="LIN/xlintstrf"+letter+" < "+dtests[0][dtest]+".in > "+dtests[2][dtest]+".out"
cmdbase="xlintstrf"+letter+" < "+dtests[0][dtest]+".in > "+dtests[2][dtest]+".out"
elif dtest==20:
# DMD EIG TESTS
cmdbase="xdmdeigtst"+letter+" < "+dtests[0][dtest]+".in > "+dtests[2][dtest]+".out"
else: else:
# EIG TESTS # EIG TESTS
cmdbase="EIG/xeigtst"+letter+" < "+dtests[0][dtest]+".in > "+dtests[2][dtest]+".out"
if (not just_errors and not short_summary):
cmdbase="xeigtst"+letter+" < "+dtests[0][dtest]+".in > "+dtests[2][dtest]+".out"
if not just_errors and not short_summary:
print("Testing "+name+" "+dtests[1][dtest]+"-"+cmdbase, end=' ') print("Testing "+name+" "+dtests[1][dtest]+"-"+cmdbase, end=' ')
# Run the process: either to read the file or run the LAPACK testing # Run the process: either to read the file or run the LAPACK testing
nb_test = run_summary_test(f, cmdbase, short_summary) nb_test = run_summary_test(f, cmdbase, short_summary)
@@ -277,19 +278,19 @@ for dtype in range_prec:
list_results[3][dtype]+=nb_test[3] list_results[3][dtype]+=nb_test[3]
got_error=nb_test[1]+nb_test[2]+nb_test[3] got_error=nb_test[1]+nb_test[2]+nb_test[3]


if (not short_summary):
if (nb_test[0]>0 and just_errors==0):
if not short_summary:
if nb_test[0] > 0 and not just_errors:
print("passed: "+str(nb_test[0])) print("passed: "+str(nb_test[0]))
if (nb_test[1]>0):
if nb_test[1] > 0:
print("failing to pass the threshold: "+str(nb_test[1])) print("failing to pass the threshold: "+str(nb_test[1]))
if (nb_test[2]>0):
if nb_test[2] > 0:
print("Illegal Error: "+str(nb_test[2])) print("Illegal Error: "+str(nb_test[2]))
if (nb_test[3]>0):
if nb_test[3] > 0:
print("Info Error: "+str(nb_test[3])) print("Info Error: "+str(nb_test[3]))
if (got_error>0 and just_errors==1):
if got_error > 0 and just_errors:
print("ERROR IS LOCATED IN "+name+" "+dtests[1][dtest]+" [ "+cmdbase+" ]") print("ERROR IS LOCATED IN "+name+" "+dtests[1][dtest]+" [ "+cmdbase+" ]")
print("") print("")
if (just_errors==0):
if not just_errors:
print("") print("")
# elif (got_error>0): # elif (got_error>0):
# print dtests[2][dtest]+".out \t"+str(nb_test[1])+"\t"+str(nb_test[2])+"\t"+str(nb_test[3]) # print dtests[2][dtest]+".out \t"+str(nb_test[1])+"\t"+str(nb_test[2])+"\t"+str(nb_test[3])
@@ -307,7 +308,7 @@ for dtype in range_prec:
list_results[2][4]+=list_results[2][dtype] list_results[2][4]+=list_results[2][dtype]
list_results[3][4]+=list_results[3][dtype] list_results[3][4]+=list_results[3][dtype]


if only_numbers==1:
if only_numbers:
print(str(list_results[1][4])+"\n"+str(list_results[2][4]+list_results[3][4])) print(str(list_results[1][4])+"\n"+str(list_results[2][4]+list_results[3][4]))
else: else:
print(summary) print(summary)


Write Preview
Loading…
Cancel
Save