diff --git a/BLAS/SRC/CMakeLists.txt b/BLAS/SRC/CMakeLists.txt
index d585729e6b..96364ebc8a 100644
--- a/BLAS/SRC/CMakeLists.txt
+++ b/BLAS/SRC/CMakeLists.txt
@@ -105,3 +105,7 @@ set_target_properties(
   SOVERSION ${LAPACK_MAJOR_VERSION}
   )
 lapack_install_library(${BLASLIB})
+
+if( TEST_FORTRAN_COMPILER )
+  add_dependencies( ${BLASLIB} run_test_zcomplexabs run_test_zcomplexdiv run_test_zcomplexmult run_test_zminMax )
+endif()
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 07df064d23..61eb03bc3f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -41,6 +41,40 @@ if(_is_coverage_build)
   find_package(codecov)
 endif()
 
+# By default test Fortran compiler complex abs and complex division
+option(TEST_FORTRAN_COMPILER "Test Fortran compiler complex abs and complex division" ON)
+if( TEST_FORTRAN_COMPILER )
+
+  add_executable( test_zcomplexabs ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexabs.f )
+  add_custom_target( run_test_zcomplexabs
+    COMMAND test_zcomplexabs 2> test_zcomplexabs.err
+    WORKING_DIRECTORY ${LAPACK_BINARY_DIR}/INSTALL
+    COMMENT "Running test_zcomplexabs in ${LAPACK_BINARY_DIR}/INSTALL with stderr: test_zcomplexabs.err"
+    SOURCES ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexabs.f )
+
+  add_executable( test_zcomplexdiv ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexdiv.f )
+  add_custom_target( run_test_zcomplexdiv
+    COMMAND test_zcomplexdiv 2> test_zcomplexdiv.err
+    WORKING_DIRECTORY ${LAPACK_BINARY_DIR}/INSTALL
+    COMMENT "Running test_zcomplexdiv in ${LAPACK_BINARY_DIR}/INSTALL with stderr: test_zcomplexdiv.err"
+    SOURCES ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexdiv.f )
+
+  add_executable( test_zcomplexmult ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexmult.f )
+  add_custom_target( run_test_zcomplexmult
+    COMMAND test_zcomplexmult 2> test_zcomplexmult.err
+    WORKING_DIRECTORY ${LAPACK_BINARY_DIR}/INSTALL
+    COMMENT "Running test_zcomplexmult in ${LAPACK_BINARY_DIR}/INSTALL with stderr: test_zcomplexmult.err"
+    SOURCES ${LAPACK_SOURCE_DIR}/INSTALL/test_zcomplexmult.f )
+
+  add_executable( test_zminMax ${LAPACK_SOURCE_DIR}/INSTALL/test_zminMax.f )
+  add_custom_target( run_test_zminMax
+    COMMAND test_zminMax 2> test_zminMax.err
+    WORKING_DIRECTORY ${LAPACK_BINARY_DIR}/INSTALL
+    COMMENT "Running test_zminMax in ${LAPACK_BINARY_DIR}/INSTALL with stderr: test_zminMax.err"
+    SOURCES ${LAPACK_SOURCE_DIR}/INSTALL/test_zminMax.f )
+    
+endif()
+
 # By default static library
 option(BUILD_SHARED_LIBS "Build shared libraries" OFF)
 
diff --git a/INSTALL/CMakeLists.txt b/INSTALL/CMakeLists.txt
index b6f6e838a0..7a693f686e 100644
--- a/INSTALL/CMakeLists.txt
+++ b/INSTALL/CMakeLists.txt
@@ -7,3 +7,16 @@ add_executable(secondtst_INT_ETIME second_INT_ETIME.f secondtst.f)
 add_executable(secondtst_INT_CPU_TIME second_INT_CPU_TIME.f secondtst.f)
 add_executable(testieee tstiee.f)
 add_executable(testversion ilaver.f LAPACK_version.f)
+
+if( NOT TARGET test_zcomplexabs )
+  add_executable( test_zcomplexabs test_zcomplexabs.f )
+endif()
+if( NOT TARGET test_zcomplexdiv )
+  add_executable( test_zcomplexdiv test_zcomplexdiv.f )
+endif()
+if( NOT TARGET test_zcomplexmult )
+  add_executable( test_zcomplexmult test_zcomplexmult.f )
+endif()
+if( NOT TARGET test_zminMax )
+  add_executable( test_zminMax test_zminMax.f )
+endif()
diff --git a/INSTALL/Makefile b/INSTALL/Makefile
index b5bdc5d4dc..1af431b23a 100644
--- a/INSTALL/Makefile
+++ b/INSTALL/Makefile
@@ -1,8 +1,8 @@
 TOPSRCDIR = ..
 include $(TOPSRCDIR)/make.inc
 
-.PHONY: all testlsame testslamch testdlamch testsecond testdsecnd testieee testversion
-all: testlsame testslamch testdlamch testsecond testdsecnd testieee testversion
+.PHONY: all testlsame testslamch testdlamch testsecond testdsecnd testieee testversion test_zcomplexabs test_zcomplexdiv test_zcomplexmult test_zminMax
+all: testlsame testslamch testdlamch testsecond testdsecnd testieee testversion test_zcomplexabs test_zcomplexdiv test_zcomplexmult test_zminMax
 
 testlsame: lsame.o lsametst.o
 	$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
@@ -27,6 +27,18 @@ testieee: tstiee.o $(TOPSRCDIR)/SRC/ieeeck.o $(TOPSRCDIR)/SRC/ilaenv.o $(TOPSRCD
 testversion: ilaver.o LAPACK_version.o
 	$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
 
+test_zcomplexabs: test_zcomplexabs.o
+	$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
+
+test_zcomplexdiv: test_zcomplexdiv.o
+	$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
+
+test_zcomplexmult: test_zcomplexmult.o
+	$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
+
+test_zminMax: test_zminMax.o
+	$(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^
+
 .PHONY: run
 run: all
 	./testlsame
@@ -36,13 +48,17 @@ run: all
 	./testdsecnd
 	./testieee
 	./testversion
+	./test_zcomplexabs 2> test_zcomplexabs.err
+	./test_zcomplexdiv 2> test_zcomplexdiv.err
+	./test_zcomplexmult 2> test_zcomplexmult.err
+	./test_zminMax 2> test_zminMax.err
 
 .PHONY: clean cleanobj cleanexe cleantest
 clean: cleanobj cleanexe cleantest
 cleanobj:
 	rm -f *.o
 cleanexe:
-	rm -f test*
+	rm -f testlsame testslamch testdlamch testsecond testdsecnd testieee testversion test_zcomplexabs test_zcomplexdiv test_zcomplexmult test_zminMax
 cleantest:
 	rm -f core
 
diff --git a/INSTALL/test_zcomplexabs.f b/INSTALL/test_zcomplexabs.f
new file mode 100644
index 0000000000..8dedcf18dc
--- /dev/null
+++ b/INSTALL/test_zcomplexabs.f
@@ -0,0 +1,303 @@
+*> \brief zabs tests the robustness and precision of the intrinsic ABS for double complex
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \author Weslley S. Pereira, University of Colorado Denver, U.S.
+*
+*> \verbatim
+*>
+*> Real values for test:
+*> (1) x = 2**m, where m = MINEXPONENT-DIGITS, ..., MINEXPONENT-1. Stop on the first success.
+*>     Mind that not all platforms might implement subnormal numbers.
+*> (2) x = 2**m, where m = MINEXPONENT, ..., 0. Stop on the first success.
+*> (3) x = OV, where OV is the overflow threshold. OV^2 overflows but the norm is OV.
+*> (4) x = 2**m, where m = MAXEXPONENT-1, ..., 1. Stop on the first success.
+*>
+*> Tests:
+*> (a) y = x + 0 * I, |y| = x
+*> (b) y = 0 + x * I, |y| = x
+*> (c) y = (3/4)*x + x * I, |y| = (5/4)*x whenever (3/4)*x and (5/4)*x can be exactly stored
+*> (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2) whenever (1/2)*x can be exactly stored
+*>
+*> Special cases:
+*>
+*> (i) Inf propagation
+*>    (1) y = Inf + 0 * I, |y| is Inf.
+*>    (2) y =-Inf + 0 * I, |y| is Inf.
+*>    (3) y = 0 + Inf * I, |y| is Inf.
+*>    (4) y = 0 - Inf * I, |y| is Inf.
+*>    (5) y = Inf + Inf * I, |y| is Inf.
+*>
+*> (n) NaN propagation
+*>    (1) y = NaN + 0 * I, |y| is NaN.
+*>    (2) y = 0 + NaN * I, |y| is NaN.
+*>    (3) y = NaN + NaN * I, |y| is NaN.
+*>
+*> \endverbatim
+*
+*> \ingroup auxOTHERauxiliary
+*
+*  =====================================================================
+      program zabs
+*
+*  -- LAPACK test routine --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+
+*     ..
+*     .. Local parameters ..
+      logical           debug
+      parameter       ( debug = .false. )
+      integer           N, nNaN, nInf
+      parameter       ( N = 4, nNaN = 3, nInf = 5 )
+      double precision  threeFourth, fiveFourth, oneHalf
+      parameter       ( threeFourth = 3.0d0 / 4,
+     $                  fiveFourth = 5.0d0 / 4,
+     $                  oneHalf = 1.0d0 / 2 )
+*     ..
+*     .. Local Variables ..
+      integer           i, min, Max, m, subnormalTreatedAs0,
+     $                  caseAFails, caseBFails, caseCFails, caseDFails
+      double precision  X( N ), R, answerC,
+     $                  answerD, aInf, aNaN, relDiff, b,
+     $                  eps, blueMin, blueMax, Xj, stepX(N), limX(N)
+      double complex    Y, cInf( nInf ), cNaN( nNaN )
+*
+*     .. Intrinsic Functions ..
+      intrinsic         ABS, DBLE, RADIX, CEILING, TINY, DIGITS, SQRT,
+     $                  MAXEXPONENT, MINEXPONENT, FLOOR, HUGE, DCMPLX,
+     $                  EPSILON
+
+*
+*     .. Initialize error counts ..
+      subnormalTreatedAs0 = 0
+      caseAFails = 0
+      caseBFails = 0
+      caseCFails = 0
+      caseDFails = 0
+*
+*     .. Initialize machine constants ..
+      min = MINEXPONENT(0.0d0)
+      Max = MAXEXPONENT(0.0d0)
+      m = DIGITS(0.0d0)
+      b = DBLE(RADIX(0.0d0))
+      eps = EPSILON(0.0d0)
+      blueMin = b**CEILING( (min - 1) * 0.5d0 )
+      blueMax = b**FLOOR( (Max - m + 1) * 0.5d0 )
+*
+*     .. Vector X ..
+      X(1) = TINY(0.0d0) * b**( DBLE(1-m) )
+      X(2) = TINY(0.0d0)
+      X(3) = HUGE(0.0d0)
+      X(4) = b**( DBLE(Max-1) )
+*
+*     .. Then modify X using the step ..
+      stepX(1) = 2.0
+      stepX(2) = 2.0
+      stepX(3) = 0.0
+      stepX(4) = 0.5
+*
+*     .. Up to the value ..
+      limX(1) = X(2)
+      limX(2) = 1.0
+      limX(3) = 0.0
+      limX(4) = 2.0
+*
+*     .. Inf entries ..
+      aInf = X(3) * 2
+      cInf(1) = DCMPLX( aInf, 0.0d0 )
+      cInf(2) = DCMPLX(-aInf, 0.0d0 )
+      cInf(3) = DCMPLX( 0.0d0, aInf )
+      cInf(4) = DCMPLX( 0.0d0,-aInf )
+      cInf(5) = DCMPLX( aInf,  aInf )
+*
+*     .. NaN entries ..
+      aNaN = aInf / aInf
+      cNaN(1) = DCMPLX( aNaN, 0.0d0 )
+      cNaN(2) = DCMPLX( 0.0d0, aNaN )
+      cNaN(3) = DCMPLX( aNaN,  aNaN )
+
+*
+*     .. Tests ..
+*
+      if( debug ) then
+        print *, '# X :=', X
+        print *, '# Blue min constant :=', blueMin
+        print *, '# Blue max constant :=', blueMax
+      endif
+*
+      Xj = X(1)
+      if( Xj .eq. 0.0d0 ) then
+        subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+        if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+            print *, "!! fl( subnormal ) may be 0"
+        endif
+      else
+        do 100 i = 1, N
+            Xj = X(i)
+            if( Xj .eq. 0.0d0 ) then
+                subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+                if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                    print *, "!! fl( subnormal ) may be 0"
+                endif
+            endif
+ 100    continue
+      endif
+*
+*     Test (a) y = x + 0 * I, |y| = x
+      do 10 i = 1, N
+        Xj = X(i)
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [a] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                Y = DCMPLX( Xj, 0.0d0 )
+                R = ABS( Y )
+                if( R .ne. Xj ) then
+                    caseAFails = caseAFails + 1
+                    if( caseAFails .eq. 1 ) then
+                        print *, "!! Some ABS(x+0*I) differ from ABS(x)"
+                    endif
+                    WRITE( 0, FMT = 9999 ) 'a',i, Xj, '(1+0*I)', R, Xj
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif
+  10  continue
+*
+*     Test (b) y = 0 + x * I, |y| = x
+      do 20 i = 1, N
+        Xj = X(i)
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [b] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                Y = DCMPLX( 0.0d0, Xj )
+                R = ABS( Y )
+                if( R .ne. Xj ) then
+                    caseBFails = caseBFails + 1
+                    if( caseBFails .eq. 1 ) then
+                        print *, "!! Some ABS(0+x*I) differ from ABS(x)"
+                    endif
+                    WRITE( 0, FMT = 9999 ) 'b',i, Xj, '(0+1*I)', R, Xj
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif
+  20  continue
+*
+*     Test (c) y = (3/4)*x + x * I, |y| = (5/4)*x
+      do 30 i = 1, N
+        if( i .eq. 3 ) go to 30
+        if( i .eq. 1 ) then
+            Xj = 4*X(i)
+        else
+            Xj = X(i)
+        endif
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [c] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                answerC = fiveFourth * Xj
+                Y = DCMPLX( threeFourth * Xj, Xj )
+                R = ABS( Y )
+                if( R .ne. answerC ) then
+                    caseCFails = caseCFails + 1
+                    if( caseCFails .eq. 1 ) then
+                        print *, 
+     $              "!! Some ABS(x*(3/4+I)) differ from (5/4)*ABS(x)"
+                    endif
+                    WRITE( 0, FMT = 9999 ) 'c',i, Xj, '(3/4+I)', R,
+     $                                      answerC
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif
+  30  continue
+*
+*     Test (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2)
+      do 40 i = 1, N
+        if( i .eq. 1 ) then
+            Xj = 2*X(i)
+        else
+            Xj = X(i)
+        endif
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [d] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                answerD = (oneHalf * Xj) * SQRT(2.0d0)
+                if( answerD .eq. 0.0d0 ) then
+                    subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+                    if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                        print *, "!! [d] fl( subnormal ) may be 0"
+                    endif
+                else
+                    Y = DCMPLX( oneHalf * Xj, oneHalf * Xj )
+                    R = ABS( Y )
+                    relDiff = ABS(R-answerD)/answerD
+                    if( relDiff .ge. (0.5*eps) ) then
+                        caseDFails = caseDFails + 1
+                        if( caseDFails .eq. 1 ) then
+                            print *, 
+     $                "!! Some ABS(x*(1+I)) differ from sqrt(2)*ABS(x)"
+                        endif
+                        WRITE( 0, FMT = 9999 ) 'd',i, (oneHalf*Xj),
+     $                                  '(1+1*I)', R, answerD
+                    endif
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif 
+  40  continue
+*
+*     Test (e) Infs
+      do 50 i = 1, nInf
+        Y = cInf(i)
+        R = ABS( Y )
+        if( .not.(R .gt. HUGE(0.0d0)) ) then
+            WRITE( *, FMT = 9997 ) 'i',i, Y, R
+        endif
+  50  continue
+*
+*     Test (f) NaNs
+      do 60 i = 1, nNaN
+        Y = cNaN(i)
+        R = ABS( Y )
+        if( R .eq. R ) then
+            WRITE( *, FMT = 9998 ) 'n',i, Y, R
+        endif
+  60  continue
+*
+*     If anything was written to stderr, print the message
+      if( (caseAFails .gt. 0) .or. (caseBFails .gt. 0) .or.
+     $    (caseCFails .gt. 0) .or. (caseDFails .gt. 0) )
+     $      print *, "# Please check the failed ABS(a+b*I) in [stderr]"
+*
+*     .. Formats ..
+ 9997 FORMAT( '[',A1,I1, '] ABS(', (ES8.1,SP,ES8.1,"*I"), ' ) = ',
+     $        ES8.1, ' differs from Inf' )
+*
+ 9998 FORMAT( '[',A1,I1, '] ABS(', (ES8.1,SP,ES8.1,"*I"), ' ) = ',
+     $        ES8.1, ' differs from NaN' )
+*
+ 9999 FORMAT( '[',A1,I1, '] ABS(', ES24.16E3, ' * ', A7, ' ) = ',
+     $         ES24.16E3, ' differs from ', ES24.16E3 )
+*
+*     End of zabs
+*
+      END
diff --git a/INSTALL/test_zcomplexdiv.f b/INSTALL/test_zcomplexdiv.f
new file mode 100644
index 0000000000..fa3ea73b8e
--- /dev/null
+++ b/INSTALL/test_zcomplexdiv.f
@@ -0,0 +1,370 @@
+*> \brief zdiv tests the robustness and precision of the double complex division
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Authors:
+*  ========
+*
+*> \author Weslley S. Pereira, University of Colorado Denver, U.S.
+*
+*> \verbatim
+*>
+*> Real values for test:
+*> (1) x = 2**m, where m = MINEXPONENT-DIGITS, ..., MINEXPONENT-1.
+*>     Mind that not all platforms might implement subnormal numbers.
+*> (2) x = 2**m, where m = MINEXPONENT, ..., 0.
+*> (3) x = OV, where OV is the overflow threshold. OV^2 overflows but the norm is OV.
+*> (4) x = 2**m, where m = MAXEXPONENT-1, ..., 1.
+*>
+*> Tests:
+*> (a) y = x + 0 * I, y/y = 1
+*> (b) y = 0 + x * I, y/y = 1
+*> (c) y = x + x * I, y/y = 1
+*> (d) y1 = 0 + x * I, y2 = x + 0 * I, y1/y2 = I
+*> (e) y1 = 0 + x * I, y2 = x + 0 * I, y2/y1 = -I
+*> (f) y = x + x * I, y/conj(y) = I
+*>
+*> Special cases:
+*>
+*> (i) Inf inputs:
+*>    (1) y = ( Inf + 0   * I)
+*>    (2) y = ( 0   + Inf * I)
+*>    (3) y = (-Inf + 0   * I)
+*>    (4) y = ( 0   - Inf * I)
+*>    (5) y = ( Inf + Inf * I)
+*> Tests:
+*>    (a) 0 / y is either 0 or NaN.
+*>    (b) 1 / y is either 0 or NaN.
+*>    (c) y / y is NaN.
+*>
+*> (n) NaN inputs:
+*>    (1) y = (NaN + 0   * I)
+*>    (2) y = (0   + NaN * I)
+*>    (3) y = (NaN + NaN * I)
+*> Tests:
+*>    (a) 0 / y is NaN.
+*>    (b) 1 / y is NaN.
+*>    (c) y / y is NaN.
+*>
+*> \endverbatim
+*
+*> \ingroup auxOTHERauxiliary
+*
+*  =====================================================================
+      program zdiv
+*
+*  -- LAPACK test routine --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+
+*     ..
+*     .. Local parameters ..
+      logical           debug
+      parameter       ( debug = .false. )
+      integer           N, nNaN, nInf
+      parameter       ( N = 4, nNaN = 3, nInf = 5 )
+      double precision  threeFourth, fiveFourth
+      parameter       ( threeFourth = 3.0d0 / 4,
+     $                  fiveFourth = 5.0d0 / 4 )
+      double complex    czero, cone
+      parameter       ( czero = DCMPLX( 0.0d0, 0.0d0 ),
+     $                  cone  = DCMPLX( 1.0d0, 0.0d0 ) )
+*     ..
+*     .. Local Variables ..
+      integer           i, min, Max, m,
+     $                  subnormalTreatedAs0, caseAFails, caseBFails,
+     $                  caseCFails, caseDFails, caseEFails, caseFFails
+      double precision  X( N ), aInf, aNaN, b,
+     $                  eps, blueMin, blueMax, OV, Xj, stepX(N), limX(N)
+      double complex    Y, Y2, R, cInf( nInf ), cNaN( nNaN )
+*
+*     .. Intrinsic Functions ..
+      intrinsic         DCONJG, DBLE, RADIX, CEILING, TINY, DIGITS,
+     $                  MAXEXPONENT, MINEXPONENT, FLOOR, HUGE, DCMPLX,
+     $                  EPSILON
+
+*
+*     .. Initialize error counts ..
+      subnormalTreatedAs0 = 0
+      caseAFails = 0
+      caseBFails = 0
+      caseCFails = 0
+      caseDFails = 0
+      caseEFails = 0
+      caseFFails = 0
+*
+*     .. Initialize machine constants ..
+      min = MINEXPONENT(0.0d0)
+      Max = MAXEXPONENT(0.0d0)
+      m = DIGITS(0.0d0)
+      b = DBLE(RADIX(0.0d0))
+      eps = EPSILON(0.0d0)
+      blueMin = b**CEILING( (min - 1) * 0.5d0 )
+      blueMax = b**FLOOR( (Max - m + 1) * 0.5d0 )
+      OV = HUGE(0.0d0)
+*
+*     .. Vector X ..
+      X(1) = TINY(0.0d0) * b**( DBLE(1-m) )
+      X(2) = TINY(0.0d0)
+      X(3) = OV
+      X(4) = b**( DBLE(Max-1) )
+*
+*     .. Then modify X using the step ..
+      stepX(1) = 2.0
+      stepX(2) = 2.0
+      stepX(3) = 0.0
+      stepX(4) = 0.5
+*
+*     .. Up to the value ..
+      limX(1) = X(2)
+      limX(2) = 1.0
+      limX(3) = 0.0
+      limX(4) = 2.0
+*
+*     .. Inf entries ..
+      aInf = OV * 2
+      cInf(1) = DCMPLX( aInf, 0.0d0 )
+      cInf(2) = DCMPLX(-aInf, 0.0d0 )
+      cInf(3) = DCMPLX( 0.0d0, aInf )
+      cInf(4) = DCMPLX( 0.0d0,-aInf )
+      cInf(5) = DCMPLX( aInf,  aInf )
+*
+*     .. NaN entries ..
+      aNaN = aInf / aInf
+      cNaN(1) = DCMPLX( aNaN, 0.0d0 )
+      cNaN(2) = DCMPLX( 0.0d0, aNaN )
+      cNaN(3) = DCMPLX( aNaN,  aNaN )
+
+*
+*     .. Tests ..
+*
+      if( debug ) then
+        print *, '# X :=', X
+        print *, '# Blue min constant :=', blueMin
+        print *, '# Blue max constant :=', blueMax
+      endif
+*
+      Xj = X(1)
+      if( Xj .eq. 0.0d0 ) then
+        subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+        if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+            print *, "!! fl( subnormal ) may be 0"
+        endif
+      else
+        do 100 i = 1, N
+            Xj = X(i)
+            if( Xj .eq. 0.0d0 ) then
+                subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+                if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                    print *, "!! fl( subnormal ) may be 0"
+                endif
+            endif
+ 100    continue
+      endif
+*
+*     Test (a) y = x + 0 * I, y/y = 1
+      do 10 i = 1, N
+        Xj = X(i)
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [a] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                Y = DCMPLX( Xj, 0.0d0 )
+                R = Y / Y
+                if( R .ne. 1.0D0 ) then
+                    caseAFails = caseAFails + 1
+                    if( caseAFails .eq. 1 ) then
+                        print *, "!! Some (x+0*I)/(x+0*I) differ from 1"
+                    endif
+                    WRITE( 0, FMT = 9999 ) 'a',i, Xj,
+     $                                     '(x+0*I)/(x+0*I)', R, 1.0D0
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif
+  10  continue
+*
+*     Test (b) y = 0 + x * I, y/y = 1
+      do 20 i = 1, N
+        Xj = X(i)
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [b] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                Y = DCMPLX( 0.0d0, Xj )
+                R = Y / Y
+                if( R .ne. 1.0D0 ) then
+                    caseBFails = caseBFails + 1
+                    if( caseBFails .eq. 1 ) then
+                        print *, "!! Some (0+x*I)/(0+x*I) differ from 1"
+                    endif
+                    WRITE( 0, FMT = 9999 ) 'b',i, Xj,
+     $                                     '(0+x*I)/(0+x*I)', R, 1.0D0
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif
+  20  continue
+*
+*     Test (c) y = x + x * I, y/y = 1
+      do 30 i = 1, N
+        Xj = X(i)
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [c] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                Y = DCMPLX( Xj, Xj )
+                R = Y / Y
+                if( R .ne. 1.0D0 ) then
+                    caseCFails = caseCFails + 1
+                    if( caseCFails .eq. 1 ) then
+                        print *, "!! Some (x+x*I)/(x+x*I) differ from 1"
+                    endif
+                    WRITE( 0, FMT = 9999 ) 'c',i, Xj,
+     $                                     '(x+x*I)/(x+x*I)', R, 1.0D0
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif
+  30  continue
+*
+*     Test (d) y1 = 0 + x * I, y2 = x + 0 * I, y1/y2 = I
+      do 40 i = 1, N
+        Xj = X(i)
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [d] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                Y  = DCMPLX( 0.0d0, Xj )
+                Y2 = DCMPLX( Xj, 0.0d0 )
+                R = Y / Y2
+                if( R .ne. DCMPLX(0.0D0,1.0D0) ) then
+                    caseDFails = caseDFails + 1
+                    if( caseDFails .eq. 1 ) then
+                        print *, "!! Some (0+x*I)/(x+0*I) differ from I"
+                    endif
+                    WRITE( 0, FMT = 9999 ) 'd',i, Xj, '(0+x*I)/(x+0*I)',
+     $                                      R, DCMPLX(0.0D0,1.0D0)
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif
+  40  continue
+*
+*     Test (e) y1 = 0 + x * I, y2 = x + 0 * I, y2/y1 = -I
+      do 50 i = 1, N
+        Xj = X(i)
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [e] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                Y  = DCMPLX( 0.0d0, Xj )
+                Y2 = DCMPLX( Xj, 0.0d0 )
+                R = Y2 / Y
+                if( R .ne. DCMPLX(0.0D0,-1.0D0) ) then
+                    caseEFails = caseEFails + 1
+                    if( caseEFails .eq. 1 ) then
+                        print *,"!! Some (x+0*I)/(0+x*I) differ from -I"
+                    endif
+                    WRITE( 0, FMT = 9999 ) 'e',i, Xj, '(x+0*I)/(0+x*I)',
+     $                                      R, DCMPLX(0.0D0,-1.0D0)
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif
+  50  continue
+*
+*     Test (f) y = x + x * I, y/conj(y) = I
+      do 60 i = 1, N
+        Xj = X(i)
+        if( Xj .eq. 0.0d0 ) then
+            subnormalTreatedAs0 = subnormalTreatedAs0 + 1
+            if( debug .or. subnormalTreatedAs0 .eq. 1 ) then
+                print *, "!! [f] fl( subnormal ) may be 0"
+            endif
+        else
+            do while( Xj .ne. limX(i) )
+                Y  = DCMPLX( Xj, Xj )
+                R = Y / DCONJG( Y )
+                if( R .ne. DCMPLX(0.0D0,1.0D0) ) then
+                    caseFFails = caseFFails + 1
+                    if( caseFFails .eq. 1 ) then
+                        print *, "!! Some (x+x*I)/(x-x*I) differ from I"
+                    endif
+                    WRITE( 0, FMT = 9999 ) 'f',i, Xj, '(x+x*I)/(x-x*I)',
+     $                                      R, DCMPLX(0.0D0,1.0D0)
+                endif
+                Xj = Xj * stepX(i)
+            end do
+        endif
+  60  continue
+*
+*     Test (g) Infs
+      do 70 i = 1, nInf
+          Y = cInf(i)
+          R = czero / Y
+          if( (R .ne. czero) .and. (R .eq. R) ) then
+              WRITE( *, FMT = 9998 ) 'ia',i, czero, Y, R, 'NaN and 0'
+          endif
+          R = cone / Y
+          if( (R .ne. czero) .and. (R .eq. R) ) then
+              WRITE( *, FMT = 9998 ) 'ib',i, cone, Y, R, 'NaN and 0'
+          endif
+          R = Y / Y
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'ic',i, Y, Y, R, 'NaN'
+          endif
+  70  continue
+*
+*     Test (h) NaNs
+      do 80 i = 1, nNaN
+          Y = cNaN(i)
+          R = czero / Y
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'na',i, czero, Y, R, 'NaN'
+          endif
+          R = cone / Y
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'nb',i, cone, Y, R, 'NaN'
+          endif
+          R = Y / Y
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'nc',i, Y, Y, R, 'NaN'
+          endif
+  80  continue
+*
+*     If anything was written to stderr, print the message
+      if( (caseAFails .gt. 0) .or. (caseBFails .gt. 0) .or.
+     $    (caseCFails .gt. 0) .or. (caseDFails .gt. 0) .or.
+     $    (caseEFails .gt. 0) .or. (caseFFails .gt. 0) )
+     $      print *, "# Please check the failed divisions in [stderr]"
+*
+*     .. Formats ..
+ 9998 FORMAT( '[',A2,I1, '] ', (ES24.16E3,SP,ES24.16E3,"*I"), ' * ',
+     $         (ES24.16E3,SP,ES24.16E3,"*I"), ' = ',
+     $         (ES24.16E3,SP,ES24.16E3,"*I"), ' differs from ', A10 )
+*
+ 9999 FORMAT( '[',A2,I1, '] X = ', ES24.16E3, ' : ', A15, ' = ',
+     $         (ES24.16E3,SP,ES24.16E3,"*I"), ' differs from ',
+     $         (ES24.16E3,SP,ES24.16E3,"*I") )
+*
+*     End of zdiv
+*
+      END
\ No newline at end of file
diff --git a/INSTALL/test_zcomplexmult.f b/INSTALL/test_zcomplexmult.f
new file mode 100644
index 0000000000..6fce1ef8fc
--- /dev/null
+++ b/INSTALL/test_zcomplexmult.f
@@ -0,0 +1,136 @@
+*> \brief zmul tests the robustness and precision of the double complex multiplication
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Authors:
+*  ========
+*
+*> \author Weslley S. Pereira, University of Colorado Denver, U.S.
+*
+*> \verbatim
+*>
+*> Tests:
+*>
+*> (a) Inf inputs:
+*>    (1) y = ( Inf + 0   * I)
+*>    (2) y = (-Inf + 0   * I)
+*>    (3) y = ( 0   + Inf * I)
+*>    (4) y = ( 0   - Inf * I)
+*>    (5) y = ( Inf + Inf * I)
+*> Tests:
+*>    (a) 0 * y is NaN.
+*>    (b) 1 * y is y is either y or NaN.
+*>    (c) y * y is either  Inf or NaN (cases 1 and 3),
+*>                 either -Inf or NaN (cases 2 and 4),
+*>                 NaN (case 5).
+*>
+*> (b) NaN inputs:
+*>    (1) y = (NaN + 0   * I)
+*>    (2) y = (0   + NaN * I)
+*>    (3) y = (NaN + NaN * I)
+*> Tests:
+*>    (a) 0 * y is NaN.
+*>    (b) 1 * y is NaN.
+*>    (c) y * y is NaN.
+*>
+*> \endverbatim
+*
+*> \ingroup auxOTHERauxiliary
+*
+*  =====================================================================
+      program zmul
+*
+*  -- LAPACK test routine --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+
+*     ..
+*     .. Constants ..
+      integer           nNaN, nInf
+      parameter       ( nNaN = 3, nInf = 5 )
+      double complex    czero, cone
+      parameter       ( czero = DCMPLX( 0.0d0, 0.0d0 ),
+     $                  cone  = DCMPLX( 1.0d0, 0.0d0 ) )
+*     ..
+*     .. Local Variables ..
+      integer           i
+      double precision  aInf, aNaN, OV
+      double complex    Y, R, cInf( nInf ), cNaN( nNaN )
+*
+*     .. Intrinsic Functions ..
+      intrinsic         HUGE, DCMPLX
+
+*
+*     .. Inf entries ..
+      OV = HUGE(0.0d0)
+      aInf = OV * 2
+      cInf(1) = DCMPLX( aInf, 0.0d0 )
+      cInf(2) = DCMPLX(-aInf, 0.0d0 )
+      cInf(3) = DCMPLX( 0.0d0, aInf )
+      cInf(4) = DCMPLX( 0.0d0,-aInf )
+      cInf(5) = DCMPLX( aInf,  aInf )
+*
+*     .. NaN entries ..
+      aNaN = aInf / aInf
+      cNaN(1) = DCMPLX( aNaN, 0.0d0 )
+      cNaN(2) = DCMPLX( 0.0d0, aNaN )
+      cNaN(3) = DCMPLX( aNaN,  aNaN )
+
+*
+*     .. Tests ..
+*
+*     Test (a) Infs
+      do 10 i = 1, nInf
+          Y = cInf(i)
+          R = czero * Y
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'ia',i, czero, Y, R, 'NaN'
+          endif
+          R = cone * Y
+          if( (R .ne. Y) .and. (R .eq. R) ) then
+              WRITE( *, FMT = 9998 ) 'ib',i, cone, Y, R,
+     $                               'the input and NaN'
+          endif
+          R = Y * Y
+          if( (i.eq.1) .or. (i.eq.2) ) then
+              if( (R .ne. cInf(1)) .and. (R .eq. R) ) then
+                  WRITE( *, FMT = 9998 ) 'ic',i, Y, Y, R, 'Inf and NaN'
+              endif
+          else if( (i.eq.3) .or. (i.eq.4) ) then
+              if( (R .ne. cInf(2)) .and. (R .eq. R) ) then
+                  WRITE( *, FMT = 9998 ) 'ic',i, Y, Y, R, '-Inf and NaN'
+              endif
+          else 
+              if( R .eq. R ) then
+                  WRITE( *, FMT = 9998 ) 'ic',i, Y, Y, R, 'NaN'
+              endif
+          endif
+  10  continue
+*
+*     Test (b) NaNs
+      do 20 i = 1, nNaN
+          Y = cNaN(i)
+          R = czero * Y
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'na',i, czero, Y, R, 'NaN'
+          endif
+          R = cone * Y
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'nb',i, cone, Y, R, 'NaN'
+          endif
+          R = Y * Y
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'nc',i, Y, Y, R, 'NaN'
+          endif
+  20  continue
+*
+*     .. Formats ..
+ 9998 FORMAT( '[',A2,I1, '] (', (ES24.16E3,SP,ES24.16E3,"*I"), ') * (',
+     $         (ES24.16E3,SP,ES24.16E3,"*I"), ') = (',
+     $         (ES24.16E3,SP,ES24.16E3,"*I"), ') differs from ', A17 )
+*
+*     End of zmul
+*
+      END
\ No newline at end of file
diff --git a/INSTALL/test_zminMax.f b/INSTALL/test_zminMax.f
new file mode 100644
index 0000000000..b98dd3c34f
--- /dev/null
+++ b/INSTALL/test_zminMax.f
@@ -0,0 +1,100 @@
+*> \brief zminMax tests the robustness and precision of the double-valued intrinsic operators MIN and MAX
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Authors:
+*  ========
+*
+*> \author Weslley S. Pereira, University of Colorado Denver, U.S.
+*
+*> \verbatim
+*>
+*> Tests with pairs of numbers (x,y):
+*> Inf inputs where x < y:
+*>    (1) (-Inf,   0)
+*>    (2) ( 0  , Inf)
+*>    (3) (-Inf, Inf)
+*> Inf inputs where x > y:
+*>    (4) ( 0  ,-Inf)
+*>    (5) ( Inf,   0)
+*>    (6) ( Inf,-Inf)
+*> NaN inputs to test NaN propagation:
+*>    (7) ( 0  , NaN)
+*>    (8) ( NaN,   0)
+*> The program tests MIN(x,y) and MAX(x,y) for every pair
+*>
+*> \endverbatim
+*
+*> \ingroup auxOTHERauxiliary
+*
+*  =====================================================================
+      program zmul
+*
+*  -- LAPACK test routine --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+
+*     ..
+*     .. Parameters ..
+      integer           n
+      parameter       ( n = 8 )
+      double precision  zero
+      parameter       ( zero = 0.0d0 )
+*     ..
+*     .. Local Variables ..
+      integer           i
+      double precision  aInf, aNaN, OV, R, X(n), Y(n)
+*
+*     .. Intrinsic Functions ..
+      intrinsic         HUGE, MIN, MAX
+
+*
+*     .. Inf and NaN entries ..
+      OV = HUGE(0.0d0)
+      aInf = OV * 2
+      aNaN = aInf / aInf
+      X = (/ -aInf, zero, -aInf,  zero, aInf,  aInf, zero, aNaN /)
+      Y = (/  zero, aInf,  aInf, -aInf, zero, -aInf, aNaN, zero /)
+
+*
+*     .. Tests ..
+*
+      do 10 i = 1, 3
+          R = MIN( X(i), Y(i) )
+          if( R .ne. X(i) ) then
+              WRITE( *, FMT = 9998 ) 'i',i, 'MIN', X(i), Y(i), R
+          endif
+          R = MAX( X(i), Y(i) )
+          if( R .ne. Y(i) ) then
+              WRITE( *, FMT = 9998 ) 'i',i, 'MAX', X(i), Y(i), R
+          endif
+  10  continue
+      do 20 i = 4, 6
+          R = MIN( X(i), Y(i) )
+          if( R .ne. Y(i) ) then
+              WRITE( *, FMT = 9998 ) 'i',i, 'MIN', X(i), Y(i), R
+          endif
+          R = MAX( X(i), Y(i) )
+          if( R .ne. X(i) ) then
+              WRITE( *, FMT = 9998 ) 'i',i, 'MAX', X(i), Y(i), R
+          endif
+  20  continue
+      do 30 i = 7, 8
+          R = MIN( X(i), Y(i) )
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'i',i, 'MIN', X(i), Y(i), R
+          endif
+          R = MAX( X(i), Y(i) )
+          if( R .eq. R ) then
+              WRITE( *, FMT = 9998 ) 'i',i, 'MAX', X(i), Y(i), R
+          endif
+  30  continue
+*
+*     .. Formats ..
+ 9998 FORMAT( '[',A1,I1, '] ', A3, '(', F5.0, ',', F5.0, ') = ', F5.0 )
+*
+*     End of zmul
+*
+      END
\ No newline at end of file
diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt
index ada0325f93..8cb3de29bc 100644
--- a/SRC/CMakeLists.txt
+++ b/SRC/CMakeLists.txt
@@ -513,6 +513,10 @@ set_target_properties(
   SOVERSION ${LAPACK_MAJOR_VERSION}
   )
 
+if( TEST_FORTRAN_COMPILER )
+  add_dependencies( ${LAPACKLIB} run_test_zcomplexabs run_test_zcomplexdiv run_test_zcomplexmult run_test_zminMax )
+endif()
+
 if(USE_XBLAS)
   target_link_libraries(${LAPACKLIB} PRIVATE ${XBLAS_LIBRARY})
 endif()