Problem in new function ?SYTRS2 was found: wrong answer is got for a matrix A(LDA,N) and for B(LDB,N). If we take LDA=N, LDB=N then everything is OK.
The bug is in calling ?TRSM and ?SCAL functions:
CALL ?TRSM('L','L','T','U',N,NRHS,ONE,A,N,B,N)
CALL ?SCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
The calling should be
CALL ?TRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,N)
CALL ?SCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
The problem vanishes after these corrections

