62 use iso_fortran_env,
only: real32
64 integer,
parameter :: WP = real32
67 REAL(KIND=wp),
PARAMETER :: one = 1.0_wp
68 REAL(KIND=wp),
PARAMETER :: zero = 0.0_wp
70 REAL(KIND=wp),
ALLOCATABLE,
DIMENSION(:,:) :: &
71 a, ac, eiga, lambda, lambdaq, f, f1, f2,&
72 z, z1, s, au, w, va, x, x0, y, y0, y1
73 REAL(KIND=wp),
ALLOCATABLE,
DIMENSION(:) :: &
74 da, dl, dr, reig, reiga, reigq, ieig, &
75 ieiga, ieigq, res, res1, resex, singvx,&
77 INTEGER ,
ALLOCATABLE,
DIMENSION(:) :: IWORK
78 REAL(KIND=wp) :: ab(2,2), wdummy(2)
79 INTEGER :: IDUMMY(2), ISEED(4), RJOBDATA(8)
80 REAL(KIND=wp) :: anorm, cond, condl, condr, dmax, eps, &
81 tol, tol2, svdiff, tmp, tmp_au, &
82 tmp_fqr, tmp_rez, tmp_rezq, tmp_zxw, &
85 INTEGER :: K, KQ, LDF, LDS, LDA, LDAU, LDW, LDX, LDY, &
86 ldz, liwork, lwork, m, n, l, lloop, nrnk
87 INTEGER :: i, iJOBREF, iJOBZ, iSCALE, INFO, KDIFF, &
88 nfail, nfail_au, nfail_f_qr, nfail_rez, &
89 nfail_rezq, nfail_svdiff, nfail_total, nfailq_total, &
90 nfail_z_xv, mode, model, moder, whtsvd
91 INTEGER iNRNK, iWHTSVD, K_TRAJ, LWMINOPT
92 CHARACTER(LEN=1) GRADE, JOBREF, JOBZ, PIVTNG, RSIGN, &
93 scale, resids, wantq, wantr
109 INTRINSIC abs, int, min, max
140 WRITE(*,*)
'L Loop Index = ', lloop
152 IF ( ( min(m,n) == 0 ) .OR. ( m < n ) )
THEN
153 WRITE(*,*)
'Bad dimensions. Required: M >= N > 0.'
187 ALLOCATE( ac(lda,m) )
190 ALLOCATE( f(ldf,n+1) )
191 ALLOCATE( f1(ldf,n+1) )
192 ALLOCATE( f2(ldf,n+1) )
194 ALLOCATE( x0(ldx,n) )
195 ALLOCATE( singvx(n) )
196 ALLOCATE( singvqx(n) )
197 ALLOCATE( y(ldy,n+1) )
198 ALLOCATE( y0(ldy,n+1) )
199 ALLOCATE( y1(m,n+1) )
201 ALLOCATE( z1(ldz,n) )
211 ALLOCATE( va(lda,m) )
212 ALLOCATE( lambda(n,2) )
213 ALLOCATE( lambdaq(n,2) )
214 ALLOCATE( eiga(m,2) )
216 ALLOCATE( au(ldau,n) )
242 ALLOCATE( iwork(2*m) )
244 CALL slatmr( m, m,
'S', iseed,
'N', da, mode, cond, &
245 dmax, rsign, grade, dl, model, condl, &
246 dr, moder, condr, pivtng, iwork, m, m, &
247 zero, -one,
'N', a, lda, iwork(m+1), info )
252 ALLOCATE(work(lwork))
254 CALL sgeev(
'N',
'V', m, ac, m, reiga, ieiga, va, m, &
255 va, m, work, lwork, info )
261 tmp = max( tmp, sqrt(reiga(i)**2+ieiga(i)**2))
265 CALL slascl(
'G', 0, 0, tmp, one, m, m, a, m, info )
266 CALL slascl(
'G', 0, 0, tmp, one, m, 2, eiga, m, info )
269 anorm =
slange(
'F', n, n, a, m, wdummy )
271 IF ( k_traj == 2 )
THEN
273 CALL slarnv(2, iseed, m, f1(1,1) )
274 f1(1:m,1) = 1.0e-10*f1(1:m,1)
276 CALL sgemv(
'N', m, m, one, a, m, f1(1,i), 1, zero, &
279 x0(1:m,1:n/2) = f1(1:m,1:n/2)
280 y0(1:m,1:n/2) = f1(1:m,2:n/2+1)
282 CALL slarnv(2, iseed, m, f1(1,1) )
284 CALL sgemv(
'N', m, m, one, a, m, f1(1,i), 1, zero, &
287 x0(1:m,n/2+1:n) = f1(1:m,1:n-n/2)
288 y0(1:m,n/2+1:n) = f1(1:m,2:n-n/2+1)
291 CALL slarnv(2, iseed, m, f(1,1) )
293 CALL sgemv(
'N', m, m, one, a, m, f(1,i), 1, zero, &
296 x0(1:m,1:n) = f(1:m,1:n)
297 y0(1:m,1:n) = f(1:m,2:n+1)
300 xnorm =
slange(
'F', m, n, x0, ldx, wdummy )
301 ynorm =
slange(
'F', m, n, y0, ldx, wdummy )
306 SELECT CASE ( ijobz )
323 SELECT CASE ( ijobref )
334 SELECT CASE ( iscale )
345 DO inrnk = -1, -2, -1
361 x(1:m,1:n) = x0(1:m,1:n)
362 y(1:m,1:n) = y0(1:m,1:n)
365 CALL sgedmd( scale, jobz, resids, jobref, whtsvd, m, &
366 n, x, ldx, y, ldy, nrnk, tol, k, reig, ieig, z, &
367 ldz, res, au, ldau, w, ldw, s, lds, wdummy, -1, &
371 ALLOCATE( iwork(liwork) )
372 lwork = int(wdummy(lwminopt))
373 ALLOCATE( work(lwork) )
376 CALL sgedmd( scale, jobz, resids, jobref, whtsvd, m, &
377 n, x, ldx, y, ldy, nrnk, tol, k, reig, ieig, z, &
378 ldz, res, au, ldau, w, ldw, s, lds, work, lwork,&
379 iwork, liwork, info )
381 singvx(1:n) = work(1:n)
384 IF ( lsame(jobz,
'V') )
THEN
390 CALL sgemm(
'N',
'N', m, k, k, one, x, ldx, w, ldw, &
394 CALL saxpy( m, -one, z(1,i), 1, z1(1,i), 1)
395 tmp = max(tmp,
snrm2( m, z1(1,i), 1 ) )
397 tmp_zxw = max(tmp_zxw, tmp )
399 IF ( tmp_zxw > 10*m*eps )
THEN
400 nfail_z_xv = nfail_z_xv + 1
406 IF ( lsame(jobref,
'R') )
THEN
414 CALL sgemm(
'N',
'N', m, k, m, one, a, lda, x, ldx, &
418 CALL saxpy( m, -one, au(1,i), 1, z1(1,i), 1)
419 tmp = max( tmp,
snrm2( m, z1(1,i),1 ) * &
420 singvx(k)/(anorm*singvx(1)) )
422 tmp_au = max( tmp_au, tmp )
424 IF ( tmp > tol2 )
THEN
425 nfail_au = nfail_au + 1
428 ELSEIF ( lsame(jobref,
'E') )
THEN
436 CALL sgemm(
'N',
'N', m, k, m, one, a, lda, au, ldau, zero, y1, m
439 IF ( ieig(i) == zero )
THEN
441 CALL saxpy( m, -reig(i), au(1,i), 1, y1(1,i), 1 )
442 resex(i) =
snrm2( m, y1(1,i), 1) /
snrm2(m,au(1,i),1)
456 CALL sgemm(
'N',
'N', m, 2, 2, -one, au(1,i), &
457 m, ab, 2, one, y1(1,i), m )
458 resex(i) =
slange(
'F', m, 2, y1(1,i), m, &
459 work )/
slange(
'F', m, 2, au(1,i), m, &
461 resex(i+1) = resex(i)
469 IF ( lsame(resids,
'R') )
THEN
473 CALL sgemm(
'N',
'N', m, k, m, one, a, lda, z, ldz, zero, y1, m
479 IF ( ieig(i) == zero )
THEN
481 CALL saxpy( m, -reig(i), z(1,i), 1, y1(1,i), 1 )
482 res1(i) =
snrm2( m, y1(1,i), 1)
496 CALL sgemm(
'N',
'N', m, 2, 2, -one, z(1,i), &
497 m, ab, 2, one, y1(1,i), m )
498 res1(i) =
slange(
'F', m, 2, y1(1,i), m, &
506 tmp = max( tmp, abs(res(i) - res1(i)) * &
507 singvx(k)/(anorm*singvx(1)) )
509 tmp_rez = max( tmp_rez, tmp )
511 IF ( tmp > tol2 )
THEN
512 nfail_rez = nfail_rez + 1
515 IF ( lsame(jobref,
'E') )
THEN
518 tmp = max( tmp, abs(res1(i) - resex(i))/(res1(i)+resex(i)) )
520 tmp_ex = max(tmp_ex,tmp)
527 lambda(i,1) = reig(i)
528 lambda(i,2) = ieig(i)
537 IF ( test_qrdmd .AND. (k_traj == 1) )
THEN
542 CALL sgedmdq( scale, jobz, resids, wantq, wantr, &
543 jobref, whtsvd, m, n+1, f1, ldf, x, ldx, y, &
544 ldy, nrnk, tol, kq, reigq, ieigq, z, ldz, &
545 res, au, ldau, w, ldw, s, lds, wdummy, &
546 -1, idummy, -1, info )
548 ALLOCATE( iwork(liwork) )
549 lwork = int(wdummy(lwminopt))
550 ALLOCATE(work(lwork))
553 CALL sgedmdq( scale, jobz, resids, wantq, wantr, &
554 jobref, whtsvd, m, n+1, f1, ldf, x, ldx, y, &
555 ldy, nrnk, tol, kq, reigq, ieigq, z, ldz, &
556 res, au, ldau, w, ldw, s, lds, &
557 work, lwork, iwork, liwork, info )
559 singvqx(1:kq) = work(min(m,n+1)+1: min(m,n+1)+kq)
568 tmp = max(tmp, abs(singvx(i)-singvqx(i)) / &
571 svdiff = max( svdiff, tmp )
572 IF ( tmp > m*n*eps )
THEN
573 nfail_svdiff = nfail_svdiff + 1
577 IF ( lsame(wantq,
'Q') .AND. lsame(wantr,
'R') )
THEN
582 CALL sgemm(
'N',
'N', m, n+1, min(m,n+1), -one, f1, &
583 ldf, y, ldy, one, f2, ldf )
584 tmp_fqr =
slange(
'F', m, n+1, f2, ldf, work ) / &
585 slange(
'F', m, n+1, f, ldf, work )
586 IF ( tmp_fqr > tol2 )
THEN
587 nfail_f_qr = nfail_f_qr + 1
592 IF ( lsame(resids,
'R') )
THEN
596 CALL sgemm(
'N',
'N', m, kq, m, one, a, m, z, m, zero, y1,
602 IF ( ieigq(i) == zero )
THEN
604 CALL saxpy( m, -reigq(i), z(1,i), 1, y1(1,i), 1 )
606 res1(i) =
snrm2( m, y1(1,i), 1)
620 CALL sgemm(
'N',
'N', m, 2, 2, -one, z(1,i), &
621 m, ab, 2, one, y1(1,i), m )
623 res1(i) =
slange(
'F', m, 2, y1(1,i), m, &
631 tmp = max( tmp, abs(res(i) - res1(i)) * &
632 singvqx(k)/(anorm*singvqx(1)) )
634 tmp_rezq = max( tmp_rezq, tmp )
635 IF ( tmp > tol2 )
THEN
636 nfail_rezq = nfail_rezq + 1
642 lambdaq(i,1) = reigq(i)
643 lambdaq(i,2) = ieigq(i)
707 WRITE(*,*)
'>>>>>>>>>>>>>>>>>>>>>>>>>>'
708 WRITE(*,*)
' Test summary for SGEDMD :'
709 WRITE(*,*)
'>>>>>>>>>>>>>>>>>>>>>>>>>>'
711 IF ( nfail_z_xv == 0 )
THEN
712 WRITE(*,*)
'>>>> Z - U*V test PASSED.'
714 WRITE(*,*)
'Z - U*V test FAILED ', nfail_z_xv,
' time(s)'
715 WRITE(*,*)
'Max error ||Z-U*V||_F was ', tmp_zxw
716 nfail_total = nfail_total + nfail_z_xv
718 IF ( nfail_au == 0 )
THEN
719 WRITE(*,*)
'>>>> A*U test PASSED. '
721 WRITE(*,*)
'A*U test FAILED ', nfail_au,
' time(s)'
722 WRITE(*,*)
'Max A*U test adjusted error measure was ', tmp_au
723 WRITE(*,*)
'It should be up to O(M*N) times EPS, EPS = ', eps
724 nfail_total = nfail_total + nfail_au
727 IF ( nfail_rez == 0 )
THEN
728 WRITE(*,*)
'>>>> Rezidual computation test PASSED.'
730 WRITE(*,*)
'Rezidual computation test FAILED ', nfail_rez,
'time(s)'
731 WRITE(*,*)
'Max residual computing test adjusted error measure was '
732 WRITE(*,*)
'It should be up to O(M*N) times EPS, EPS = ', eps
733 nfail_total = nfail_total + nfail_rez
736 IF ( nfail_total == 0 )
THEN
737 WRITE(*,*)
'>>>> SGEDMD :: ALL TESTS PASSED.'
739 WRITE(*,*) nfail_total,
'FAILURES!'
740 WRITE(*,*)
'>>>>>>>>>>>>>> SGEDMD :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
743 IF ( test_qrdmd )
THEN
745 WRITE(*,*)
'>>>>>>>>>>>>>>>>>>>>>>>>>>'
746 WRITE(*,*)
' Test summary for SGEDMDQ :'
747 WRITE(*,*)
'>>>>>>>>>>>>>>>>>>>>>>>>>>'
750 IF ( nfail_svdiff == 0 )
THEN
751 WRITE(*,*)
'>>>> SGEDMD and SGEDMDQ computed singular &
752 &values test PASSED.'
754 WRITE(*,*)
'SGEDMD and SGEDMDQ discrepancies in &
755 &the singular values unacceptable ', &
756 nfail_svdiff,
' times. Test FAILED.'
757 WRITE(*,*)
'The maximal discrepancy in the singular values (relative to the norm) was '
758 WRITE(*,*)
'It should be up to O(M*N) times EPS, EPS = ', eps
759 nfailq_total = nfailq_total + nfail_svdiff
762 IF ( nfail_f_qr == 0 )
THEN
763 WRITE(*,*)
'>>>> F - Q*R test PASSED.'
765 WRITE(*,*)
'F - Q*R test FAILED ', nfail_f_qr,
' time(s)'
766 WRITE(*,*)
'The largest relative residual was ', tmp_fqr
767 WRITE(*,*)
'It should be up to O(M*N) times EPS, EPS = ', eps
768 nfailq_total = nfailq_total + nfail_f_qr
771 IF ( nfail_rezq == 0 )
THEN
772 WRITE(*,*)
'>>>> Rezidual computation test PASSED.'
774 WRITE(*,*)
'Rezidual computation test FAILED ', nfail_rezq,
'time(s)'
775 WRITE(*,*)
'Max residual computing test adjusted error measure was '
776 WRITE(*,*)
'It should be up to O(M*N) times EPS, EPS = ', eps
777 nfailq_total = nfailq_total + nfail_rezq
780 IF ( nfailq_total == 0 )
THEN
781 WRITE(*,*)
'>>>>>>> SGEDMDQ :: ALL TESTS PASSED.'
783 WRITE(*,*) nfailq_total,
'FAILURES!'
784 WRITE(*,*)
'>>>>>>> SGEDMDQ :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
790 WRITE(*,*)
'Test completed.'
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
subroutine sgedmd(jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, w, ldw, s, lds, work, lwork, iwork, liwork, info)
SGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
subroutine sgedmdq(jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, reig, imeig, z, ldz, res, b, ldb, v, ldv, s, lds, work, lwork, iwork, liwork, info)
SGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
subroutine sgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
SGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
logical function lsame(ca, cb)
LSAME
real(wp) function snrm2(n, x, incx)
SNRM2
subroutine slatmr(m, n, dist, iseed, sym, d, mode, cond, dmax, rsign, grade, dl, model, condl, dr, moder, condr, pivtng, ipivot, kl, ku, sparse, anorm, pack, a, lda, iwork, info)
SLATMR