LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zchkdmd.f90
Go to the documentation of this file.
1! This is a test program for checking the implementations of
2! the implementations of the following subroutines
3!
4! ZGEDMD, for computation of the
5! Dynamic Mode Decomposition (DMD)
6! ZGEDMDQ, for computation of a
7! QR factorization based compressed DMD
8!
9! Developed and supported by:
10! ===========================
11! Developed and coded by Zlatko Drmac, Faculty of Science,
12! University of Zagreb; drmac@math.hr
13! In cooperation with
14! AIMdyn Inc., Santa Barbara, CA.
15! ========================================================
16! How to run the code (compiler, link info)
17! ========================================================
18! Compile as FORTRAN 90 (or later) and link with BLAS and
19! LAPACK libraries.
20! NOTE: The code is developed and tested on top of the
21! Intel MKL library (versions 2022.0.3 and 2022.2.0),
22! using the Intel Fortran compiler.
23!
24! For developers of the C++ implementation
25! ========================================================
26! See the LAPACK++ and Template Numerical Toolkit (TNT)
27!
28! Note on a development of the GPU HP implementation
29! ========================================================
30! Work in progress. See CUDA, MAGMA, SLATE.
31! NOTE: The four SVD subroutines used in this code are
32! included as a part of R&D and for the completeness.
33! This was also an opportunity to test those SVD codes.
34! If the scaling option is used all four are essentially
35! equally good. For implementations on HP platforms,
36! one can use whichever SVD is available.
37!............................................................
38
39!............................................................
40!............................................................
41!
42 PROGRAM dmd_test
43 use iso_fortran_env, only: real64
44 IMPLICIT NONE
45 integer, parameter :: WP = real64
46
47!............................................................
48 REAL(KIND=wp), PARAMETER :: one = 1.0_wp
49 REAL(KIND=wp), PARAMETER :: zero = 0.0_wp
50
51 COMPLEX(KIND=WP), PARAMETER :: ZONE = ( 1.0_wp, 0.0_wp )
52 COMPLEX(KIND=WP), PARAMETER :: ZZERO = ( 0.0_wp, 0.0_wp )
53!............................................................
54 REAL(KIND=wp), ALLOCATABLE, DIMENSION(:) :: res, &
55 res1, resex, singvx, singvqx, work
56 INTEGER , ALLOCATABLE, DIMENSION(:) :: IWORK
57 REAL(KIND=wp) :: wdummy(2)
58 INTEGER :: IDUMMY(4), ISEED(4)
59 REAL(KIND=wp) :: anorm, cond, condl, condr, eps, &
60 tol, tol2, svdiff, tmp, tmp_au, &
61 tmp_fqr, tmp_rez, tmp_rezq, tmp_zxw, &
62 tmp_ex
63
64!............................................................
65 COMPLEX(KIND=WP) :: ZMAX
66 INTEGER :: LZWORK
67 COMPLEX(KIND=WP), ALLOCATABLE, DIMENSION(:,:) :: ZA, ZAC, &
68 zau, zf, zf0, zf1, zs, zw, &
69 zx, zx0, zy, zy0, zy1, zz, zz1
70 COMPLEX(KIND=WP), ALLOCATABLE, DIMENSION(:) :: ZDA, ZDR, &
71 zdl, zeigs, zeigsa, zwork
72 COMPLEX(KIND=WP) :: ZDUMMY(22), ZDUM2X2(2,2)
73!............................................................
74 INTEGER :: K, KQ, LDF, LDS, LDA, LDAU, LDW, LDX, LDY, &
75 ldz, liwork, lwork, m, n, lloop, nrnk, nrnksp
76 INTEGER :: i, iJOBREF, iJOBZ, iSCALE, INFO, j, &
77 nfail, nfail_au, nfail_f_qr, nfail_rez, &
78 nfail_rezq, nfail_svdiff, nfail_total, nfailq_total, &
79 nfail_z_xv, mode, model, moder, whtsvd, &
80 whtsvdsp
81 INTEGER :: iNRNK, iWHTSVD, K_TRAJ, LWMINOPT
82 CHARACTER :: GRADE, JOBREF, JOBZ, PIVTNG, RSIGN, &
83 scale, resids, wantq, wantr
84 LOGICAL :: TEST_QRDMD
85
86!.....external subroutines (BLAS and LAPACK)
87 EXTERNAL daxpy, dgeev, dgemm, dgemv, dlacpy, dlascl
88 EXTERNAL zgeev, zgemv, zlascl
89 EXTERNAL zlarnv, zlatmr
90 EXTERNAL zaxpy, zgemm
91!.....external subroutines DMD package, part 1
92! subroutines under test
93 EXTERNAL zgedmd, zgedmdq
94!.....external functions (BLAS and LAPACK)
95 EXTERNAL dlamch, dznrm2
96 REAL(KIND=wp) :: dlamch, dznrm2
97 REAL(KIND=wp) :: zlange
98 EXTERNAL izamax
99 INTEGER IZAMAX
100 EXTERNAL lsame
101 LOGICAL LSAME
102
103 INTRINSIC abs, int, min, max, sign
104!............................................................
105
106 ! The test is always in pairs : ( ZGEDMD and ZGEDMDQ )
107 ! because the test includes comparing the results (in pairs).
108!.....................................................................................
109 test_qrdmd = .true. ! This code by default performs tests on ZGEDMDQ
110 ! Since the QR factorizations based algorithm is designed for
111 ! single trajectory data, only single trajectory tests will
112 ! be performed with xGEDMDQ.
113 wantq = 'Q'
114 wantr = 'R'
115!.................................................................................
116
117 eps = dlamch( 'P' ) ! machine precision DP
118
119 ! Global counters of failures of some particular tests
120 nfail = 0
121 nfail_rez = 0
122 nfail_rezq = 0
123 nfail_z_xv = 0
124 nfail_f_qr = 0
125 nfail_au = 0
126 nfail_svdiff = 0
127 nfail_total = 0
128 nfailq_total = 0
129
130 DO lloop = 1, 4
131
132 WRITE(*,*) 'L Loop Index = ', lloop
133
134 ! Set the dimensions of the problem ...
135 WRITE(*,*) 'M = '
136 READ(*,*) m
137 WRITE(*,*) m
138 ! ... and the number of snapshots.
139 WRITE(*,*) 'N = '
140 READ(*,*) n
141 WRITE(*,*) n
142
143 ! ... Test the dimensions
144 IF ( ( min(m,n) == 0 ) .OR. ( m < n ) ) THEN
145 WRITE(*,*) 'Bad dimensions. Required: M >= N > 0.'
146 stop
147 END IF
148!.............
149 ! The seed inside the LLOOP so that each pass can be reproduced easily.
150 iseed(1) = 4
151 iseed(2) = 3
152 iseed(3) = 2
153 iseed(4) = 1
154
155 lda = m
156 ldf = m
157 ldx = m
158 ldy = m
159 ldw = n
160 ldz = m
161 ldau = m
162 lds = n
163
164 tmp_zxw = zero
165 tmp_au = zero
166 tmp_rez = zero
167 tmp_rezq = zero
168 svdiff = zero
169 tmp_ex = zero
170
171 ALLOCATE( za(lda,m) )
172 ALLOCATE( zac(lda,m) )
173 ALLOCATE( zf(ldf,n+1) )
174 ALLOCATE( zf0(ldf,n+1) )
175 ALLOCATE( zf1(ldf,n+1) )
176 ALLOCATE( zx(ldx,n) )
177 ALLOCATE( zx0(ldx,n) )
178 ALLOCATE( zy(ldy,n+1) )
179 ALLOCATE( zy0(ldy,n+1) )
180 ALLOCATE( zy1(ldy,n+1) )
181 ALLOCATE( zau(ldau,n) )
182 ALLOCATE( zw(ldw,n) )
183 ALLOCATE( zs(lds,n) )
184 ALLOCATE( zz(ldz,n) )
185 ALLOCATE( zz1(ldz,n) )
186 ALLOCATE( res(n) )
187 ALLOCATE( res1(n) )
188 ALLOCATE( resex(n) )
189 ALLOCATE( zeigs(n) )
190 ALLOCATE( singvx(n) )
191 ALLOCATE( singvqx(n) )
192
193 tol = m*eps
194 ! This mimics O(M*N)*EPS bound for accumulated roundoff error.
195 ! The factor 10 is somewhat arbitrary.
196 tol2 = 10*m*n*eps
197
198!.............
199
200 DO k_traj = 1, 2
201 ! Number of intial conditions in the simulation/trajectories (1 or 2)
202
203 cond = 1.0d4
204 zmax = (1.0d1,1.0d1)
205 rsign = 'F'
206 grade = 'N'
207 model = 6
208 condl = 1.0d1
209 moder = 6
210 condr = 1.0d1
211 pivtng = 'N'
212
213 ! Loop over all parameter MODE values for ZLATMR (+1,..,+6)
214 DO mode = 1, 6
215
216 ALLOCATE( iwork(2*m) )
217 ALLOCATE( zda(m) )
218 ALLOCATE( zdl(m) )
219 ALLOCATE( zdr(m) )
220
221 CALL zlatmr( m, m, 'N', iseed, 'N', zda, mode, cond, &
222 zmax, rsign, grade, zdl, model, condl, &
223 zdr, moder, condr, pivtng, iwork, m, m, &
224 zero, -one, 'N', za, lda, iwork(m+1), info )
225 DEALLOCATE( zdr )
226 DEALLOCATE( zdl )
227 DEALLOCATE( zda )
228 DEALLOCATE( iwork )
229
230 lzwork = max(1,2*m)
231 ALLOCATE( zeigsa(m) )
232 ALLOCATE( zwork(lzwork) )
233 ALLOCATE( work(2*m) )
234 zac(1:m,1:m) = za(1:m,1:m)
235 CALL zgeev( 'N','N', m, zac, lda, zeigsa, zdum2x2, 2, &
236 zdum2x2, 2, zwork, lzwork, work, info ) ! LAPACK CALL
237 DEALLOCATE(work)
238 DEALLOCATE(zwork)
239
240 tmp = abs(zeigsa(izamax(m, zeigsa, 1))) ! The spectral radius of ZA
241 ! Scale the matrix ZA to have unit spectral radius.
242 CALL zlascl( 'G',0, 0, tmp, one, m, m, &
243 za, lda, info )
244 CALL zlascl( 'G',0, 0, tmp, one, m, 1, &
245 zeigsa, m, info )
246 anorm = zlange( 'F', m, m, za, lda, wdummy )
247
248 IF ( k_traj == 2 ) THEN
249 ! generate data as two trajectories
250 ! with two inital conditions
251 CALL zlarnv(2, iseed, m, zf(1,1) )
252 DO i = 1, n/2
253 CALL zgemv( 'N', m, m, zone, za, lda, zf(1,i), 1, &
254 zzero, zf(1,i+1), 1 )
255 END DO
256 zx0(1:m,1:n/2) = zf(1:m,1:n/2)
257 zy0(1:m,1:n/2) = zf(1:m,2:n/2+1)
258
259 CALL zlarnv(2, iseed, m, zf(1,1) )
260 DO i = 1, n-n/2
261 CALL zgemv( 'N', m, m, zone, za, lda, zf(1,i), 1, &
262 zzero, zf(1,i+1), 1 )
263 END DO
264 zx0(1:m,n/2+1:n) = zf(1:m,1:n-n/2)
265 zy0(1:m,n/2+1:n) = zf(1:m,2:n-n/2+1)
266 ELSE
267 CALL zlarnv(2, iseed, m, zf(1,1) )
268 DO i = 1, n
269 CALL zgemv( 'N', m, m, zone, za, m, zf(1,i), 1, &
270 zzero, zf(1,i+1), 1 )
271 END DO
272 zf0(1:m,1:n+1) = zf(1:m,1:n+1)
273 zx0(1:m,1:n) = zf0(1:m,1:n)
274 zy0(1:m,1:n) = zf0(1:m,2:n+1)
275 END IF
276
277 DEALLOCATE( zeigsa )
278!........................................................................
279
280 DO ijobz = 1, 4
281
282 SELECT CASE ( ijobz )
283 CASE(1)
284 jobz = 'V'
285 resids = 'R'
286 CASE(2)
287 jobz = 'V'
288 resids = 'N'
289 CASE(3)
290 jobz = 'F'
291 resids = 'N'
292 CASE(4)
293 jobz = 'N'
294 resids = 'N'
295 END SELECT
296
297 DO ijobref = 1, 3
298
299 SELECT CASE ( ijobref )
300 CASE(1)
301 jobref = 'R'
302 CASE(2)
303 jobref = 'E'
304 CASE(3)
305 jobref = 'N'
306 END SELECT
307
308 DO iscale = 1, 4
309
310 SELECT CASE ( iscale )
311 CASE(1)
312 scale = 'S'
313 CASE(2)
314 scale = 'C'
315 CASE(3)
316 scale = 'Y'
317 CASE(4)
318 scale = 'N'
319 END SELECT
320
321 DO inrnk = -1, -2, -1
322 nrnk = inrnk
323 nrnksp = inrnk
324
325 DO iwhtsvd = 1, 3
326 ! Check all four options to compute the POD basis
327 ! via the SVD.
328 whtsvd = iwhtsvd
329 whtsvdsp = iwhtsvd
330
331 DO lwminopt = 1, 2
332 ! Workspace query for the minimal (1) and for the optimal
333 ! (2) workspace lengths determined by workspace query.
334
335 ! ZGEDMD is always tested and its results are also used for
336 ! comparisons with ZGEDMDQ.
337
338 zx(1:m,1:n) = zx0(1:m,1:n)
339 zy(1:m,1:n) = zy0(1:m,1:n)
340
341 CALL zgedmd( scale, jobz, resids, jobref, whtsvd, &
342 m, n, zx, ldx, zy, ldy, nrnk, tol, &
343 k, zeigs, zz, ldz, res, zau, ldau, &
344 zw, ldw, zs, lds, zdummy, -1, &
345 wdummy, -1, idummy, -1, info )
346 IF ( (info .EQ. 2) .OR. ( info .EQ. 3 ) &
347 .OR. ( info < 0 ) ) THEN
348 WRITE(*,*) 'Call to ZGEDMD workspace query failed. &
349 &Check the calling sequence and the code.'
350 WRITE(*,*) 'The error code is ', info
351 WRITE(*,*) 'The input parameters were ', &
352 scale, jobz, resids, jobref, whtsvd, &
353 m, n, ldx, ldy, nrnk, tol, ldz, ldau, ldw, lds
354 stop
355 END IF
356
357 lzwork = int(zdummy(lwminopt))
358 lwork = int(wdummy(1))
359 liwork = idummy(1)
360
361 ALLOCATE(zwork(lzwork))
362 ALLOCATE( work(lwork))
363 ALLOCATE(iwork(liwork))
364
365 CALL zgedmd( scale, jobz, resids, jobref, whtsvd, &
366 m, n, zx, ldx, zy, ldy, nrnk, tol, &
367 k, zeigs, zz, ldz, res, zau, ldau, &
368 zw, ldw, zs, lds, zwork, lzwork, &
369 work, lwork, iwork, liwork, info )
370
371 IF ( info /= 0 ) THEN
372 WRITE(*,*) 'Call to ZGEDMD failed. &
373 &Check the calling sequence and the code.'
374 WRITE(*,*) 'The error code is ', info
375 WRITE(*,*) 'The input parameters were ',&
376 scale, jobz, resids, jobref, whtsvd, &
377 m, n, ldx, ldy, nrnk, tol
378 stop
379 END IF
380
381 singvx(1:n) = work(1:n)
382
383 !...... ZGEDMD check point
384 IF ( lsame(jobz,'V') ) THEN
385 ! Check that Z = X*W, on return from ZGEDMD
386 ! This checks that the returned eigenvectors in Z are
387 ! the product of the SVD'POD basis returned in X
388 ! and the eigenvectors of the rayleigh quotient
389 ! returned in W
390 CALL zgemm( 'N', 'N', m, k, k, zone, zx, ldx, zw, ldw, &
391 zzero, zz1, ldz )
392 tmp = zero
393 DO i = 1, k
394 CALL zaxpy( m, -zone, zz(1,i), 1, zz1(1,i), 1)
395 tmp = max(tmp, dznrm2( m, zz1(1,i), 1 ) )
396 END DO
397 tmp_zxw = max(tmp_zxw, tmp )
398 IF ( tmp_zxw <= 10*m*eps ) THEN
399 !WRITE(*,*) ' :) .... OK .........ZGEDMD PASSED.'
400 ELSE
401 nfail_z_xv = nfail_z_xv + 1
402 WRITE(*,*) ':( .................ZGEDMD FAILED!', &
403 'Check the code for implementation errors.'
404 WRITE(*,*) 'The input parameters were ',&
405 scale, jobz, resids, jobref, whtsvd, &
406 m, n, ldx, ldy, nrnk, tol
407 END IF
408 END IF
409
410
411 !...... ZGEDMD check point
412 IF ( lsame(jobref,'R') ) THEN
413 ! The matrix A*U is returned for computing refined Ritz vectors.
414 ! Check that A*U is computed correctly using the formula
415 ! A*U = Y * V * inv(SIGMA). This depends on the
416 ! accuracy in the computed singular values and vectors of X.
417 ! See the paper for an error analysis.
418 ! Note that the left singular vectors of the input matrix X
419 ! are returned in the array X.
420 CALL zgemm( 'N', 'N', m, k, m, zone, za, lda, zx, ldx, &
421 zzero, zz1, ldz )
422 tmp = zero
423 DO i = 1, k
424 CALL zaxpy( m, -zone, zau(1,i), 1, zz1(1,i), 1)
425 tmp = max( tmp, dznrm2( m, zz1(1,i),1 ) * &
426 singvx(k)/(anorm*singvx(1)) )
427 END DO
428 tmp_au = max( tmp_au, tmp )
429 IF ( tmp <= tol2 ) THEN
430 !WRITE(*,*) ':) .... OK .........ZGEDMD PASSED.'
431 ELSE
432 nfail_au = nfail_au + 1
433 WRITE(*,*) ':( .................ZGEDMD FAILED!', &
434 'Check the code for implementation errors.'
435 WRITE(*,*) 'The input parameters were ',&
436 scale, jobz, resids, jobref, whtsvd, &
437 m, n, ldx, ldy, nrnk, tol
438 END IF
439 ELSEIF ( lsame(jobref,'E') ) THEN
440 ! The unscaled vectors of the Exact DMD are computed.
441 ! This option is included for the sake of completeness,
442 ! for users who prefer the Exact DMD vectors. The
443 ! returned vectors are in the real form, in the same way
444 ! as the Ritz vectors. Here we just save the vectors
445 ! and test them separately using a Matlab script.
446
447
448 CALL zgemm( 'N', 'N', m, k, m, zone, za, lda, zau, ldau, zzero, zy1, ldy )
449
450 DO i=1, k
451 ! have a real eigenvalue with real eigenvector
452 CALL zaxpy( m, -zeigs(i), zau(1,i), 1, zy1(1,i), 1 )
453 resex(i) = dznrm2( m, zy1(1,i), 1) / dznrm2(m,zau(1,i),1)
454 END DO
455 END IF
456 !...... ZGEDMD check point
457
458 IF ( lsame(resids, 'R') ) THEN
459 ! Compare the residuals returned by ZGEDMD with the
460 ! explicitly computed residuals using the matrix A.
461 ! Compute explicitly Y1 = A*Z
462 CALL zgemm( 'N', 'N', m, k, m, zone, za, lda, zz, ldz, zzero, zy1, ldy )
463 ! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
464 ! of the invariant subspaces that correspond to complex conjugate
465 ! pairs of eigencalues. (See the description of Z in ZGEDMD,)
466
467 DO i=1, k
468 ! have a real eigenvalue with real eigenvector
469 CALL zaxpy( m, -zeigs(i), zz(1,i), 1, zy1(1,i), 1 )
470 res1(i) = dznrm2( m, zy1(1,i), 1)
471 END DO
472 tmp = zero
473 DO i = 1, k
474 tmp = max( tmp, abs(res(i) - res1(i)) * &
475 singvx(k)/(anorm*singvx(1)) )
476 END DO
477 tmp_rez = max( tmp_rez, tmp )
478 IF ( tmp <= tol2 ) THEN
479 !WRITE(*,*) ':) .... OK ..........ZGEDMD PASSED.'
480 ELSE
481 nfail_rez = nfail_rez + 1
482 WRITE(*,*) ':( ..................ZGEDMD FAILED!', &
483 'Check the code for implementation errors.'
484 WRITE(*,*) 'The input parameters were ',&
485 scale, jobz, resids, jobref, whtsvd, &
486 m, n, ldx, ldy, nrnk, tol
487 END IF
488
489
490 IF ( lsame(jobref,'E') ) THEN
491 tmp = zero
492 DO i = 1, k
493 tmp = max( tmp, abs(res1(i) - resex(i))/(res1(i)+resex(i)) )
494 END DO
495 tmp_ex = max(tmp_ex,tmp)
496 END IF
497
498 END IF
499
500 DEALLOCATE(zwork)
501 DEALLOCATE(work)
502 DEALLOCATE(iwork)
503
504 IF ( test_qrdmd .AND. (k_traj == 1) ) THEN
505
506 zf(1:m,1:n+1) = zf0(1:m,1:n+1)
507
508 CALL zgedmdq( scale, jobz, resids, wantq, wantr, jobref, &
509 whtsvd, m, n+1, zf, ldf, zx, ldx, zy, ldy, &
510 nrnk, tol, k, zeigs, zz, ldz, res, zau, &
511 ldau, zw, ldw, zs, lds, zdummy, -1, &
512 wdummy, -1, idummy, -1, info )
513
514 lzwork = int(zdummy(lwminopt))
515 ALLOCATE( zwork(lzwork) )
516 liwork = idummy(1)
517 ALLOCATE(iwork(liwork))
518 lwork = int(wdummy(1))
519 ALLOCATE(work(lwork))
520
521 CALL zgedmdq( scale, jobz, resids, wantq, wantr, jobref, &
522 whtsvd, m, n+1, zf, ldf, zx, ldx, zy, ldy, &
523 nrnk, tol, kq, zeigs, zz, ldz, res, zau, &
524 ldau, zw, ldw, zs, lds, zwork, lzwork, &
525 work, lwork, iwork, liwork, info )
526
527 IF ( info /= 0 ) THEN
528 WRITE(*,*) 'Call to ZGEDMDQ failed. &
529 &Check the calling sequence and the code.'
530 WRITE(*,*) 'The error code is ', info
531 WRITE(*,*) 'The input parameters were ',&
532 scale, jobz, resids, wantq, wantr, whtsvd, &
533 m, n, ldx, ldy, nrnk, tol
534 stop
535 END IF
536 singvqx(1:n) = work(1:n)
537
538 !..... ZGEDMDQ check point
539
540 IF ( 1 == 0 ) THEN
541 ! Comparison of ZGEDMD and ZGEDMDQ singular values disabled
542 tmp = zero
543 DO i = 1, min(k, kq)
544 tmp = max(tmp, abs(singvx(i)-singvqx(i)) / &
545 singvx(1) )
546 END DO
547 svdiff = max( svdiff, tmp )
548 IF ( tmp > m*n*eps ) THEN
549 WRITE(*,*) 'FAILED! Something was wrong with the run.'
550 nfail_svdiff = nfail_svdiff + 1
551 DO j =1, 3
552 write(*,*) j, singvx(j), singvqx(j)
553 read(*,*)
554 END DO
555
556 END IF
557 END IF
558
559 !..... ZGEDMDQ check point
560 IF ( lsame(wantq,'Q') .AND. lsame(wantr,'R') ) THEN
561 ! Check that the QR factors are computed and returned
562 ! as requested. The residual ||F-Q*R||_F / ||F||_F
563 ! is compared to M*N*EPS.
564 zf1(1:m,1:n+1) = zf0(1:m,1:n+1)
565 CALL zgemm( 'N', 'N', m, n+1, min(m,n+1), -zone, zf, &
566 ldf, zy, ldy, zone, zf1, ldf )
567 tmp_fqr = zlange( 'F', m, n+1, zf1, ldf, work ) / &
568 zlange( 'F', m, n+1, zf0, ldf, work )
569 IF ( tmp_fqr > tol2 ) THEN
570 WRITE(*,*) 'FAILED! Something was wrong with the run.'
571 nfail_f_qr = nfail_f_qr + 1
572 ELSE
573 !WRITE(*,*) '........ PASSED.'
574 END IF
575 END IF
576
577 !..... ZGEDMDQ check point
578 IF ( lsame(resids, 'R') ) THEN
579 ! Compare the residuals returned by ZGEDMDQ with the
580 ! explicitly computed residuals using the matrix A.
581 ! Compute explicitly Y1 = A*Z
582 CALL zgemm( 'N', 'N', m, kq, m, zone, za, lda, zz, ldz, zzero, zy1, ldy )
583 ! ... and then A*Z(:,i) - LAMBDA(i)*Z(:,i), using the real forms
584 ! of the invariant subspaces that correspond to complex conjugate
585 ! pairs of eigencalues. (See the description of Z in ZGEDMDQ)
586
587 DO i=1, kq
588 ! have a real eigenvalue with real eigenvector
589 CALL zaxpy( m, -zeigs(i), zz(1,i), 1, zy1(1,i), 1 )
590 ! Y(1:M,i) = Y(1:M,i) - REIG(i)*Z(1:M,i)
591 res1(i) = dznrm2( m, zy1(1,i), 1)
592 END DO
593 tmp = zero
594 DO i = 1, kq
595 tmp = max( tmp, abs(res(i) - res1(i)) * &
596 singvqx(kq)/(anorm*singvqx(1)) )
597 END DO
598 tmp_rezq = max( tmp_rezq, tmp )
599 IF ( tmp <= tol2 ) THEN
600 !WRITE(*,*) '.... OK ........ ZGEDMDQ PASSED.'
601 ELSE
602 nfail_rezq = nfail_rezq + 1
603 WRITE(*,*) '................ ZGEDMDQ FAILED!', &
604 'Check the code for implementation errors.'
605 stop
606 END IF
607
608 END IF
609
610 DEALLOCATE( zwork )
611 DEALLOCATE( work )
612 DEALLOCATE( iwork )
613
614 END IF ! ZGEDMDQ
615
616!.......................................................................................................
617
618 END DO ! LWMINOPT
619 !write(*,*) 'LWMINOPT loop completed'
620 END DO ! iWHTSVD
621 !write(*,*) 'WHTSVD loop completed'
622 END DO ! iNRNK -2:-1
623 !write(*,*) 'NRNK loop completed'
624 END DO ! iSCALE 1:4
625 !write(*,*) 'SCALE loop completed'
626 END DO
627 !write(*,*) 'JOBREF loop completed'
628 END DO ! iJOBZ
629 !write(*,*) 'JOBZ loop completed'
630
631 END DO ! MODE -6:6
632 !write(*,*) 'MODE loop completed'
633 END DO ! 1 or 2 trajectories
634 !write(*,*) 'trajectories loop completed'
635
636 DEALLOCATE( za )
637 DEALLOCATE( zac )
638 DEALLOCATE( zz )
639 DEALLOCATE( zf )
640 DEALLOCATE( zf0 )
641 DEALLOCATE( zf1 )
642 DEALLOCATE( zx )
643 DEALLOCATE( zx0 )
644 DEALLOCATE( zy )
645 DEALLOCATE( zy0 )
646 DEALLOCATE( zy1 )
647 DEALLOCATE( zau )
648 DEALLOCATE( zw )
649 DEALLOCATE( zs )
650 DEALLOCATE( zz1 )
651 DEALLOCATE( res )
652 DEALLOCATE( res1 )
653 DEALLOCATE( resex )
654 DEALLOCATE( zeigs )
655 DEALLOCATE( singvx )
656 DEALLOCATE( singvqx )
657
658 END DO ! LLOOP
659
660 WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
661 WRITE(*,*) ' Test summary for ZGEDMD :'
662 WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
663 WRITE(*,*)
664 IF ( nfail_z_xv == 0 ) THEN
665 WRITE(*,*) '>>>> Z - U*V test PASSED.'
666 ELSE
667 WRITE(*,*) 'Z - U*V test FAILED ', nfail_z_xv, ' time(s)'
668 WRITE(*,*) 'Max error ||Z-U*V||_F was ', tmp_zxw
669 nfail_total = nfail_total + nfail_z_xv
670 END IF
671 IF ( nfail_au == 0 ) THEN
672 WRITE(*,*) '>>>> A*U test PASSED. '
673 ELSE
674 WRITE(*,*) 'A*U test FAILED ', nfail_au, ' time(s)'
675 WRITE(*,*) 'Max A*U test adjusted error measure was ', tmp_au
676 WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
677 nfail_total = nfail_total + nfail_au
678 END IF
679
680 IF ( nfail_rez == 0 ) THEN
681 WRITE(*,*) '>>>> Rezidual computation test PASSED.'
682 ELSE
683 WRITE(*,*) 'Rezidual computation test FAILED ', nfail_rez, 'time(s)'
684 WRITE(*,*) 'Max residual computing test adjusted error measure was ', tmp_rez
685 WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
686 nfail_total = nfail_total + nfail_rez
687 END IF
688
689 IF ( nfail_total == 0 ) THEN
690 WRITE(*,*) '>>>> ZGEDMD :: ALL TESTS PASSED.'
691 ELSE
692 WRITE(*,*) nfail_total, 'FAILURES!'
693 WRITE(*,*) '>>>>>>>>>>>>>> ZGEDMD :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
694 END IF
695
696 IF ( test_qrdmd ) THEN
697 WRITE(*,*)
698 WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
699 WRITE(*,*) ' Test summary for ZGEDMDQ :'
700 WRITE(*,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>'
701 WRITE(*,*)
702
703 IF ( nfail_svdiff == 0 ) THEN
704 WRITE(*,*) '>>>> ZGEDMD and ZGEDMDQ computed singular &
705 &values test PASSED.'
706 ELSE
707 WRITE(*,*) 'ZGEDMD and ZGEDMDQ discrepancies in &
708 &the singular values unacceptable ', &
709 nfail_svdiff, ' times. Test FAILED.'
710 WRITE(*,*) 'The maximal discrepancy in the singular values (relative to the norm) was ', svdiff
711 WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
712 nfailq_total = nfailq_total + nfail_svdiff
713 END IF
714
715 IF ( nfail_f_qr == 0 ) THEN
716 WRITE(*,*) '>>>> F - Q*R test PASSED.'
717 ELSE
718 WRITE(*,*) 'F - Q*R test FAILED ', nfail_f_qr, ' time(s)'
719 WRITE(*,*) 'The largest relative residual was ', tmp_fqr
720 WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
721 nfailq_total = nfailq_total + nfail_f_qr
722 END IF
723
724 IF ( nfail_rezq == 0 ) THEN
725 WRITE(*,*) '>>>> Rezidual computation test PASSED.'
726 ELSE
727 WRITE(*,*) 'Rezidual computation test FAILED ', nfail_rezq, 'time(s)'
728 WRITE(*,*) 'Max residual computing test adjusted error measure was ', tmp_rezq
729 WRITE(*,*) 'It should be up to O(M*N) times EPS, EPS = ', eps
730 nfailq_total = nfailq_total + nfail_rezq
731 END IF
732
733 IF ( nfailq_total == 0 ) THEN
734 WRITE(*,*) '>>>>>>> ZGEDMDQ :: ALL TESTS PASSED.'
735 ELSE
736 WRITE(*,*) nfailq_total, 'FAILURES!'
737 WRITE(*,*) '>>>>>>> ZGEDMDQ :: TESTS FAILED. CHECK THE IMPLEMENTATION.'
738 END IF
739
740 END IF
741
742 WRITE(*,*)
743 WRITE(*,*) 'Test completed.'
744 stop
745 END
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zgedmdq
subroutine zgedmdq(jobs, jobz, jobr, jobq, jobt, jobf, whtsvd, m, n, f, ldf, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, v, ldv, s, lds, zwork, lzwork, work, lwork, iwork, liwork, info)
ZGEDMDQ computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
Definition zgedmdq.f90:556
zgedmd
subroutine zgedmd(jobs, jobz, jobr, jobf, whtsvd, m, n, x, ldx, y, ldy, nrnk, tol, k, eigs, z, ldz, res, b, ldb, w, ldw, s, lds, zwork, lzwork, rwork, lrwork, iwork, liwork, info)
ZGEDMD computes the Dynamic Mode Decomposition (DMD) for a pair of data snapshot matrices.
Definition zgedmd.f90:498
dgeev
subroutine dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition dgeev.f:192
zgeev
subroutine zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
ZGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition zgeev.f:180
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
izamax
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlange
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlange.f:115
zlarnv
subroutine zlarnv(idist, iseed, n, x)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition zlarnv.f:99
zlascl
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition zlascl.f:143
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dznrm2
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition dznrm2.f90:90
zlatmr
subroutine zlatmr(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)
ZLATMR
Definition zlatmr.f:490