Merge branch 'Reference-LAPACK:master' into xorgqr

Author: jprhyne

BLAS/SRC/CMakeLists.txt

@@ -32,15 +32,15 @@
 set(SBLAS1 isamax.f sasum.f saxpy.f scopy.f sdot.f snrm2.f90
 	srot.f srotg.f90 sscal.f sswap.f sdsdot.f srotmg.f srotm.f)
 
-set(CBLAS1 scabs1.f scasum.f scnrm2.f90 icamax.f caxpy.f ccopy.f
+set(CBLAS1 scabs1.f scasum.f scnrm2.f90 icamax.f90 caxpy.f ccopy.f
 	cdotc.f cdotu.f csscal.f crotg.f90 cscal.f cswap.f csrot.f)
 
 set(DBLAS1 idamax.f dasum.f daxpy.f dcopy.f ddot.f dnrm2.f90
 	drot.f drotg.f90 dscal.f dsdot.f dswap.f drotmg.f drotm.f)
 
 set(DB1AUX sscal.f isamax.f)
 
-set(ZBLAS1 dcabs1.f dzasum.f dznrm2.f90 izamax.f zaxpy.f zcopy.f
+set(ZBLAS1 dcabs1.f dzasum.f dznrm2.f90 izamax.f90 zaxpy.f zcopy.f
 	zdotc.f zdotu.f zdscal.f zrotg.f90 zscal.f zswap.f zdrot.f)
 
 set(CB1AUX
@@ -49,7 +49,7 @@ set(CB1AUX
     sswap.f)
 
 set(ZB1AUX
-    icamax.f idamax.f
+    icamax.f90 idamax.f
     cgemm.f cherk.f cscal.f ctrsm.f
     dasum.f daxpy.f dcopy.f ddot.f dgemm.f dgemv.f dnrm2.f90 drot.f dscal.f
     dswap.f

BLAS/SRC/icamax.f renamed to BLAS/SRC/DEPRECATED/icamax.f

@@ -35,6 +35,16 @@
 *>
 *> alpha and beta are scalars, and A, B and C are matrices, with op( A )
 *> an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
+*>
+*> Note: if alpha and/or beta is zero, some parts of the matrix-matrix
+*>  operations are not performed. This results in the following NaN/Inf
+*>  propagation quirks:
+*>
+*>  1. If alpha is zero, NaNs or Infs in A or B do not affect the result.
+*>  2. If both alpha and beta are zero, then a zero matrix is returned in C,
+*>   irrespective of any NaNs or Infs in A, B or C.
+*>  3. If only beta is zero, alpha*op( A )*op( B ) is returned, irrespective
+*>   of any NaNs or Infs in C.
 *> \endverbatim
 *
 *  Arguments:
@@ -92,7 +102,9 @@
 *> \param[in] ALPHA
 *> \verbatim
 *>          ALPHA is COMPLEX
-*>           On entry, ALPHA specifies the scalar alpha.
+*>           On entry, ALPHA specifies the scalar alpha. If ALPHA is zero the
+*>           values in A and B do not affect the result. This also means that
+*>           NaN/Inf propagation from A and B is inhibited if ALPHA is zero.
 *> \endverbatim
 *>
 *> \param[in] A
@@ -102,7 +114,10 @@
 *>           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
 *>           part of the array  A  must contain the matrix  A,  otherwise
 *>           the leading  k by m  part of the array  A  must contain  the
-*>           matrix A.
+*>           matrix A, except if ALPHA is zero.
+*>           If ALPHA is zero, none of the values in A affect the result, even
+*>           if they are NaN/Inf. This also implies that if ALPHA is zero,
+*>           the matrix elements of A need not be initialized by the caller.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -121,7 +136,10 @@
 *>           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
 *>           part of the array  B  must contain the matrix  B,  otherwise
 *>           the leading  n by k  part of the array  B  must contain  the
-*>           matrix B.
+*>           matrix B, except if ALPHA is zero.
+*>           If ALPHA is zero, none of the values in B affect the result, even
+*>           if they are NaN/Inf. This also implies that if ALPHA is zero,
+*>           the matrix elements of B need not be initialized by the caller.
 *> \endverbatim
 *>
 *> \param[in] LDB
@@ -136,16 +154,19 @@
 *> \param[in] BETA
 *> \verbatim
 *>          BETA is COMPLEX
-*>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
-*>           supplied as zero then C need not be set on input.
+*>           On entry,  BETA  specifies the scalar  beta. If BETA is zero the
+*>           values in C do not affect the result. This also means that
+*>           NaN/Inf propagation from C is inhibited if BETA is zero.
 *> \endverbatim
 *>
 *> \param[in,out] C
 *> \verbatim
 *>          C is COMPLEX array, dimension ( LDC, N )
 *>           Before entry, the leading  m by n  part of the array  C must
-*>           contain the matrix  C,  except when  beta  is zero, in which
-*>           case C need not be set on entry.
+*>           contain the matrix  C, except if beta is zero.
+*>           If beta is zero, none of the values in C affect the result, even
+*>           if they are NaN/Inf. This also implies that if beta is zero,
+*>           the matrix elements of C need not be initialized by the caller.
 *>           On exit, the array  C  is overwritten by the  m by n  matrix
 *>           ( alpha*op( A )*op( B ) + beta*C ).
 *> \endverbatim

BLAS/SRC/izamax.f renamed to BLAS/SRC/DEPRECATED/izamax.f

@@ -35,6 +35,16 @@
 *>
 *> alpha and beta are scalars, and A, B and C are matrices, with op( A )
 *> an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
+*>
+*> Note: if alpha and/or beta is zero, some parts of the matrix-matrix
+*>  operations are not performed. This results in the following NaN/Inf
+*>  propagation quirks:
+*>
+*>  1. If alpha is zero, NaNs or Infs in A or B do not affect the result.
+*>  2. If both alpha and beta are zero, then a zero matrix is returned in C,
+*>   irrespective of any NaNs or Infs in A, B or C.
+*>  3. If only beta is zero, alpha*op( A )*op( B ) is returned, irrespective
+*>   of any NaNs or Infs in C.
 *> \endverbatim
 *
 *  Arguments:
@@ -51,6 +61,9 @@
 *>              TRANSA = 'T' or 't',  op( A ) = A**T.
 *>
 *>              TRANSA = 'C' or 'c',  op( A ) = A**T.
+*>
+*>           Note: TRANSA = 'C' is supported for the sake of API consistency
+*>           between all ?GEMM variants.
 *> \endverbatim
 *>
 *> \param[in] TRANSB
@@ -64,6 +77,9 @@
 *>              TRANSB = 'T' or 't',  op( B ) = B**T.
 *>
 *>              TRANSB = 'C' or 'c',  op( B ) = B**T.
+*>
+*>           Note: TRANSB = 'C' is supported for the sake of API consistency
+*>           between all ?GEMM variants.
 *> \endverbatim
 *>
 *> \param[in] M
@@ -92,7 +108,9 @@
 *> \param[in] ALPHA
 *> \verbatim
 *>          ALPHA is DOUBLE PRECISION.
-*>           On entry, ALPHA specifies the scalar alpha.
+*>           On entry, ALPHA specifies the scalar alpha. If ALPHA is zero the
+*>           values in A and B do not affect the result. This also means that
+*>           NaN/Inf propagation from A and B is inhibited if ALPHA is zero.
 *> \endverbatim
 *>
 *> \param[in] A
@@ -102,7 +120,10 @@
 *>           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
 *>           part of the array  A  must contain the matrix  A,  otherwise
 *>           the leading  k by m  part of the array  A  must contain  the
-*>           matrix A.
+*>           matrix A, except if ALPHA is zero.
+*>           If ALPHA is zero, none of the values in A affect the result, even
+*>           if they are NaN/Inf. This also implies that if ALPHA is zero,
+*>           the matrix elements of A need not be initialized by the caller.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -121,7 +142,10 @@
 *>           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
 *>           part of the array  B  must contain the matrix  B,  otherwise
 *>           the leading  n by k  part of the array  B  must contain  the
-*>           matrix B.
+*>           matrix B, except if ALPHA is zero.
+*>           If ALPHA is zero, none of the values in B affect the result, even
+*>           if they are NaN/Inf. This also implies that if ALPHA is zero,
+*>           the matrix elements of B need not be initialized by the caller.
 *> \endverbatim
 *>
 *> \param[in] LDB
@@ -136,16 +160,19 @@
 *> \param[in] BETA
 *> \verbatim
 *>          BETA is DOUBLE PRECISION.
-*>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
-*>           supplied as zero then C need not be set on input.
+*>           On entry,  BETA  specifies the scalar  beta.  If BETA is zero the
+*>           values in C do not affect the result. This also means that
+*>           NaN/Inf propagation from C is inhibited if BETA is zero.
 *> \endverbatim
 *>
 *> \param[in,out] C
 *> \verbatim
 *>          C is DOUBLE PRECISION array, dimension ( LDC, N )
 *>           Before entry, the leading  m by n  part of the array  C must
-*>           contain the matrix  C,  except when  beta  is zero, in which
-*>           case C need not be set on entry.
+*>           contain the matrix  C, except if beta is zero.
+*>           If beta is zero, none of the values in C affect the result, even
+*>           if they are NaN/Inf. This also implies that if beta is zero,
+*>           the matrix elements of C need not be initialized by the caller.
 *>           On exit, the array  C  is overwritten by the  m by n  matrix
 *>           ( alpha*op( A )*op( B ) + beta*C ).
 *> \endverbatim

BLAS/SRC/cgemm.f

@@ -89,7 +89,7 @@ function DNRM2( n, x, incx )
    integer, parameter :: wp = kind(1.d0)
    real(wp) :: DNRM2
 !
-!  -- Reference BLAS level1 routine (version 3.9.1) --
+!  -- Reference BLAS level1 routine --
 !  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 !  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 !     March 2021

BLAS/SRC/dgemm.f

@@ -90,7 +90,7 @@ function DZNRM2( n, x, incx )
    integer, parameter :: wp = kind(1.d0)
    real(wp) :: DZNRM2
 !
-!  -- Reference BLAS level1 routine (version 3.9.1) --
+!  -- Reference BLAS level1 routine --
 !  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 !  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 !     March 2021