Skip to content

Commit 14f3af6

Browse files
langoulangou
authored
Merge pull request #1234 from nakatamaho/cleanup/csdts
cleanup: separate real and complex zero/one constants in {c,z}csdts
2 parents 1c493c4 + c21b6a6 commit 14f3af6

2 files changed

Lines changed: 136 additions & 137 deletions

File tree

‎TESTING/EIG/ccsdts.f‎

Lines changed: 68 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -246,10 +246,10 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
246246
* =====================================================================
247247
*
248248
* .. Parameters ..
249-
REAL REALONE, REALZERO
250-
PARAMETER ( REALONE = 1.0E0, REALZERO = 0.0E0 )
251-
COMPLEX ZERO, ONE
252-
PARAMETER ( ZERO = (0.0E0,0.0E0), ONE = (1.0E0,0.0E0) )
249+
REAL ZERO, ONE
250+
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
251+
COMPLEX CZERO, CONE
252+
PARAMETER ( CZERO = (0.0E0,0.0E0), CONE = (1.0E0,0.0E0) )
253253
REAL PIOVER2
254254
PARAMETER ( PIOVER2 = 1.57079632679489661923132169163975144210E0 )
255255
* ..
@@ -271,13 +271,13 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
271271
* .. Executable Statements ..
272272
*
273273
ULP = SLAMCH( 'Precision' )
274-
ULPINV = REALONE / ULP
274+
ULPINV = ONE / ULP
275275
*
276276
* The first half of the routine checks the 2-by-2 CSD
277277
*
278-
CALL CLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
279-
CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -REALONE,
280-
$ X, LDX, REALONE, WORK, LDX )
278+
CALL CLASET( 'Full', M, M, CZERO, CONE, WORK, LDX )
279+
CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE,
280+
$ X, LDX, ONE, WORK, LDX )
281281
IF (M.GT.0) THEN
282282
EPS2 = MAX( ULP,
283283
$ CLANGE( '1', M, M, WORK, LDX, RWORK ) / REAL( M ) )
@@ -301,64 +301,64 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
301301
*
302302
CALL CLACPY( 'Full', M, M, X, LDX, XF, LDX )
303303
*
304-
CALL CGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
305-
$ XF, LDX, V1T, LDV1T, ZERO, WORK, LDX )
304+
CALL CGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, CONE,
305+
$ XF, LDX, V1T, LDV1T, CZERO, WORK, LDX )
306306
*
307-
CALL CGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
308-
$ U1, LDU1, WORK, LDX, ZERO, XF, LDX )
307+
CALL CGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, CONE,
308+
$ U1, LDU1, WORK, LDX, CZERO, XF, LDX )
309309
*
310310
DO I = 1, MIN(P,Q)-R
311-
XF(I,I) = XF(I,I) - ONE
311+
XF(I,I) = XF(I,I) - CONE
312312
END DO
313313
DO I = 1, R
314314
XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
315-
$ XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - CMPLX( COS(THETA(I)),
316-
$ 0.0E0 )
315+
$ XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) -
316+
$ CMPLX( COS(THETA(I)), ZERO )
317317
END DO
318318
*
319319
CALL CGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
320-
$ ONE, XF(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
320+
$ CONE, XF(1,Q+1), LDX, V2T, LDV2T, CZERO, WORK, LDX )
321321
*
322322
CALL CGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
323-
$ ONE, U1, LDU1, WORK, LDX, ZERO, XF(1,Q+1), LDX )
323+
$ CONE, U1, LDU1, WORK, LDX, CZERO, XF(1,Q+1), LDX )
324324
*
325325
DO I = 1, MIN(P,M-Q)-R
326-
XF(P-I+1,M-I+1) = XF(P-I+1,M-I+1) + ONE
326+
XF(P-I+1,M-I+1) = XF(P-I+1,M-I+1) + CONE
327327
END DO
328328
DO I = 1, R
329329
XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
330330
$ XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
331-
$ CMPLX( SIN(THETA(R-I+1)), 0.0E0 )
331+
$ CMPLX( SIN(THETA(R-I+1)), ZERO )
332332
END DO
333333
*
334-
CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
335-
$ XF(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
334+
CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q,
335+
$ CONE, XF(P+1,1), LDX, V1T, LDV1T, CZERO, WORK, LDX )
336336
*
337337
CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
338-
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,1), LDX )
338+
$ CONE, U2, LDU2, WORK, LDX, CZERO, XF(P+1,1), LDX )
339339
*
340340
DO I = 1, MIN(M-P,Q)-R
341-
XF(M-I+1,Q-I+1) = XF(M-I+1,Q-I+1) - ONE
341+
XF(M-I+1,Q-I+1) = XF(M-I+1,Q-I+1) - CONE
342342
END DO
343343
DO I = 1, R
344344
XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
345345
$ XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
346-
$ CMPLX( SIN(THETA(R-I+1)), 0.0E0 )
346+
$ CMPLX( SIN(THETA(R-I+1)), ZERO )
347347
END DO
348348
*
349349
CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
350-
$ ONE, XF(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
350+
$ CONE, XF(P+1,Q+1), LDX, V2T, LDV2T, CZERO, WORK, LDX )
351351
*
352352
CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
353-
$ ONE, U2, LDU2, WORK, LDX, ZERO, XF(P+1,Q+1), LDX )
353+
$ CONE, U2, LDU2, WORK, LDX, CZERO, XF(P+1,Q+1), LDX )
354354
*
355355
DO I = 1, MIN(M-P,M-Q)-R
356-
XF(P+I,Q+I) = XF(P+I,Q+I) - ONE
356+
XF(P+I,Q+I) = XF(P+I,Q+I) - CONE
357357
END DO
358358
DO I = 1, R
359359
XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
360360
$ XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
361-
$ CMPLX( COS(THETA(I)), 0.0E0 )
361+
$ CMPLX( COS(THETA(I)), ZERO )
362362
END DO
363363
*
364364
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
@@ -383,9 +383,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
383383
*
384384
* Compute I - U1'*U1
385385
*
386-
CALL CLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
387-
CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -REALONE,
388-
$ U1, LDU1, REALONE, WORK, LDU1 )
386+
CALL CLASET( 'Full', P, P, CZERO, CONE, WORK, LDU1 )
387+
CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE,
388+
$ U1, LDU1, ONE, WORK, LDU1 )
389389
*
390390
* Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
391391
*
@@ -394,9 +394,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
394394
*
395395
* Compute I - U2'*U2
396396
*
397-
CALL CLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
398-
CALL CHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -REALONE,
399-
$ U2, LDU2, REALONE, WORK, LDU2 )
397+
CALL CLASET( 'Full', M-P, M-P, CZERO, CONE, WORK, LDU2 )
398+
CALL CHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE,
399+
$ U2, LDU2, ONE, WORK, LDU2 )
400400
*
401401
* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
402402
*
@@ -405,9 +405,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
405405
*
406406
* Compute I - V1T*V1T'
407407
*
408-
CALL CLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
409-
CALL CHERK( 'Upper', 'No transpose', Q, Q, -REALONE,
410-
$ V1T, LDV1T, REALONE, WORK, LDV1T )
408+
CALL CLASET( 'Full', Q, Q, CZERO, CONE, WORK, LDV1T )
409+
CALL CHERK( 'Upper', 'No transpose', Q, Q, -ONE,
410+
$ V1T, LDV1T, ONE, WORK, LDV1T )
411411
*
412412
* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
413413
*
@@ -416,9 +416,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
416416
*
417417
* Compute I - V2T*V2T'
418418
*
419-
CALL CLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
420-
CALL CHERK( 'Upper', 'No transpose', M-Q, M-Q, -REALONE,
421-
$ V2T, LDV2T, REALONE, WORK, LDV2T )
419+
CALL CLASET( 'Full', M-Q, M-Q, CZERO, CONE, WORK, LDV2T )
420+
CALL CHERK( 'Upper', 'No transpose', M-Q, M-Q, -ONE,
421+
$ V2T, LDV2T, ONE, WORK, LDV2T )
422422
*
423423
* Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
424424
*
@@ -427,9 +427,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
427427
*
428428
* Check sorting
429429
*
430-
RESULT( 9 ) = REALZERO
430+
RESULT( 9 ) = ZERO
431431
DO I = 1, R
432-
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
432+
IF( THETA(I).LT.ZERO .OR. THETA(I).GT.PIOVER2 ) THEN
433433
RESULT( 9 ) = ULPINV
434434
END IF
435435
IF( I.GT.1) THEN
@@ -441,9 +441,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
441441
*
442442
* The second half of the routine checks the 2-by-1 CSD
443443
*
444-
CALL CLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDX )
445-
CALL CHERK( 'Upper', 'Conjugate transpose', Q, M, -REALONE,
446-
$ X, LDX, REALONE, WORK, LDX )
444+
CALL CLASET( 'Full', Q, Q, CZERO, CONE, WORK, LDX )
445+
CALL CHERK( 'Upper', 'Conjugate transpose', Q, M, -ONE,
446+
$ X, LDX, ONE, WORK, LDX )
447447
IF (M.GT.0) THEN
448448
EPS2 = MAX( ULP,
449449
$ CLANGE( '1', Q, Q, WORK, LDX, RWORK ) / REAL( M ) )
@@ -464,34 +464,34 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
464464
*
465465
* Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
466466
*
467-
CALL CGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
468-
$ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
467+
CALL CGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, CONE,
468+
$ X, LDX, V1T, LDV1T, CZERO, WORK, LDX )
469469
*
470-
CALL CGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
471-
$ U1, LDU1, WORK, LDX, ZERO, X, LDX )
470+
CALL CGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, CONE,
471+
$ U1, LDU1, WORK, LDX, CZERO, X, LDX )
472472
*
473473
DO I = 1, MIN(P,Q)-R
474-
X(I,I) = X(I,I) - ONE
474+
X(I,I) = X(I,I) - CONE
475475
END DO
476476
DO I = 1, R
477477
X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
478-
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - CMPLX( COS(THETA(I)),
479-
$ 0.0E0 )
478+
$ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) -
479+
$ CMPLX( COS(THETA(I)), ZERO )
480480
END DO
481481
*
482-
CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
483-
$ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
482+
CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q,
483+
$ CONE, X(P+1,1), LDX, V1T, LDV1T, CZERO, WORK, LDX )
484484
*
485485
CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
486-
$ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
486+
$ CONE, U2, LDU2, WORK, LDX, CZERO, X(P+1,1), LDX )
487487
*
488488
DO I = 1, MIN(M-P,Q)-R
489-
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
489+
X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - CONE
490490
END DO
491491
DO I = 1, R
492492
X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
493493
$ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
494-
$ CMPLX( SIN(THETA(R-I+1)), 0.0E0 )
494+
$ CMPLX( SIN(THETA(R-I+1)), ZERO )
495495
END DO
496496
*
497497
* Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
@@ -506,9 +506,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
506506
*
507507
* Compute I - U1'*U1
508508
*
509-
CALL CLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
510-
CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -REALONE,
511-
$ U1, LDU1, REALONE, WORK, LDU1 )
509+
CALL CLASET( 'Full', P, P, CZERO, CONE, WORK, LDU1 )
510+
CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE,
511+
$ U1, LDU1, ONE, WORK, LDU1 )
512512
*
513513
* Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
514514
*
@@ -517,9 +517,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
517517
*
518518
* Compute I - U2'*U2
519519
*
520-
CALL CLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
521-
CALL CHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -REALONE,
522-
$ U2, LDU2, REALONE, WORK, LDU2 )
520+
CALL CLASET( 'Full', M-P, M-P, CZERO, CONE, WORK, LDU2 )
521+
CALL CHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE,
522+
$ U2, LDU2, ONE, WORK, LDU2 )
523523
*
524524
* Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
525525
*
@@ -528,9 +528,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
528528
*
529529
* Compute I - V1T*V1T'
530530
*
531-
CALL CLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
532-
CALL CHERK( 'Upper', 'No transpose', Q, Q, -REALONE,
533-
$ V1T, LDV1T, REALONE, WORK, LDV1T )
531+
CALL CLASET( 'Full', Q, Q, CZERO, CONE, WORK, LDV1T )
532+
CALL CHERK( 'Upper', 'No transpose', Q, Q, -ONE,
533+
$ V1T, LDV1T, ONE, WORK, LDV1T )
534534
*
535535
* Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
536536
*
@@ -539,9 +539,9 @@ SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
539539
*
540540
* Check sorting
541541
*
542-
RESULT( 15 ) = REALZERO
542+
RESULT( 15 ) = ZERO
543543
DO I = 1, R
544-
IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
544+
IF( THETA(I).LT.ZERO .OR. THETA(I).GT.PIOVER2 ) THEN
545545
RESULT( 15 ) = ULPINV
546546
END IF
547547
IF( I.GT.1) THEN

0 commit comments

Comments
 (0)