Fix sorcsd

[weslleyspereira]

Closes #695.

Revert "SORCSD2BY1: remove dead code"
This reverts commit d245b4f.

Revert "SORCSD: fix documentation on matrix dimensions"
This reverts commit bdcd890.

Keep corrections in the documentation, which make sense.

[weslleyspereira]

The double precision tests still show: 39 failing tests. See https://github.com/Reference-LAPACK/lapack/actions/runs/2806513051

My personal machine found 16 tests failing for the CS decomposition in double precision:

Command: "/usr/bin/cmake" "-DTEST=/home/weslleyp/storage/lapack/build_RelWithDebInfo/bin/xeigtstd" "-DINPUT=/home/weslleyp/storage/lapack/TESTING/csd.in" "-DOUTPUT=/home/weslleyp/storage/lapack/build_RelWithDebInfo/TESTING/dcsd.out" "-DINTDIR=." "-P" "/home/weslleyp/storage/lapack/TESTING/runtest.cmake"
Directory: /home/weslleyp/storage/lapack/build_RelWithDebInfo/TESTING
"LAPACK-xeigtstd_csd_in" start time: Aug 05 17:04 MDT
Output:
----------------------------------------------------------
Running: /home/weslleyp/storage/lapack/build_RelWithDebInfo/bin/xeigtstd
ARGS= OUTPUT_FILE;/home/weslleyp/storage/lapack/build_RelWithDebInfo/TESTING/dcsd.out;ERROR_FILE;/home/weslleyp/storage/lapack/build_RelWithDebInfo/TESTING/dcsd.out.err;INPUT_FILE;/home/weslleyp/storage/lapack/TESTING/csd.in
Test OUTPUT:

 Tests of the CS Decomposition routines

 LAPACK VERSION 3.*.1

 The following parameter values will be used:
    M:         0    10    10    10    10    21    24    30    22    32

    P:         0     4     4     0    10     9    10    20    12    12

    N:         0     0    10     4     4    15    12     8    20     8


 Relative machine underflow is taken to be    0.222507-307
 Relative machine overflow  is taken to be    0.179769+309
 Relative machine precision is taken to be    0.111022D-15

 Routines pass computational tests if test ratio is less than   30.00


 CSD routines passed the tests of the error exits (  8 tests done)

 CSD: CS Decomposition
 Matrix types: 
    1: Random orthogonal matrix (Haar measure)
    2: Nearly orthogonal matrix with uniformly distributed angles atan2( S, C ) in CS decomposition
    3: Random orthogonal matrix with clustered angles atan2( S, C ) in CS decomposition
 Test ratios: 
   2-by-2 CSD
    1: norm( U1' * X11 * V1 - C ) / ( max(  P,  Q) * max(norm(I-X'*X),EPS) )
    2: norm( U1' * X12 * V2-(-S)) / ( max(  P,M-Q) * max(norm(I-X'*X),EPS) )
    3: norm( U2' * X21 * V1 - S ) / ( max(M-P,  Q) * max(norm(I-X'*X),EPS) )
    4: norm( U2' * X22 * V2 - C ) / ( max(M-P,M-Q) * max(norm(I-X'*X),EPS) )
    5: norm( I - U1'*U1 ) / (   P   * EPS )
    6: norm( I - U2'*U2 ) / ( (M-P) * EPS )
    7: norm( I - V1'*V1 ) / (   Q   * EPS )
    8: norm( I - V2'*V2 ) / ( (M-Q) * EPS )
    9: principal angle ordering ( 0 or ULP )
   2-by-1 CSD
   10: norm( U1' * X11 * V1 - C ) / ( max(  P,  Q) * max(norm(I-X'*X),EPS) )
   11: norm( U2' * X21 * V1 - S ) / ( max(  M-P,Q) * max(norm(I-X'*X),EPS) )
   12: norm( I - U1'*U1 ) / (   P   * EPS )
   13: norm( I - U2'*U2 ) / ( (M-P) * EPS )
   14: norm( I - V1'*V1 ) / (   Q   * EPS )
   15: principal angle ordering ( 0 or ULP )
 M=  21 P=   9, Q=  15, type  2, test 10, ratio= 0.802656E+11
 M=  21 P=   9, Q=  15, type  2, test 11, ratio= 0.988228E+11
 M=  21 P=   9, Q=  15, type  4, test 11, ratio= 0.600480E+15
 M=  24 P=  10, Q=  12, type  4, test 10, ratio= 0.375300E+15
 M=  24 P=  10, Q=  12, type  4, test 11, ratio= 0.643371E+15
 M=  30 P=  20, Q=   8, type  2, test 10, ratio= 0.280538E+12
 M=  30 P=  20, Q=   8, type  2, test 11, ratio= 0.267436E+12
 M=  22 P=  12, Q=  20, type  1, test 10, ratio= 0.217202E+14
 M=  22 P=  12, Q=  20, type  1, test 11, ratio= 0.260562E+14
 M=  22 P=  12, Q=  20, type  4, test 11, ratio= 0.450360E+15
 M=  32 P=  12, Q=   8, type  1, test 10, ratio= 0.896914E+15
 M=  32 P=  12, Q=   8, type  1, test 11, ratio= 0.491797E+15
 M=  32 P=  12, Q=   8, type  3, test 10, ratio= 0.211857E+15
 M=  32 P=  12, Q=   8, type  3, test 11, ratio= 0.171219E+14
 M=  32 P=  12, Q=   8, type  4, test 10, ratio= 0.375300E+15
 M=  32 P=  12, Q=   8, type  4, test 11, ratio= 0.225180E+15
 CSD:     16 out of    600 tests failed to pass the threshold


 End of tests
 Total time used =         0.01 seconds


Test ERROR:

Test /home/weslleyp/storage/lapack/build_RelWithDebInfo/bin/xeigtstd returned 0

Any idea where they could be, @christoph-conrads ?

• Maybe related to that: I found an inconsistency between single and double precision codes, and proposed deb3d85 as a fix. I had the same failing tests.

[martin-frbg]

my tests in #695 showed you'd need to restore dorbdb6.f to fix the 39 errors.

[weslleyspereira]

my tests in #697 showed you'd need to restore dorbdb6.f to fix the 39 errors.

Thanks, @martin-frbg! I completely bypassed your comment: #695 (comment). I will apply those changes too.

[martin-frbg]

Seems what breaks the tests is the changes around line 266 (257 in the 3.10.1 version) where the conditional for the early RETURN was changed and elements set to zero that were previously left alone. (Probably the latter is the actual cause, still testing right now)

[martin-frbg]

Yes minimal fix is just

       IF( NORM_NEW .LE. N * EPS * NORM ) THEN
-         DO IX = 1, 1 + (M1-1)*INCX1, INCX1
-           X1( IX ) = ZERO
-         END DO
-         DO IX = 1, 1 + (M2-1)*INCX2, INCX2
-           X2( IX ) = ZERO
-         END DO
          RETURN
       END IF

but I do not see any benefit from any of the other changes in dorbdb6.f - maybe I'd need to be a mathematician rather than a lowly computational chemist for that ?

[weslleyspereira]

@martin-frbg, your minimal fix makes my tests pass indeed. Thanks!

However, I understand that setting X1 and X2 to this code is important for dorbdb5. See:

lapack/SRC/dorbdb5.f

Lines 218 to 226 in 3381a0e

	CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, LDQ2,
	$ WORK, LWORK, CHILDINFO )
	*
	* If the projection is nonzero, then return
	*
	IF( DNRM2(M1,X1,INCX1) .NE. ZERO
	$ .OR. DNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
	RETURN
	END IF

The point here, as I understand, was to identify zero projections. If a zero projection is found, use a different projection:

lapack/SRC/dorbdb5.f

Lines 228 to 264 in 3381a0e

	* Project each standard basis vector e_1,...,e_M1 in turn, stopping
	* when a nonzero projection is found
	*
	DO I = 1, M1
	DO J = 1, M1
	X1(J) = ZERO
	END DO
	X1(I) = ONE
	DO J = 1, M2
	X2(J) = ZERO
	END DO
	CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
	$ LDQ2, WORK, LWORK, CHILDINFO )
	IF( DNRM2(M1,X1,INCX1) .NE. ZERO
	$ .OR. DNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
	RETURN
	END IF
	END DO
	*
	* Project each standard basis vector e_(M1+1),...,e_(M1+M2) in turn,
	* stopping when a nonzero projection is found
	*
	DO I = 1, M2
	DO J = 1, M1
	X1(J) = ZERO
	END DO
	DO J = 1, M2
	X2(J) = ZERO
	END DO
	X2(I) = ONE
	CALL DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
	$ LDQ2, WORK, LWORK, CHILDINFO )
	IF( DNRM2(M1,X1,INCX1) .NE. ZERO
	$ .OR. DNRM2(M2,X2,INCX2) .NE. ZERO ) THEN
	RETURN
	END IF
	END DO

I will let @christoph-conrads comment on that before making additional changes.

[christoph-conrads]

The reason for modifying the xORBDB6 return value can be found in #634. xORBDB5 is calling xORBDB6 to determine if a vector x lies in the range of a matrix A or not. xORBDB6 is a procedure (a function with return type void in C). Hence, the only way to determine if x lies in ran(A) was to compute the norm of x after xORBDB6 returned and check how close it is to zero. xORBDB5 and xORBDB6 used different criteria for the zero check leading to a bug in xORCSD2BY1. To fix the issue and after consulting with the original author Brian D. Sutton, it was decided to set x explicitly to zero to save the caller from having to re-implement parts of the xORBDB6 logic.

Assuming setting x explicitly to zero is a problem, one could turn xORBDB6 into a function returning a boolean. Modifying the entire vector x and then computing its norm only to figure out the result of xORBDB6 is a tad bit too complicated anyway but this solution did not break the API.

[weslleyspereira]

I have rebased with the current master ( that includes #702 ), and then all tests pass in my machine. I am now running the tests on Github Actions. Thanks @mjacobse for noticing the issue with the types!

[weslleyspereira]

I will merge this PR since the failing tests in double precision were fixed by #702.

[Edit:] Actually, can someone review my changes so I can merge it? Thanks!

[angsch]

@weslleyspereira Thanks for fixing the issues. The changes make sense and I support merging your PR in.