re PR fortran/29383 (Fortran 2003/F95[TR15580:1999]: Floating point exception (IEEE) support)

	PR fortran/29383

gcc/fortran/
	* gfortran.h (gfc_simplify_ieee_selected_real_kind): New prototype.
	* libgfortran.h (GFC_FPE_*): Use simple integer values, valid in
	both C and Fortran.
	* expr.c (gfc_check_init_expr): Simplify IEEE_SELECTED_REAL_KIND.
	* simplify.c (gfc_simplify_ieee_selected_real_kind): New function.
	* module.c (mio_symbol): Keep track of symbols which came from
	intrinsic modules.
	(gfc_use_module): Keep track of the IEEE modules.
	* trans-decl.c (gfc_get_symbol_decl): Adjust code since
	we have new intrinsic modules.
	(gfc_build_builtin_function_decls): Build decls for
	ieee_procedure_entry and ieee_procedure_exit.
	(is_from_ieee_module, is_ieee_module_used, save_fp_state,
	restore_fp_state): New functions.
	(gfc_generate_function_code): Save and restore floating-point
	state on procedure entry/exit, when IEEE modules are used.
	* intrinsic.texi: Document the IEEE modules.

libgfortran/
	* configure.host: Add checks for IEEE support, rework priorities.
	* configure.ac: Define IEEE_SUPPORT, check for fpsetsticky and
	fpresetsticky.
	* configure: Regenerate.
	* Makefile.am: Build new ieee files, install IEEE_* modules.
	* Makefile.in: Regenerate.
	* gfortran.map (GFORTRAN_1.6): Add new symbols.
	* libgfortran.h (get_fpu_trap_exceptions, set_fpu_trap_exceptions,
	support_fpu_trap, set_fpu_except_flags, support_fpu_flag,
	support_fpu_rounding_mode, get_fpu_state, set_fpu_state): New
	prototypes.
	* config/fpu-*.h (get_fpu_trap_exceptions,
	set_fpu_trap_exceptions, support_fpu_trap, set_fpu_except_flags,
	support_fpu_flag, support_fpu_rounding_mode, get_fpu_state,
	set_fpu_state): New functions.
	* ieee/ieee_features.F90: New file.
	* ieee/ieee_exceptions.F90: New file.
	* ieee/ieee_arithmetic.F90: New file.
	* ieee/ieee_helper.c: New file.

gcc/testsuite/
	* lib/target-supports.exp (check_effective_target_fortran_ieee): 
	New function.
	* gfortran.dg/ieee/ieee.exp: New file.
	* gfortran.dg/ieee/ieee_1.F90: New file.
	* gfortran.dg/ieee/ieee_2.f90: New file.
	* gfortran.dg/ieee/ieee_3.f90: New file.
	* gfortran.dg/ieee/ieee_4.f90: New file.
	* gfortran.dg/ieee/ieee_5.f90: New file.
	* gfortran.dg/ieee/ieee_6.f90: New file.
	* gfortran.dg/ieee/ieee_7.f90: New file.
	* gfortran.dg/ieee/ieee_rounding_1.f90: New file.

From-SVN: r212102

19 files changed, 1600 insertions, 16 deletions

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index f1ac532..a5f6f9d5 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,3 +1,24 @@
+2014-06-28  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
+
+	PR fortran/29383
+	* gfortran.h (gfc_simplify_ieee_selected_real_kind): New prototype.
+	* libgfortran.h (GFC_FPE_*): Use simple integer values, valid in
+	both C and Fortran.
+	* expr.c (gfc_check_init_expr): Simplify IEEE_SELECTED_REAL_KIND.
+	* simplify.c (gfc_simplify_ieee_selected_real_kind): New function.
+	* module.c (mio_symbol): Keep track of symbols which came from
+	intrinsic modules.
+	(gfc_use_module): Keep track of the IEEE modules.
+	* trans-decl.c (gfc_get_symbol_decl): Adjust code since
+	we have new intrinsic modules.
+	(gfc_build_builtin_function_decls): Build decls for
+	ieee_procedure_entry and ieee_procedure_exit.
+	(is_from_ieee_module, is_ieee_module_used, save_fp_state,
+	restore_fp_state): New functions.
+	(gfc_generate_function_code): Save and restore floating-point
+	state on procedure entry/exit, when IEEE modules are used.
+	* intrinsic.texi: Document the IEEE modules.
+
 2014-06-25  Tobias Burnus  <burnus@net-b.de>
 
 	* interface.c (check_intents): Fix diagnostic with
diff --git a/gcc/fortran/expr.c b/gcc/fortran/expr.c
index feb089e..3e3a664 100644
--- a/gcc/fortran/expr.c
+++ b/gcc/fortran/expr.c
@@ -2460,9 +2460,23 @@ gfc_check_init_expr (gfc_expr *e)
 
       {
 	gfc_intrinsic_sym* isym;
-	gfc_symbol* sym;
+	gfc_symbol* sym = e->symtree->n.sym;
+
+	/* Special case for IEEE_SELECTED_REAL_KIND from the intrinsic
+	   module IEEE_ARITHMETIC, which is allowed in initialization
+	   expressions.  */
+	if (!strcmp(sym->name, "ieee_selected_real_kind")
+	    && sym->from_intmod == INTMOD_IEEE_ARITHMETIC)
+	  {
+	    gfc_expr *new_expr = gfc_simplify_ieee_selected_real_kind (e);
+	    if (new_expr)
+	      {
+		gfc_replace_expr (e, new_expr);
+		t = true;
+		break;
+	      }
+	  }
 
-	sym = e->symtree->n.sym;
 	if (!gfc_is_intrinsic (sym, 0, e->where)
 	    || (m = gfc_intrinsic_func_interface (e, 0)) != MATCH_YES)
 	  {
diff --git a/gcc/fortran/gfortran.h b/gcc/fortran/gfortran.h
index 1c4638f..3481319 100644
--- a/gcc/fortran/gfortran.h
+++ b/gcc/fortran/gfortran.h
@@ -678,7 +678,8 @@ iso_c_binding_symbol;
 
 typedef enum
 {
-  INTMOD_NONE = 0, INTMOD_ISO_FORTRAN_ENV, INTMOD_ISO_C_BINDING
+  INTMOD_NONE = 0, INTMOD_ISO_FORTRAN_ENV, INTMOD_ISO_C_BINDING,
+  INTMOD_IEEE_FEATURES, INTMOD_IEEE_EXCEPTIONS, INTMOD_IEEE_ARITHMETIC
 }
 intmod_id;
 
@@ -2870,6 +2871,8 @@ gfc_formal_arglist *gfc_sym_get_dummy_args (gfc_symbol *);
 /* intrinsic.c -- true if working in an init-expr, false otherwise.  */
 extern bool gfc_init_expr_flag;
 
+gfc_expr *gfc_simplify_ieee_selected_real_kind (gfc_expr *);
+
 /* Given a symbol that we have decided is intrinsic, mark it as such
    by placing it into a special module that is otherwise impossible to
    read or write.  */
diff --git a/gcc/fortran/intrinsic.texi b/gcc/fortran/intrinsic.texi
index 202063f..87f6478 100644
--- a/gcc/fortran/intrinsic.texi
+++ b/gcc/fortran/intrinsic.texi
@@ -13155,6 +13155,7 @@ Fortran 95 elemental function: @ref{IEOR}
 @menu
 * ISO_FORTRAN_ENV::
 * ISO_C_BINDING::
+* IEEE modules::
 * OpenMP Modules OMP_LIB and OMP_LIB_KINDS::
 @end menu
 
@@ -13366,6 +13367,35 @@ Moreover, the following two named constants are defined:
 
 Both are equivalent to the value @code{NULL} in C.
 
+
+
+@node IEEE modules
+@section IEEE modules: @code{IEEE_EXCEPTIONS}, @code{IEEE_ARITHMETIC}, and @code{IEEE_FEATURES}
+@table @asis
+@item @emph{Standard}:
+Fortran 2003 and later
+@end table
+
+The @code{IEEE_EXCEPTIONS}, @code{IEEE_ARITHMETIC}, and @code{IEEE_FEATURES}
+intrinsic modules provide support for exceptions and IEEE arithmetic, as
+defined in Fortran 2003 and later standards, and the IEC 60559:1989 standard
+(@emph{Binary floating-point arithmetic for microprocessor systems}). These
+modules are only provided on the following supported platforms:
+
+@itemize @bullet
+@item i386 and x86_64 processors
+@item platforms which use the GNU C Library (glibc)
+@item platforms with support for SysV/386 routines for floating point
+interface (including Solaris and BSDs)
+@item platforms with the AIX OS
+@end itemize
+
+For full compliance with the Fortran standards, code using the
+@code{IEEE_EXCEPTIONS} or @code{IEEE_ARITHMETIC} modules should be compiled
+with the following options: @code{-fno-unsafe-math-optimizations
+-frounding-math -fsignaling-nans}.
+
+
 @node OpenMP Modules OMP_LIB and OMP_LIB_KINDS
 @section OpenMP Modules @code{OMP_LIB} and @code{OMP_LIB_KINDS}
 @table @asis
diff --git a/gcc/fortran/libgfortran.h b/gcc/fortran/libgfortran.h
index 230b638..1f8616f 100644
--- a/gcc/fortran/libgfortran.h
+++ b/gcc/fortran/libgfortran.h
@@ -35,13 +35,14 @@ along with GCC; see the file COPYING3.  If not see
 					   obsolescent in later standards.  */
 
 
-/* Bitmasks for the various FPE that can be enabled.  */
-#define GFC_FPE_INVALID    (1<<0)
-#define GFC_FPE_DENORMAL   (1<<1)
-#define GFC_FPE_ZERO       (1<<2)
-#define GFC_FPE_OVERFLOW   (1<<3)
-#define GFC_FPE_UNDERFLOW  (1<<4)
-#define GFC_FPE_INEXACT    (1<<5)
+/* Bitmasks for the various FPE that can be enabled.  These need to be straight integers
+   e.g., 8 instead of (1<<3), because they will be included in Fortran source.  */
+#define GFC_FPE_INVALID      1
+#define GFC_FPE_DENORMAL     2
+#define GFC_FPE_ZERO         4
+#define GFC_FPE_OVERFLOW     8
+#define GFC_FPE_UNDERFLOW   16
+#define GFC_FPE_INEXACT     32
 
 /* Defines for floating-point rounding modes.  */
 #define GFC_FPE_DOWNWARD   1
@@ -49,6 +50,10 @@ along with GCC; see the file COPYING3.  If not see
 #define GFC_FPE_TOWARDZERO 3
 #define GFC_FPE_UPWARD     4
 
+/* Size of the buffer required to store FPU state for any target.
+   In particular, this has to be larger than fenv_t on all glibc targets.
+   Currently, the winner is x86_64 with 32 bytes.  */
+#define GFC_FPE_STATE_BUFFER_SIZE 32
 
 /* Bitmasks for the various runtime checks that can be enabled.  */
 #define GFC_RTCHECK_BOUNDS      (1<<0)
diff --git a/gcc/fortran/module.c b/gcc/fortran/module.c
index ec67960..bd7da1c 100644
--- a/gcc/fortran/module.c
+++ b/gcc/fortran/module.c
@@ -190,6 +190,9 @@ static gzFile module_fp;
 static const char *module_name;
 static gfc_use_list *module_list;
 
+/* If we're reading an intrinsic module, this is its ID.  */
+static intmod_id current_intmod;
+
 /* Content of module.  */
 static char* module_content;
 
@@ -4096,7 +4099,10 @@ mio_symbol (gfc_symbol *sym)
   else
     {
       mio_integer (&intmod);
-      sym->from_intmod = (intmod_id) intmod;
+      if (current_intmod)
+	sym->from_intmod = current_intmod;
+      else
+	sym->from_intmod = (intmod_id) intmod;
     }
   
   mio_integer (&(sym->intmod_sym_id));
@@ -6733,6 +6739,7 @@ gfc_use_module (gfc_use_list *module)
   module_name = module->module_name;
   gfc_rename_list = module->rename;
   only_flag = module->only_flag;
+  current_intmod = INTMOD_NONE;
 
   filename = XALLOCAVEC (char, strlen (module_name) + strlen (MODULE_EXTENSION)
 			       + 1);
@@ -6777,6 +6784,26 @@ gfc_use_module (gfc_use_list *module)
       if (module_fp == NULL && module->intrinsic)
 	gfc_fatal_error ("Can't find an intrinsic module named '%s' at %C",
 			 module_name);
+
+      /* Check for the IEEE modules, so we can mark their symbols
+	 accordingly when we read them.  */
+      if (strcmp (module_name, "ieee_features") == 0
+	  && gfc_notify_std (GFC_STD_F2003, "IEEE_FEATURES module at %C"))
+	{
+	  current_intmod = INTMOD_IEEE_FEATURES;
+	}
+      else if (strcmp (module_name, "ieee_exceptions") == 0
+	       && gfc_notify_std (GFC_STD_F2003,
+				  "IEEE_EXCEPTIONS module at %C"))
+	{
+	  current_intmod = INTMOD_IEEE_EXCEPTIONS;
+	}
+      else if (strcmp (module_name, "ieee_arithmetic") == 0
+	       && gfc_notify_std (GFC_STD_F2003,
+				  "IEEE_ARITHMETIC module at %C"))
+	{
+	  current_intmod = INTMOD_IEEE_ARITHMETIC;
+	}
     }
 
   if (module_fp == NULL)
diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index d18bc08..60d8593 100644
--- a/gcc/fortran/simplify.c
+++ b/gcc/fortran/simplify.c
@@ -5460,12 +5460,13 @@ gfc_simplify_selected_real_kind (gfc_expr *p, gfc_expr *q, gfc_expr *rdx)
       if (gfc_real_kinds[i].range >= range)
 	found_range = 1;
 
-      if (gfc_real_kinds[i].radix >= radix)
+      if (radix == 0 || gfc_real_kinds[i].radix == radix)
 	found_radix = 1;
 
       if (gfc_real_kinds[i].precision >= precision
 	  && gfc_real_kinds[i].range >= range
-	  && gfc_real_kinds[i].radix >= radix && gfc_real_kinds[i].kind < kind)
+	  && (radix == 0 || gfc_real_kinds[i].radix == radix)
+	  && gfc_real_kinds[i].kind < kind)
 	kind = gfc_real_kinds[i].kind;
     }
 
@@ -5488,6 +5489,87 @@ gfc_simplify_selected_real_kind (gfc_expr *p, gfc_expr *q, gfc_expr *rdx)
 
 
 gfc_expr *
+gfc_simplify_ieee_selected_real_kind (gfc_expr *expr)
+{
+  gfc_actual_arglist *arg = expr->value.function.actual;
+  gfc_expr *p = arg->expr, *r = arg->next->expr,
+	   *rad = arg->next->next->expr;
+  int precision, range, radix, res;
+  int found_precision, found_range, found_radix, i;
+
+  if (p)
+  {
+    if (p->expr_type != EXPR_CONSTANT
+	|| gfc_extract_int (p, &precision) != NULL)
+      return NULL;
+  }
+  else
+    precision = 0;
+
+  if (r)
+  {
+    if (r->expr_type != EXPR_CONSTANT
+	|| gfc_extract_int (r, &range) != NULL)
+      return NULL;
+  }
+  else
+    range = 0;
+
+  if (rad)
+  {
+    if (rad->expr_type != EXPR_CONSTANT
+	|| gfc_extract_int (rad, &radix) != NULL)
+      return NULL;
+  }
+  else
+    radix = 0;
+
+  res = INT_MAX;
+  found_precision = 0;
+  found_range = 0;
+  found_radix = 0;
+
+  for (i = 0; gfc_real_kinds[i].kind != 0; i++)
+    {
+      /* We only support the target's float and double types.  */
+      if (!gfc_real_kinds[i].c_float && !gfc_real_kinds[i].c_double)
+	continue;
+
+      if (gfc_real_kinds[i].precision >= precision)
+	found_precision = 1;
+
+      if (gfc_real_kinds[i].range >= range)
+	found_range = 1;
+
+      if (radix == 0 || gfc_real_kinds[i].radix == radix)
+	found_radix = 1;
+
+      if (gfc_real_kinds[i].precision >= precision
+	  && gfc_real_kinds[i].range >= range
+	  && (radix == 0 || gfc_real_kinds[i].radix == radix)
+	  && gfc_real_kinds[i].kind < res)
+	res = gfc_real_kinds[i].kind;
+    }
+
+  if (res == INT_MAX)
+    {
+      if (found_radix && found_range && !found_precision)
+	res = -1;
+      else if (found_radix && found_precision && !found_range)
+	res = -2;
+      else if (found_radix && !found_precision && !found_range)
+	res = -3;
+      else if (found_radix)
+	res = -4;
+      else
+	res = -5;
+    }
+
+  return gfc_get_int_expr (gfc_default_integer_kind, &expr->where, res);
+}
+
+
+gfc_expr *
 gfc_simplify_set_exponent (gfc_expr *x, gfc_expr *i)
 {
   gfc_expr *result;
diff --git a/gcc/fortran/trans-decl.c b/gcc/fortran/trans-decl.c
index 291dd1f..cbcd52d 100644
--- a/gcc/fortran/trans-decl.c
+++ b/gcc/fortran/trans-decl.c
@@ -90,6 +90,9 @@ static stmtblock_t caf_init_block;
 tree gfc_static_ctors;
 
 
+/* Whether we've seen a symbol from an IEEE module in the namespace.  */
+static int seen_ieee_symbol;
+
 /* Function declarations for builtin library functions.  */
 
 tree gfor_fndecl_pause_numeric;
@@ -118,6 +121,8 @@ tree gfor_fndecl_in_unpack;
 tree gfor_fndecl_associated;
 tree gfor_fndecl_system_clock4;
 tree gfor_fndecl_system_clock8;
+tree gfor_fndecl_ieee_procedure_entry;
+tree gfor_fndecl_ieee_procedure_exit;
 
 
 /* Coarray run-time library function decls.  */
@@ -1376,8 +1381,9 @@ gfc_get_symbol_decl (gfc_symbol * sym)
 
   /* Special case for array-valued named constants from intrinsic
      procedures; those are inlined.  */
-  if (sym->attr.use_assoc && sym->from_intmod
-      && sym->attr.flavor == FL_PARAMETER)
+  if (sym->attr.use_assoc && sym->attr.flavor == FL_PARAMETER
+      && (sym->from_intmod == INTMOD_ISO_FORTRAN_ENV
+	  || sym->from_intmod == INTMOD_ISO_C_BINDING))
     intrinsic_array_parameter = true;
 
   /* If use associated compilation, use the module
@@ -3269,6 +3275,14 @@ gfc_build_builtin_function_decls (void)
 	get_identifier (PREFIX("set_fpe")),
 	void_type_node, 1, integer_type_node);
 
+  gfor_fndecl_ieee_procedure_entry = gfc_build_library_function_decl (
+	get_identifier (PREFIX("ieee_procedure_entry")),
+	void_type_node, 1, pvoid_type_node);
+
+  gfor_fndecl_ieee_procedure_exit = gfc_build_library_function_decl (
+	get_identifier (PREFIX("ieee_procedure_exit")),
+	void_type_node, 1, pvoid_type_node);
+
   /* Keep the array dimension in sync with the call, later in this file.  */
   gfor_fndecl_set_options = gfc_build_library_function_decl_with_spec (
 	get_identifier (PREFIX("set_options")), "..R",
@@ -5530,6 +5544,55 @@ gfc_generate_return (void)
 }
 
 
+static void
+is_from_ieee_module (gfc_symbol *sym)
+{
+  if (sym->from_intmod == INTMOD_IEEE_FEATURES
+      || sym->from_intmod == INTMOD_IEEE_EXCEPTIONS
+      || sym->from_intmod == INTMOD_IEEE_ARITHMETIC)
+    seen_ieee_symbol = 1;
+}
+
+
+static int
+is_ieee_module_used (gfc_namespace *ns)
+{
+  seen_ieee_symbol = 0;
+  gfc_traverse_ns (ns, is_from_ieee_module);
+  return seen_ieee_symbol;
+}
+
+
+static tree
+save_fp_state (stmtblock_t *block)
+{
+  tree type, fpstate, tmp;
+
+  type = build_array_type (char_type_node,
+	                   build_range_type (size_type_node, size_zero_node,
+					     size_int (32)));
+  fpstate = gfc_create_var (type, "fpstate");
+  fpstate = gfc_build_addr_expr (pvoid_type_node, fpstate);
+
+  tmp = build_call_expr_loc (input_location, gfor_fndecl_ieee_procedure_entry,
+			     1, fpstate);
+  gfc_add_expr_to_block (block, tmp);
+
+  return fpstate;
+}
+
+
+static void
+restore_fp_state (stmtblock_t *block, tree fpstate)
+{
+  tree tmp;
+
+  tmp = build_call_expr_loc (input_location, gfor_fndecl_ieee_procedure_exit,
+			     1, fpstate);
+  gfc_add_expr_to_block (block, tmp);
+}
+
+
 /* Generate code for a function.  */
 
 void
@@ -5539,13 +5602,14 @@ gfc_generate_function_code (gfc_namespace * ns)
   tree old_context;
   tree decl;
   tree tmp;
+  tree fpstate = NULL_TREE;
   stmtblock_t init, cleanup;
   stmtblock_t body;
   gfc_wrapped_block try_block;
   tree recurcheckvar = NULL_TREE;
   gfc_symbol *sym;
   gfc_symbol *previous_procedure_symbol;
-  int rank;
+  int rank, ieee;
   bool is_recursive;
 
   sym = ns->proc_name;
@@ -5636,6 +5700,12 @@ gfc_generate_function_code (gfc_namespace * ns)
       free (msg);
     }
 
+  /* Check if an IEEE module is used in the procedure.  If so, save
+     the floating point state.  */
+  ieee = is_ieee_module_used (ns);
+  if (ieee)
+    fpstate = save_fp_state (&init);
+
   /* Now generate the code for the body of this function.  */
   gfc_init_block (&body);
 
@@ -5719,6 +5789,10 @@ gfc_generate_function_code (gfc_namespace * ns)
       recurcheckvar = NULL;
     }
 
+  /* If IEEE modules are loaded, restore the floating-point state.  */
+  if (ieee)
+    restore_fp_state (&cleanup, fpstate);
+
   /* Finish the function body and add init and cleanup code.  */
   tmp = gfc_finish_block (&body);
   gfc_start_wrapped_block (&try_block, tmp);
diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog
index df79f3b..739e0aa 100644
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,3 +1,18 @@
+2014-06-28  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
+
+	PR fortran/29383
+	* lib/target-supports.exp (check_effective_target_fortran_ieee): 
+	New function.
+	* gfortran.dg/ieee/ieee.exp: New file.
+	* gfortran.dg/ieee/ieee_1.F90: New file.
+	* gfortran.dg/ieee/ieee_2.f90: New file.
+	* gfortran.dg/ieee/ieee_3.f90: New file.
+	* gfortran.dg/ieee/ieee_4.f90: New file.
+	* gfortran.dg/ieee/ieee_5.f90: New file.
+	* gfortran.dg/ieee/ieee_6.f90: New file.
+	* gfortran.dg/ieee/ieee_7.f90: New file.
+	* gfortran.dg/ieee/ieee_rounding_1.f90: New file.
+
 2014-06-28  Jonathan Wakely  <jwakely@redhat.com>
 
 	* g++.dg/cpp0x/elision_conv.C: New.
diff --git a/gcc/testsuite/gfortran.dg/ieee/ieee.exp b/gcc/testsuite/gfortran.dg/ieee/ieee.exp
new file mode 100644
index 0000000..77e63b7
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/ieee.exp
@@ -0,0 +1,59 @@
+# Copyright (C) 2013 Free Software Foundation, Inc.
+#
+# This file is part of GCC.
+#
+# GCC is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3, or (at your option)
+# any later version.
+#
+# GCC is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with GCC; see the file COPYING3.  If not see
+# <http://www.gnu.org/licenses/>.
+
+# GCC testsuite that uses the `dg.exp' driver.
+
+# Load support procs.
+load_lib gfortran-dg.exp
+load_lib target-supports.exp
+
+# Initialize `dg'.
+dg-init
+
+# Flags specified in each test
+global DEFAULT_FFLAGS
+if ![info exists DEFAULT_FFLAGS] then {
+    set DEFAULT_FFLAGS ""
+}
+
+# Flags for finding the IEEE modules
+if [info exists TOOL_OPTIONS] {
+   set specpath [get_multilibs ${TOOL_OPTIONS}]
+} else {
+   set specpath [get_multilibs]
+}
+set options "-fintrinsic-modules-path $specpath/libgfortran/"
+
+# Bail out if IEEE tests are not supported at all
+if ![check_effective_target_fortran_ieee $options ] {
+ return
+}
+
+# Add target-independent options to require IEEE compatibility
+set options "$DEFAULT_FFLAGS $options -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans"
+
+# Add target-specific options to require IEEE compatibility
+set target_options [add_options_for_ieee ""]
+set options "$options $target_options"
+
+# Main loop.
+gfortran-dg-runtest [lsort \
+      [find $srcdir/$subdir *.\[fF\]{,90,95,03,08} ] ] $options
+
+# All done.
+dg-finish
diff --git a/gcc/testsuite/gfortran.dg/ieee/ieee_1.F90 b/gcc/testsuite/gfortran.dg/ieee/ieee_1.F90
new file mode 100644
index 0000000..9c1c4e3
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/ieee_1.F90
@@ -0,0 +1,174 @@
+! { dg-do run }
+! { dg-additional-options "-ffree-line-length-none -O0" }
+!
+! Use dg-additional-options rather than dg-options to avoid overwriting the
+! default IEEE options which are passed by ieee.exp and necessary.
+
+  use ieee_features, only : ieee_datatype, ieee_denormal, ieee_divide, &
+      ieee_halting, ieee_inexact_flag, ieee_inf, ieee_invalid_flag, &
+      ieee_nan, ieee_rounding, ieee_sqrt, ieee_underflow_flag
+  use ieee_exceptions
+
+  implicit none
+
+  interface use_real
+    procedure use_real_4, use_real_8
+  end interface use_real
+
+  type(ieee_flag_type), parameter :: x(5) = &
+    [ IEEE_INVALID, IEEE_OVERFLOW, IEEE_DIVIDE_BY_ZERO, &
+      IEEE_UNDERFLOW, IEEE_INEXACT ]
+  logical :: l(5) = .false.
+  character(len=5) :: s
+
+#define FLAGS_STRING(S) \
+  call ieee_get_flag(x, l) ; \
+  write(S,"(5(A1))") merge(["I","O","Z","U","P"],[" "," "," "," "," "],l)
+
+#define CHECK_FLAGS(expected) \
+  FLAGS_STRING(s) ; \
+  if (s /= expected) then ; \
+    write (*,"(A,I0,A,A)") "Flags at line ", __LINE__, ": ", s ; \
+    call abort ; \
+  end if ; \
+  call check_flag_sub
+
+  real :: sx
+  double precision :: dx
+
+  ! This file tests IEEE_SET_FLAG and IEEE_GET_FLAG
+
+  !!!! IEEE float
+
+  ! Initial flags are all off
+  CHECK_FLAGS("     ")
+
+  ! Check we can clear them
+  call ieee_set_flag(ieee_all, .false.)
+  CHECK_FLAGS("     ")
+
+  ! Raise invalid, then clear
+  sx = -1
+  call use_real(sx)
+  sx = sqrt(sx)
+  call use_real(sx)
+  CHECK_FLAGS("I    ")
+  call ieee_set_flag(ieee_all, .false.)
+  CHECK_FLAGS("     ")
+
+  ! Raise overflow and precision
+  sx = huge(sx)
+  CHECK_FLAGS("     ")
+  sx = sx*sx
+  CHECK_FLAGS(" O  P")
+  call use_real(sx)
+
+  ! Also raise divide-by-zero
+  sx = 0
+  sx = 1 / sx
+  CHECK_FLAGS(" OZ P")
+  call use_real(sx)
+
+  ! Clear them
+  call ieee_set_flag([ieee_overflow,ieee_inexact,&
+                      ieee_divide_by_zero],[.false.,.false.,.true.])
+  CHECK_FLAGS("  Z  ")
+  call ieee_set_flag(ieee_divide_by_zero, .false.)
+  CHECK_FLAGS("     ")
+
+  ! Raise underflow
+  sx = tiny(sx)
+  CHECK_FLAGS("     ")
+  sx = sx / 10
+  call use_real(sx)
+  CHECK_FLAGS("   UP")
+
+  ! Raise everything
+  call ieee_set_flag(ieee_all, .true.)
+  CHECK_FLAGS("IOZUP")
+
+  ! And clear
+  call ieee_set_flag(ieee_all, .false.)
+  CHECK_FLAGS("     ")
+
+  !!!! IEEE double
+
+  ! Initial flags are all off
+  CHECK_FLAGS("     ")
+
+  ! Check we can clear them
+  call ieee_set_flag(ieee_all, .false.)
+  CHECK_FLAGS("     ")
+
+  ! Raise invalid, then clear
+  dx = -1
+  call use_real(dx)
+  dx = sqrt(dx)
+  call use_real(dx)
+  CHECK_FLAGS("I    ")
+  call ieee_set_flag(ieee_all, .false.)
+  CHECK_FLAGS("     ")
+
+  ! Raise overflow and precision
+  dx = huge(dx)
+  CHECK_FLAGS("     ")
+  dx = dx*dx
+  CHECK_FLAGS(" O  P")
+  call use_real(dx)
+
+  ! Also raise divide-by-zero
+  dx = 0
+  dx = 1 / dx
+  CHECK_FLAGS(" OZ P")
+  call use_real(dx)
+
+  ! Clear them
+  call ieee_set_flag([ieee_overflow,ieee_inexact,&
+                      ieee_divide_by_zero],[.false.,.false.,.true.])
+  CHECK_FLAGS("  Z  ")
+  call ieee_set_flag(ieee_divide_by_zero, .false.)
+  CHECK_FLAGS("     ")
+
+  ! Raise underflow
+  dx = tiny(dx)
+  CHECK_FLAGS("     ")
+  dx = dx / 10
+  CHECK_FLAGS("   UP")
+  call use_real(dx)
+
+  ! Raise everything
+  call ieee_set_flag(ieee_all, .true.)
+  CHECK_FLAGS("IOZUP")
+
+  ! And clear
+  call ieee_set_flag(ieee_all, .false.)
+  CHECK_FLAGS("     ")
+
+contains
+
+  subroutine check_flag_sub
+    use ieee_exceptions
+    logical :: l(5) = .false.
+    type(ieee_flag_type), parameter :: x(5) = &
+      [ IEEE_INVALID, IEEE_OVERFLOW, IEEE_DIVIDE_BY_ZERO, &
+        IEEE_UNDERFLOW, IEEE_INEXACT ]
+    call ieee_get_flag(x, l)
+
+    if (any(l)) then
+      print *, "Flags not cleared in subroutine"
+      call abort
+    end if
+  end subroutine
+
+  ! Interface to a routine that avoids calculations to be optimized out,
+  ! making it appear that we use the result
+  subroutine use_real_4(x)
+    real :: x
+    if (x == 123456.789) print *, "toto"
+  end subroutine
+  subroutine use_real_8(x)
+    double precision :: x
+    if (x == 123456.789) print *, "toto"
+  end subroutine
+
+end
diff --git a/gcc/testsuite/gfortran.dg/ieee/ieee_2.f90 b/gcc/testsuite/gfortran.dg/ieee/ieee_2.f90
new file mode 100644
index 0000000..b138061
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/ieee_2.f90
@@ -0,0 +1,413 @@
+! { dg-do run }
+
+  use, intrinsic :: ieee_features
+  use, intrinsic :: ieee_exceptions
+  use, intrinsic :: ieee_arithmetic
+  implicit none
+
+  interface check_equal
+    procedure check_equal_float, check_equal_double
+  end interface
+
+  interface check_not_equal
+    procedure check_not_equal_float, check_not_equal_double
+  end interface
+
+  real :: sx1, sx2, sx3
+  double precision :: dx1, dx2, dx3
+  type(ieee_round_type) :: mode
+
+  ! Test IEEE_COPY_SIGN
+  sx1 = 1.3
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -sx1) call abort
+  sx1 = huge(sx1)
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -sx1) call abort
+  sx1 = ieee_value(sx1, ieee_positive_inf)
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -sx1) call abort
+  sx1 = tiny(sx1)
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -sx1) call abort
+  sx1 = tiny(sx1)
+  sx1 = sx1 / 101
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -sx1) call abort
+
+  sx1 = -1.3
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= abs(sx1)) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -abs(sx1)) call abort
+  sx1 = -huge(sx1)
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= abs(sx1)) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -abs(sx1)) call abort
+  sx1 = ieee_value(sx1, ieee_negative_inf)
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= abs(sx1)) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -abs(sx1)) call abort
+  sx1 = -tiny(sx1)
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= abs(sx1)) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -abs(sx1)) call abort
+  sx1 = -tiny(sx1)
+  sx1 = sx1 / 101
+  if (ieee_copy_sign(sx1, sx1) /= sx1) call abort
+  if (ieee_copy_sign(sx1, -sx1) /= -sx1) call abort
+  if (ieee_copy_sign(sx1, 1.) /= abs(sx1)) call abort
+  if (ieee_copy_sign(sx1, -1.) /= -abs(sx1)) call abort
+
+  if (ieee_class(ieee_copy_sign(0., -1.)) /= ieee_negative_zero) call abort
+  if (ieee_class(ieee_copy_sign(-0., -1.)) /= ieee_negative_zero) call abort
+  if (ieee_class(ieee_copy_sign(0., 1.)) /= ieee_positive_zero) call abort
+  if (ieee_class(ieee_copy_sign(-0., 1.)) /= ieee_positive_zero) call abort
+
+  sx1 = ieee_value(0., ieee_quiet_nan)
+  if (ieee_class(ieee_copy_sign(sx1, 1.)) /= ieee_quiet_nan) call abort
+  if (ieee_class(ieee_copy_sign(sx1, -1.)) /= ieee_quiet_nan) call abort
+
+  dx1 = 1.3
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -1.d0) /= -dx1) call abort
+  dx1 = huge(dx1)
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.d0) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -1.) /= -dx1) call abort
+  dx1 = ieee_value(dx1, ieee_positive_inf)
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -1.d0) /= -dx1) call abort
+  dx1 = tiny(dx1)
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.d0) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -1.) /= -dx1) call abort
+  dx1 = tiny(dx1)
+  dx1 = dx1 / 101
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -1.d0) /= -dx1) call abort
+
+  dx1 = -1.3d0
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.d0) /= abs(dx1)) call abort
+  if (ieee_copy_sign(dx1, -1.) /= -abs(dx1)) call abort
+  dx1 = -huge(dx1)
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.) /= abs(dx1)) call abort
+  if (ieee_copy_sign(dx1, -1.d0) /= -abs(dx1)) call abort
+  dx1 = ieee_value(dx1, ieee_negative_inf)
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.d0) /= abs(dx1)) call abort
+  if (ieee_copy_sign(dx1, -1.) /= -abs(dx1)) call abort
+  dx1 = -tiny(dx1)
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.) /= abs(dx1)) call abort
+  if (ieee_copy_sign(dx1, -1.d0) /= -abs(dx1)) call abort
+  dx1 = -tiny(dx1)
+  dx1 = dx1 / 101
+  if (ieee_copy_sign(dx1, dx1) /= dx1) call abort
+  if (ieee_copy_sign(dx1, -dx1) /= -dx1) call abort
+  if (ieee_copy_sign(dx1, 1.d0) /= abs(dx1)) call abort
+  if (ieee_copy_sign(dx1, -1.) /= -abs(dx1)) call abort
+
+  if (ieee_class(ieee_copy_sign(0.d0, -1.)) /= ieee_negative_zero) call abort
+  if (ieee_class(ieee_copy_sign(-0.d0, -1.)) /= ieee_negative_zero) call abort
+  if (ieee_class(ieee_copy_sign(0.d0, 1.)) /= ieee_positive_zero) call abort
+  if (ieee_class(ieee_copy_sign(-0.d0, 1.)) /= ieee_positive_zero) call abort
+
+  dx1 = ieee_value(0.d0, ieee_quiet_nan)
+  if (ieee_class(ieee_copy_sign(dx1, 1.d0)) /= ieee_quiet_nan) call abort
+  if (ieee_class(ieee_copy_sign(dx1, -1.)) /= ieee_quiet_nan) call abort
+
+  ! Test IEEE_LOGB
+
+  if (ieee_logb(1.17) /= exponent(1.17) - 1) call abort
+  if (ieee_logb(-1.17) /= exponent(-1.17) - 1) call abort
+  if (ieee_logb(huge(sx1)) /= exponent(huge(sx1)) - 1) call abort
+  if (ieee_logb(-huge(sx1)) /= exponent(-huge(sx1)) - 1) call abort
+  if (ieee_logb(tiny(sx1)) /= exponent(tiny(sx1)) - 1) call abort
+  if (ieee_logb(-tiny(sx1)) /= exponent(-tiny(sx1)) - 1) call abort
+
+  if (ieee_class(ieee_logb(0.)) /= ieee_negative_inf) call abort
+  if (ieee_class(ieee_logb(-0.)) /= ieee_negative_inf) call abort
+
+  sx1 = ieee_value(sx1, ieee_positive_inf)
+  if (ieee_class(ieee_logb(sx1)) /= ieee_positive_inf) call abort
+  if (ieee_class(ieee_logb(-sx1)) /= ieee_positive_inf) call abort
+
+  sx1 = ieee_value(sx1, ieee_quiet_nan)
+  if (ieee_class(ieee_logb(sx1)) /= ieee_quiet_nan) call abort
+
+  if (ieee_logb(1.17d0) /= exponent(1.17d0) - 1) call abort
+  if (ieee_logb(-1.17d0) /= exponent(-1.17d0) - 1) call abort
+  if (ieee_logb(huge(dx1)) /= exponent(huge(dx1)) - 1) call abort
+  if (ieee_logb(-huge(dx1)) /= exponent(-huge(dx1)) - 1) call abort
+  if (ieee_logb(tiny(dx1)) /= exponent(tiny(dx1)) - 1) call abort
+  if (ieee_logb(-tiny(dx1)) /= exponent(-tiny(dx1)) - 1) call abort
+
+  if (ieee_class(ieee_logb(0.d0)) /= ieee_negative_inf) call abort
+  if (ieee_class(ieee_logb(-0.d0)) /= ieee_negative_inf) call abort
+
+  dx1 = ieee_value(dx1, ieee_positive_inf)
+  if (ieee_class(ieee_logb(dx1)) /= ieee_positive_inf) call abort
+  if (ieee_class(ieee_logb(-dx1)) /= ieee_positive_inf) call abort
+
+  dx1 = ieee_value(dx1, ieee_quiet_nan)
+  if (ieee_class(ieee_logb(dx1)) /= ieee_quiet_nan) call abort
+
+  ! Test IEEE_NEXT_AFTER
+
+  if (ieee_next_after(0.12, 1.0) /= nearest(0.12, 1.0)) call abort
+  if (ieee_next_after(0.12, -1.0) /= nearest(0.12, -1.0)) call abort
+
+  sx1 = 0.12
+  if (ieee_next_after(sx1, sx1) /= sx1) call abort
+  sx1 = -0.12
+  if (ieee_next_after(sx1, sx1) /= sx1) call abort
+  sx1 = huge(sx1)
+  if (ieee_next_after(sx1, sx1) /= sx1) call abort
+  sx1 = tiny(sx1)
+  if (ieee_next_after(sx1, sx1) /= sx1) call abort
+  sx1 = 0
+  if (ieee_next_after(sx1, sx1) /= sx1) call abort
+  sx1 = ieee_value(sx1, ieee_negative_inf)
+  if (ieee_next_after(sx1, sx1) /= sx1) call abort
+  sx1 = ieee_value(sx1, ieee_quiet_nan)
+  if (ieee_class(ieee_next_after(sx1, sx1)) /= ieee_quiet_nan) call abort
+
+  if (ieee_next_after(0., 1.0) <= 0) call abort
+  if (ieee_next_after(0., -1.0) >= 0) call abort
+  sx1 = ieee_next_after(huge(sx1), ieee_value(sx1, ieee_negative_inf))
+  if (.not. sx1 < huge(sx1)) call abort
+  sx1 = ieee_next_after(huge(sx1), ieee_value(sx1, ieee_positive_inf))
+  if (ieee_class(sx1) /= ieee_positive_inf) call abort
+  sx1 = ieee_next_after(-tiny(sx1), 1.0)
+  if (ieee_class(sx1) /= ieee_negative_denormal) call abort
+
+  if (ieee_next_after(0.12d0, 1.0d0) /= nearest(0.12d0, 1.0)) call abort
+  if (ieee_next_after(0.12d0, -1.0) /= nearest(0.12d0, -1.0)) call abort
+
+  dx1 = 0.12
+  if (ieee_next_after(dx1, dx1) /= dx1) call abort
+  dx1 = -0.12
+  if (ieee_next_after(dx1, dx1) /= dx1) call abort
+  dx1 = huge(dx1)
+  if (ieee_next_after(dx1, dx1) /= dx1) call abort
+  dx1 = tiny(dx1)
+  if (ieee_next_after(dx1, dx1) /= dx1) call abort
+  dx1 = 0
+  if (ieee_next_after(dx1, dx1) /= dx1) call abort
+  dx1 = ieee_value(dx1, ieee_negative_inf)
+  if (ieee_next_after(dx1, dx1) /= dx1) call abort
+  dx1 = ieee_value(dx1, ieee_quiet_nan)
+  if (ieee_class(ieee_next_after(dx1, dx1)) /= ieee_quiet_nan) call abort
+
+  if (ieee_next_after(0.d0, 1.0) <= 0) call abort
+  if (ieee_next_after(0.d0, -1.0d0) >= 0) call abort
+  dx1 = ieee_next_after(huge(dx1), ieee_value(dx1, ieee_negative_inf))
+  if (.not. dx1 < huge(dx1)) call abort
+  dx1 = ieee_next_after(huge(dx1), ieee_value(dx1, ieee_positive_inf))
+  if (ieee_class(dx1) /= ieee_positive_inf) call abort
+  dx1 = ieee_next_after(-tiny(dx1), 1.0d0)
+  if (ieee_class(dx1) /= ieee_negative_denormal) call abort
+
+  ! Test IEEE_REM
+
+  if (ieee_rem(4.0, 3.0) /= 1.0) call abort
+  if (ieee_rem(-4.0, 3.0) /= -1.0) call abort
+  if (ieee_rem(2.0, 3.0d0) /= -1.0d0) call abort
+  if (ieee_rem(-2.0, 3.0d0) /= 1.0d0) call abort
+  if (ieee_rem(2.0d0, 3.0d0) /= -1.0d0) call abort
+  if (ieee_rem(-2.0d0, 3.0d0) /= 1.0d0) call abort
+
+  if (ieee_class(ieee_rem(ieee_value(0., ieee_quiet_nan), 1.0)) &
+      /= ieee_quiet_nan) call abort
+  if (ieee_class(ieee_rem(1.0, ieee_value(0.d0, ieee_quiet_nan))) &
+      /= ieee_quiet_nan) call abort
+
+  if (ieee_class(ieee_rem(ieee_value(0., ieee_positive_inf), 1.0)) &
+      /= ieee_quiet_nan) call abort
+  if (ieee_class(ieee_rem(ieee_value(0.d0, ieee_negative_inf), 1.0)) &
+      /= ieee_quiet_nan) call abort
+  if (ieee_rem(-1.0, ieee_value(0., ieee_positive_inf)) &
+      /= -1.0) call abort
+  if (ieee_rem(1.0, ieee_value(0.d0, ieee_negative_inf)) &
+      /= 1.0) call abort
+
+
+  ! Test IEEE_RINT
+
+  if (ieee_support_rounding (ieee_nearest, sx1)) then
+    call ieee_get_rounding_mode (mode)
+    call ieee_set_rounding_mode (ieee_nearest)
+    sx1 = 7 / 3.
+    sx1 = ieee_rint (sx1)
+    call ieee_set_rounding_mode (mode)
+    if (sx1 /= 2) call abort
+  end if
+
+  if (ieee_support_rounding (ieee_up, sx1)) then
+    call ieee_get_rounding_mode (mode)
+    call ieee_set_rounding_mode (ieee_up)
+    sx1 = 7 / 3.
+    sx1 = ieee_rint (sx1)
+    call ieee_set_rounding_mode (mode)
+    if (sx1 /= 3) call abort
+  end if
+
+  if (ieee_support_rounding (ieee_down, sx1)) then
+    call ieee_get_rounding_mode (mode)
+    call ieee_set_rounding_mode (ieee_down)
+    sx1 = 7 / 3.
+    sx1 = ieee_rint (sx1)
+    call ieee_set_rounding_mode (mode)
+    if (sx1 /= 2) call abort
+  end if
+
+  if (ieee_support_rounding (ieee_to_zero, sx1)) then
+    call ieee_get_rounding_mode (mode)
+    call ieee_set_rounding_mode (ieee_to_zero)
+    sx1 = 7 / 3.
+    sx1 = ieee_rint (sx1)
+    call ieee_set_rounding_mode (mode)
+    if (sx1 /= 2) call abort
+  end if
+
+  if (ieee_class(ieee_rint(0.)) /= ieee_positive_zero) call abort
+  if (ieee_class(ieee_rint(-0.)) /= ieee_negative_zero) call abort
+
+  if (ieee_support_rounding (ieee_nearest, dx1)) then
+    call ieee_get_rounding_mode (mode)
+    call ieee_set_rounding_mode (ieee_nearest)
+    dx1 = 7 / 3.d0
+    dx1 = ieee_rint (dx1)
+    call ieee_set_rounding_mode (mode)
+    if (dx1 /= 2) call abort
+  end if
+
+  if (ieee_support_rounding (ieee_up, dx1)) then
+    call ieee_get_rounding_mode (mode)
+    call ieee_set_rounding_mode (ieee_up)
+    dx1 = 7 / 3.d0
+    dx1 = ieee_rint (dx1)
+    call ieee_set_rounding_mode (mode)
+    if (dx1 /= 3) call abort
+  end if
+
+  if (ieee_support_rounding (ieee_down, dx1)) then
+    call ieee_get_rounding_mode (mode)
+    call ieee_set_rounding_mode (ieee_down)
+    dx1 = 7 / 3.d0
+    dx1 = ieee_rint (dx1)
+    call ieee_set_rounding_mode (mode)
+    if (dx1 /= 2) call abort
+  end if
+
+  if (ieee_support_rounding (ieee_to_zero, dx1)) then
+    call ieee_get_rounding_mode (mode)
+    call ieee_set_rounding_mode (ieee_to_zero)
+    dx1 = 7 / 3.d0
+    dx1 = ieee_rint (dx1)
+    call ieee_set_rounding_mode (mode)
+    if (dx1 /= 2) call abort
+  end if
+
+  if (ieee_class(ieee_rint(0.d0)) /= ieee_positive_zero) call abort
+  if (ieee_class(ieee_rint(-0.d0)) /= ieee_negative_zero) call abort
+
+  ! Test IEEE_SCALB
+
+  sx1 = 1
+  if (ieee_scalb(sx1, 2) /= 4.) call abort
+  if (ieee_scalb(-sx1, 2) /= -4.) call abort
+  if (ieee_scalb(sx1, -2) /= 1/4.) call abort
+  if (ieee_scalb(-sx1, -2) /= -1/4.) call abort
+  if (ieee_class(ieee_scalb(sx1, huge(0))) /= ieee_positive_inf) call abort
+  if (ieee_class(ieee_scalb(-sx1, huge(0))) /= ieee_negative_inf) call abort
+  if (ieee_class(ieee_scalb(sx1, -huge(0))) /= ieee_positive_zero) call abort
+  if (ieee_class(ieee_scalb(-sx1, -huge(0))) /= ieee_negative_zero) call abort
+
+  sx1 = ieee_value(sx1, ieee_quiet_nan)
+  if (ieee_class(ieee_scalb(sx1, 1)) /= ieee_quiet_nan) call abort
+  sx1 = ieee_value(sx1, ieee_positive_inf)
+  if (ieee_class(ieee_scalb(sx1, -42)) /= ieee_positive_inf) call abort
+  sx1 = ieee_value(sx1, ieee_negative_inf)
+  if (ieee_class(ieee_scalb(sx1, -42)) /= ieee_negative_inf) call abort
+
+  dx1 = 1
+  if (ieee_scalb(dx1, 2) /= 4.d0) call abort
+  if (ieee_scalb(-dx1, 2) /= -4.d0) call abort
+  if (ieee_scalb(dx1, -2) /= 1/4.d0) call abort
+  if (ieee_scalb(-dx1, -2) /= -1/4.d0) call abort
+  if (ieee_class(ieee_scalb(dx1, huge(0))) /= ieee_positive_inf) call abort
+  if (ieee_class(ieee_scalb(-dx1, huge(0))) /= ieee_negative_inf) call abort
+  if (ieee_class(ieee_scalb(dx1, -huge(0))) /= ieee_positive_zero) call abort
+  if (ieee_class(ieee_scalb(-dx1, -huge(0))) /= ieee_negative_zero) call abort
+
+  dx1 = ieee_value(dx1, ieee_quiet_nan)
+  if (ieee_class(ieee_scalb(dx1, 1)) /= ieee_quiet_nan) call abort
+  dx1 = ieee_value(dx1, ieee_positive_inf)
+  if (ieee_class(ieee_scalb(dx1, -42)) /= ieee_positive_inf) call abort
+  dx1 = ieee_value(dx1, ieee_negative_inf)
+  if (ieee_class(ieee_scalb(dx1, -42)) /= ieee_negative_inf) call abort
+
+contains
+
+  subroutine check_equal_float (x, y)
+    real, intent(in) :: x, y
+    if (x /= y) then
+      print *, x, y
+      call abort
+    end if
+  end subroutine
+
+  subroutine check_equal_double (x, y)
+    double precision, intent(in) :: x, y
+    if (x /= y) then
+      print *, x, y
+      call abort
+    end if
+  end subroutine
+
+  subroutine check_not_equal_float (x, y)
+    real, intent(in) :: x, y
+    if (x == y) then
+      print *, x, y
+      call abort
+    end if
+  end subroutine
+
+  subroutine check_not_equal_double (x, y)
+    double precision, intent(in) :: x, y
+    if (x == y) then
+      print *, x, y
+      call abort
+    end if
+  end subroutine
+
+end
diff --git a/gcc/testsuite/gfortran.dg/ieee/ieee_3.f90 b/gcc/testsuite/gfortran.dg/ieee/ieee_3.f90
new file mode 100644
index 0000000..b2c7186
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/ieee_3.f90
@@ -0,0 +1,167 @@
+! { dg-do run }
+
+  use :: ieee_arithmetic
+  implicit none
+
+  real :: sx1, sx2, sx3
+  double precision :: dx1, dx2, dx3
+  integer, parameter :: s = kind(sx1), d = kind(dx1)
+  type(ieee_round_type) :: mode
+
+  ! Test IEEE_IS_FINITE
+
+  if (ieee_support_datatype(0._s)) then
+    if (.not. ieee_is_finite(0.2_s)) call abort
+    if (.not. ieee_is_finite(-0.2_s)) call abort
+    if (.not. ieee_is_finite(0._s)) call abort
+    if (.not. ieee_is_finite(-0._s)) call abort
+    if (.not. ieee_is_finite(tiny(0._s))) call abort
+    if (.not. ieee_is_finite(tiny(0._s)/100)) call abort
+    if (.not. ieee_is_finite(huge(0._s))) call abort
+    if (.not. ieee_is_finite(-huge(0._s))) call abort
+    sx1 = huge(sx1)
+    if (ieee_is_finite(2*sx1)) call abort
+    if (ieee_is_finite(2*(-sx1))) call abort
+    sx1 = ieee_value(sx1, ieee_quiet_nan)
+    if (ieee_is_finite(sx1)) call abort
+  end if
+
+  if (ieee_support_datatype(0._d)) then
+    if (.not. ieee_is_finite(0.2_d)) call abort
+    if (.not. ieee_is_finite(-0.2_d)) call abort
+    if (.not. ieee_is_finite(0._d)) call abort
+    if (.not. ieee_is_finite(-0._d)) call abort
+    if (.not. ieee_is_finite(tiny(0._d))) call abort
+    if (.not. ieee_is_finite(tiny(0._d)/100)) call abort
+    if (.not. ieee_is_finite(huge(0._d))) call abort
+    if (.not. ieee_is_finite(-huge(0._d))) call abort
+    dx1 = huge(dx1)
+    if (ieee_is_finite(2*dx1)) call abort
+    if (ieee_is_finite(2*(-dx1))) call abort
+    dx1 = ieee_value(dx1, ieee_quiet_nan)
+    if (ieee_is_finite(dx1)) call abort
+  end if
+
+  ! Test IEEE_IS_NAN
+
+  if (ieee_support_datatype(0._s)) then
+    if (ieee_is_nan(0.2_s)) call abort
+    if (ieee_is_nan(-0.2_s)) call abort
+    if (ieee_is_nan(0._s)) call abort
+    if (ieee_is_nan(-0._s)) call abort
+    if (ieee_is_nan(tiny(0._s))) call abort
+    if (ieee_is_nan(tiny(0._s)/100)) call abort
+    if (ieee_is_nan(huge(0._s))) call abort
+    if (ieee_is_nan(-huge(0._s))) call abort
+    sx1 = huge(sx1)
+    if (ieee_is_nan(2*sx1)) call abort
+    if (ieee_is_nan(2*(-sx1))) call abort
+    sx1 = ieee_value(sx1, ieee_quiet_nan)
+    if (.not. ieee_is_nan(sx1)) call abort
+    sx1 = -1
+    if (.not. ieee_is_nan(sqrt(sx1))) call abort
+  end if
+
+  if (ieee_support_datatype(0._d)) then
+    if (ieee_is_nan(0.2_d)) call abort
+    if (ieee_is_nan(-0.2_d)) call abort
+    if (ieee_is_nan(0._d)) call abort
+    if (ieee_is_nan(-0._d)) call abort
+    if (ieee_is_nan(tiny(0._d))) call abort
+    if (ieee_is_nan(tiny(0._d)/100)) call abort
+    if (ieee_is_nan(huge(0._d))) call abort
+    if (ieee_is_nan(-huge(0._d))) call abort
+    dx1 = huge(dx1)
+    if (ieee_is_nan(2*dx1)) call abort
+    if (ieee_is_nan(2*(-dx1))) call abort
+    dx1 = ieee_value(dx1, ieee_quiet_nan)
+    if (.not. ieee_is_nan(dx1)) call abort
+    dx1 = -1
+    if (.not. ieee_is_nan(sqrt(dx1))) call abort
+  end if
+
+  ! IEEE_IS_NEGATIVE
+
+  if (ieee_support_datatype(0._s)) then
+    if (ieee_is_negative(0.2_s)) call abort
+    if (.not. ieee_is_negative(-0.2_s)) call abort
+    if (ieee_is_negative(0._s)) call abort
+    if (.not. ieee_is_negative(-0._s)) call abort
+    if (ieee_is_negative(tiny(0._s))) call abort
+    if (ieee_is_negative(tiny(0._s)/100)) call abort
+    if (.not. ieee_is_negative(-tiny(0._s))) call abort
+    if (.not. ieee_is_negative(-tiny(0._s)/100)) call abort
+    if (ieee_is_negative(huge(0._s))) call abort
+    if (.not. ieee_is_negative(-huge(0._s))) call abort
+    sx1 = huge(sx1)
+    if (ieee_is_negative(2*sx1)) call abort
+    if (.not. ieee_is_negative(2*(-sx1))) call abort
+    sx1 = ieee_value(sx1, ieee_quiet_nan)
+    if (ieee_is_negative(sx1)) call abort
+    sx1 = -1
+    if (ieee_is_negative(sqrt(sx1))) call abort
+  end if
+
+  if (ieee_support_datatype(0._d)) then
+    if (ieee_is_negative(0.2_d)) call abort
+    if (.not. ieee_is_negative(-0.2_d)) call abort
+    if (ieee_is_negative(0._d)) call abort
+    if (.not. ieee_is_negative(-0._d)) call abort
+    if (ieee_is_negative(tiny(0._d))) call abort
+    if (ieee_is_negative(tiny(0._d)/100)) call abort
+    if (.not. ieee_is_negative(-tiny(0._d))) call abort
+    if (.not. ieee_is_negative(-tiny(0._d)/100)) call abort
+    if (ieee_is_negative(huge(0._d))) call abort
+    if (.not. ieee_is_negative(-huge(0._d))) call abort
+    dx1 = huge(dx1)
+    if (ieee_is_negative(2*dx1)) call abort
+    if (.not. ieee_is_negative(2*(-dx1))) call abort
+    dx1 = ieee_value(dx1, ieee_quiet_nan)
+    if (ieee_is_negative(dx1)) call abort
+    dx1 = -1
+    if (ieee_is_negative(sqrt(dx1))) call abort
+  end if
+
+  ! Test IEEE_IS_NORMAL
+
+  if (ieee_support_datatype(0._s)) then
+    if (.not. ieee_is_normal(0.2_s)) call abort
+    if (.not. ieee_is_normal(-0.2_s)) call abort
+    if (.not. ieee_is_normal(0._s)) call abort
+    if (.not. ieee_is_normal(-0._s)) call abort
+    if (.not. ieee_is_normal(tiny(0._s))) call abort
+    if (ieee_is_normal(tiny(0._s)/100)) call abort
+    if (.not. ieee_is_normal(-tiny(0._s))) call abort
+    if (ieee_is_normal(-tiny(0._s)/100)) call abort
+    if (.not. ieee_is_normal(huge(0._s))) call abort
+    if (.not. ieee_is_normal(-huge(0._s))) call abort
+    sx1 = huge(sx1)
+    if (ieee_is_normal(2*sx1)) call abort
+    if (ieee_is_normal(2*(-sx1))) call abort
+    sx1 = ieee_value(sx1, ieee_quiet_nan)
+    if (ieee_is_normal(sx1)) call abort
+    sx1 = -1
+    if (ieee_is_normal(sqrt(sx1))) call abort
+  end if
+
+  if (ieee_support_datatype(0._d)) then
+    if (.not. ieee_is_normal(0.2_d)) call abort
+    if (.not. ieee_is_normal(-0.2_d)) call abort
+    if (.not. ieee_is_normal(0._d)) call abort
+    if (.not. ieee_is_normal(-0._d)) call abort
+    if (.not. ieee_is_normal(tiny(0._d))) call abort
+    if (ieee_is_normal(tiny(0._d)/100)) call abort
+    if (.not. ieee_is_normal(-tiny(0._d))) call abort
+    if (ieee_is_normal(-tiny(0._d)/100)) call abort
+    if (.not. ieee_is_normal(huge(0._d))) call abort
+    if (.not. ieee_is_normal(-huge(0._d))) call abort
+    dx1 = huge(dx1)
+    if (ieee_is_normal(2*dx1)) call abort
+    if (ieee_is_normal(2*(-dx1))) call abort
+    dx1 = ieee_value(dx1, ieee_quiet_nan)
+    if (ieee_is_normal(dx1)) call abort
+    dx1 = -1
+    if (ieee_is_normal(sqrt(dx1))) call abort
+  end if
+
+end
diff --git a/gcc/testsuite/gfortran.dg/ieee/ieee_4.f90 b/gcc/testsuite/gfortran.dg/ieee/ieee_4.f90
new file mode 100644
index 0000000..e5f1cee
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/ieee_4.f90
@@ -0,0 +1,189 @@
+! { dg-do run }
+
+  use :: ieee_arithmetic
+  implicit none
+
+  real :: sx1, sx2, sx3
+  double precision :: dx1, dx2, dx3
+  integer, parameter :: s = kind(sx1), d = kind(dx1)
+  type(ieee_round_type) :: mode
+
+  ! Test IEEE_CLASS
+
+  if (ieee_support_datatype(0._s)) then
+    sx1 = 0.1_s
+    if (ieee_class(sx1) /= ieee_positive_normal) call abort
+    if (ieee_class(-sx1) /= ieee_negative_normal) call abort
+    sx1 = huge(sx1)
+    if (ieee_class(sx1) /= ieee_positive_normal) call abort
+    if (ieee_class(-sx1) /= ieee_negative_normal) call abort
+    if (ieee_class(2*sx1) /= ieee_positive_inf) call abort
+    if (ieee_class(2*(-sx1)) /= ieee_negative_inf) call abort
+    sx1 = tiny(sx1)
+    if (ieee_class(sx1) /= ieee_positive_normal) call abort
+    if (ieee_class(-sx1) /= ieee_negative_normal) call abort
+    if (ieee_class(sx1 / 2) /= ieee_positive_denormal) call abort
+    if (ieee_class((-sx1) / 2) /= ieee_negative_denormal) call abort
+    sx1 = -1
+    if (ieee_class(sqrt(sx1)) /= ieee_quiet_nan) call abort
+    sx1 = 0
+    if (ieee_class(sx1) /= ieee_positive_zero) call abort
+    if (ieee_class(-sx1) /= ieee_negative_zero) call abort
+  end if
+
+  if (ieee_support_datatype(0._d)) then
+    dx1 = 0.1_d
+    if (ieee_class(dx1) /= ieee_positive_normal) call abort
+    if (ieee_class(-dx1) /= ieee_negative_normal) call abort
+    dx1 = huge(dx1)
+    if (ieee_class(dx1) /= ieee_positive_normal) call abort
+    if (ieee_class(-dx1) /= ieee_negative_normal) call abort
+    if (ieee_class(2*dx1) /= ieee_positive_inf) call abort
+    if (ieee_class(2*(-dx1)) /= ieee_negative_inf) call abort
+    dx1 = tiny(dx1)
+    if (ieee_class(dx1) /= ieee_positive_normal) call abort
+    if (ieee_class(-dx1) /= ieee_negative_normal) call abort
+    if (ieee_class(dx1 / 2) /= ieee_positive_denormal) call abort
+    if (ieee_class((-dx1) / 2) /= ieee_negative_denormal) call abort
+    dx1 = -1
+    if (ieee_class(sqrt(dx1)) /= ieee_quiet_nan) call abort
+    dx1 = 0
+    if (ieee_class(dx1) /= ieee_positive_zero) call abort
+    if (ieee_class(-dx1) /= ieee_negative_zero) call abort
+  end if
+
+  ! Test IEEE_VALUE and IEEE_UNORDERED
+
+  if (ieee_support_datatype(0._s)) then
+    sx1 = ieee_value(sx1, ieee_quiet_nan)
+    if (.not. ieee_is_nan(sx1)) call abort
+    if (.not. ieee_unordered(sx1, sx1)) call abort
+    if (.not. ieee_unordered(sx1, 0._s)) call abort
+    if (.not. ieee_unordered(sx1, 0._d)) call abort
+    if (.not. ieee_unordered(0._s, sx1)) call abort
+    if (.not. ieee_unordered(0._d, sx1)) call abort
+    if (ieee_unordered(0._s, 0._s)) call abort
+
+    sx1 = ieee_value(sx1, ieee_positive_inf)
+    if (ieee_is_finite(sx1)) call abort
+    if (ieee_is_nan(sx1)) call abort
+    if (ieee_is_negative(sx1)) call abort
+    if (ieee_is_normal(sx1)) call abort
+
+    sx1 = ieee_value(sx1, ieee_negative_inf)
+    if (ieee_is_finite(sx1)) call abort
+    if (ieee_is_nan(sx1)) call abort
+    if (.not. ieee_is_negative(sx1)) call abort
+    if (ieee_is_normal(sx1)) call abort
+
+    sx1 = ieee_value(sx1, ieee_positive_normal)
+    if (.not. ieee_is_finite(sx1)) call abort
+    if (ieee_is_nan(sx1)) call abort
+    if (ieee_is_negative(sx1)) call abort
+    if (.not. ieee_is_normal(sx1)) call abort
+
+    sx1 = ieee_value(sx1, ieee_negative_normal)
+    if (.not. ieee_is_finite(sx1)) call abort
+    if (ieee_is_nan(sx1)) call abort
+    if (.not. ieee_is_negative(sx1)) call abort
+    if (.not. ieee_is_normal(sx1)) call abort
+
+    sx1 = ieee_value(sx1, ieee_positive_denormal)
+    if (.not. ieee_is_finite(sx1)) call abort
+    if (ieee_is_nan(sx1)) call abort
+    if (ieee_is_negative(sx1)) call abort
+    if (ieee_is_normal(sx1)) call abort
+    if (sx1 <= 0) call abort
+    if (sx1 >= tiny(sx1)) call abort
+
+    sx1 = ieee_value(sx1, ieee_negative_denormal)
+    if (.not. ieee_is_finite(sx1)) call abort
+    if (ieee_is_nan(sx1)) call abort
+    if (.not. ieee_is_negative(sx1)) call abort
+    if (ieee_is_normal(sx1)) call abort
+    if (sx1 >= 0) call abort
+    if (sx1 <= -tiny(sx1)) call abort
+
+    sx1 = ieee_value(sx1, ieee_positive_zero)
+    if (.not. ieee_is_finite(sx1)) call abort
+    if (ieee_is_nan(sx1)) call abort
+    if (ieee_is_negative(sx1)) call abort
+    if (.not. ieee_is_normal(sx1)) call abort
+    if (sx1 /= 0) call abort
+
+    sx1 = ieee_value(sx1, ieee_negative_zero)
+    if (.not. ieee_is_finite(sx1)) call abort
+    if (ieee_is_nan(sx1)) call abort
+    if (.not. ieee_is_negative(sx1)) call abort
+    if (.not. ieee_is_normal(sx1)) call abort
+    if (sx1 /= 0) call abort
+
+  end if
+
+  if (ieee_support_datatype(0._d)) then
+    dx1 = ieee_value(dx1, ieee_quiet_nan)
+    if (.not. ieee_is_nan(dx1)) call abort
+    if (.not. ieee_unordered(dx1, dx1)) call abort
+    if (.not. ieee_unordered(dx1, 0._s)) call abort
+    if (.not. ieee_unordered(dx1, 0._d)) call abort
+    if (.not. ieee_unordered(0._s, dx1)) call abort
+    if (.not. ieee_unordered(0._d, dx1)) call abort
+    if (ieee_unordered(0._d, 0._d)) call abort
+
+    dx1 = ieee_value(dx1, ieee_positive_inf)
+    if (ieee_is_finite(dx1)) call abort
+    if (ieee_is_nan(dx1)) call abort
+    if (ieee_is_negative(dx1)) call abort
+    if (ieee_is_normal(dx1)) call abort
+
+    dx1 = ieee_value(dx1, ieee_negative_inf)
+    if (ieee_is_finite(dx1)) call abort
+    if (ieee_is_nan(dx1)) call abort
+    if (.not. ieee_is_negative(dx1)) call abort
+    if (ieee_is_normal(dx1)) call abort
+
+    dx1 = ieee_value(dx1, ieee_positive_normal)
+    if (.not. ieee_is_finite(dx1)) call abort
+    if (ieee_is_nan(dx1)) call abort
+    if (ieee_is_negative(dx1)) call abort
+    if (.not. ieee_is_normal(dx1)) call abort
+
+    dx1 = ieee_value(dx1, ieee_negative_normal)
+    if (.not. ieee_is_finite(dx1)) call abort
+    if (ieee_is_nan(dx1)) call abort
+    if (.not. ieee_is_negative(dx1)) call abort
+    if (.not. ieee_is_normal(dx1)) call abort
+
+    dx1 = ieee_value(dx1, ieee_positive_denormal)
+    if (.not. ieee_is_finite(dx1)) call abort
+    if (ieee_is_nan(dx1)) call abort
+    if (ieee_is_negative(dx1)) call abort
+    if (ieee_is_normal(dx1)) call abort
+    if (dx1 <= 0) call abort
+    if (dx1 >= tiny(dx1)) call abort
+
+    dx1 = ieee_value(dx1, ieee_negative_denormal)
+    if (.not. ieee_is_finite(dx1)) call abort
+    if (ieee_is_nan(dx1)) call abort
+    if (.not. ieee_is_negative(dx1)) call abort
+    if (ieee_is_normal(dx1)) call abort
+    if (dx1 >= 0) call abort
+    if (dx1 <= -tiny(dx1)) call abort
+
+    dx1 = ieee_value(dx1, ieee_positive_zero)
+    if (.not. ieee_is_finite(dx1)) call abort
+    if (ieee_is_nan(dx1)) call abort
+    if (ieee_is_negative(dx1)) call abort
+    if (.not. ieee_is_normal(dx1)) call abort
+    if (dx1 /= 0) call abort
+
+    dx1 = ieee_value(dx1, ieee_negative_zero)
+    if (.not. ieee_is_finite(dx1)) call abort
+    if (ieee_is_nan(dx1)) call abort
+    if (.not. ieee_is_negative(dx1)) call abort
+    if (.not. ieee_is_normal(dx1)) call abort
+    if (dx1 /= 0) call abort
+
+  end if
+
+end
diff --git a/gcc/testsuite/gfortran.dg/ieee/ieee_5.f90 b/gcc/testsuite/gfortran.dg/ieee/ieee_5.f90
new file mode 100644
index 0000000..4ef1525
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/ieee_5.f90
@@ -0,0 +1,34 @@
+! { dg-do run }
+
+  use :: ieee_arithmetic
+  implicit none
+
+  logical mode
+
+  ! Test IEEE_SET_UNDERFLOW_MODE, IEEE_GET_UNDERFLOW_MODE,
+  ! and IEEE_SUPPORT_UNDERFLOW_CONTROL
+  !
+  ! We don't have any targets where this is supported yet, so
+  ! we just check these subroutines are present.
+
+  if (ieee_support_underflow_control() &
+      .or. ieee_support_underflow_control(0.)) then
+
+    call ieee_get_underflow_mode(mode)
+    call ieee_set_underflow_mode(.false.)
+    call ieee_set_underflow_mode(.true.)
+    call ieee_set_underflow_mode(mode)
+
+  end if
+
+  if (ieee_support_underflow_control() &
+      .or. ieee_support_underflow_control(0.d0)) then
+
+    call ieee_get_underflow_mode(mode)
+    call ieee_set_underflow_mode(.false.)
+    call ieee_set_underflow_mode(.true.)
+    call ieee_set_underflow_mode(mode)
+
+  end if
+
+end
diff --git a/gcc/testsuite/gfortran.dg/ieee/ieee_6.f90 b/gcc/testsuite/gfortran.dg/ieee/ieee_6.f90
new file mode 100644
index 0000000..a9a9517
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/ieee_6.f90
@@ -0,0 +1,78 @@
+! { dg-do run }
+!
+! This test will fail on older x86_64 glibc (< 2.20), due to this bug:
+! https://sourceware.org/bugzilla/show_bug.cgi?id=16198
+! We usually won't see it anyway, because on such systems x86_64 assembly
+! (libgfortran/config/fpu-387.h) is used.
+!
+  use :: ieee_arithmetic
+  implicit none
+
+  type(ieee_status_type) :: s1, s2
+  logical :: flags(5), halt(5)
+  type(ieee_round_type) :: mode
+  real :: x
+
+  ! Test IEEE_GET_STATUS and IEEE_SET_STATUS
+
+  call ieee_set_flag(ieee_all, .false.)
+  call ieee_set_rounding_mode(ieee_down)
+  call ieee_set_halting_mode(ieee_all, .false.)
+
+  call ieee_get_status(s1)
+  call ieee_set_status(s1)
+
+  call ieee_get_flag(ieee_all, flags)
+  if (any(flags)) call abort
+  call ieee_get_rounding_mode(mode)
+  if (mode /= ieee_down) call abort
+  call ieee_get_halting_mode(ieee_all, halt)
+  if (any(halt)) call abort
+
+  call ieee_set_rounding_mode(ieee_to_zero)
+  call ieee_set_flag(ieee_underflow, .true.)
+  call ieee_set_halting_mode(ieee_overflow, .true.)
+  x = -1
+  x = sqrt(x)
+  if (.not. ieee_is_nan(x)) call abort
+
+  call ieee_get_status(s2)
+
+  call ieee_get_flag(ieee_all, flags)
+  if (.not. (all(flags .eqv. [.false.,.false.,.true.,.true.,.false.]) &
+             .or. all(flags .eqv. [.false.,.false.,.true.,.false.,.false.]))) call abort
+  call ieee_get_rounding_mode(mode)
+  if (mode /= ieee_to_zero) call abort
+  call ieee_get_halting_mode(ieee_all, halt)
+  if ((.not. halt(1)) .or. any(halt(2:))) call abort
+
+  call ieee_set_status(s2)
+
+  call ieee_get_flag(ieee_all, flags)
+  if (.not. (all(flags .eqv. [.false.,.false.,.true.,.true.,.false.]) &
+             .or. all(flags .eqv. [.false.,.false.,.true.,.false.,.false.]))) call abort
+  call ieee_get_rounding_mode(mode)
+  if (mode /= ieee_to_zero) call abort
+  call ieee_get_halting_mode(ieee_all, halt)
+  if ((.not. halt(1)) .or. any(halt(2:))) call abort
+
+  call ieee_set_status(s1)
+
+  call ieee_get_flag(ieee_all, flags)
+  if (any(flags)) call abort
+  call ieee_get_rounding_mode(mode)
+  if (mode /= ieee_down) call abort
+  call ieee_get_halting_mode(ieee_all, halt)
+  if (any(halt)) call abort
+
+  call ieee_set_status(s2)
+
+  call ieee_get_flag(ieee_all, flags)
+  if (.not. (all(flags .eqv. [.false.,.false.,.true.,.true.,.false.]) &
+             .or. all(flags .eqv. [.false.,.false.,.true.,.false.,.false.]))) call abort
+  call ieee_get_rounding_mode(mode)
+  if (mode /= ieee_to_zero) call abort
+  call ieee_get_halting_mode(ieee_all, halt)
+  if ((.not. halt(1)) .or. any(halt(2:))) call abort
+
+end
diff --git a/gcc/testsuite/gfortran.dg/ieee/ieee_7.f90 b/gcc/testsuite/gfortran.dg/ieee/ieee_7.f90
new file mode 100644
index 0000000..a66e905
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/ieee_7.f90
@@ -0,0 +1,34 @@
+! { dg-do run }
+
+  use :: ieee_arithmetic
+  implicit none
+
+  ! Test IEEE_SELECTED_REAL_KIND in specification expressions
+
+  integer(kind=ieee_selected_real_kind()) :: i1
+  integer(kind=ieee_selected_real_kind(10)) :: i2
+  integer(kind=ieee_selected_real_kind(10,10)) :: i3
+  integer(kind=ieee_selected_real_kind(10,10,2)) :: i4
+
+  ! Test IEEE_SELECTED_REAL_KIND
+
+  if (ieee_support_datatype(0.)) then
+    if (ieee_selected_real_kind() /= kind(0.)) call abort
+    if (ieee_selected_real_kind(0) /= kind(0.)) call abort
+    if (ieee_selected_real_kind(0,0) /= kind(0.)) call abort
+    if (ieee_selected_real_kind(0,0,2) /= kind(0.)) call abort
+  end if
+
+  if (ieee_support_datatype(0.d0)) then
+    if (ieee_selected_real_kind(precision(0.)+1) /= kind(0.d0)) call abort
+    if (ieee_selected_real_kind(precision(0.),range(0.)+1) /= kind(0.d0)) call abort
+    if (ieee_selected_real_kind(precision(0.)+1,range(0.)+1) /= kind(0.d0)) call abort
+    if (ieee_selected_real_kind(precision(0.)+1,range(0.)+1,2) /= kind(0.d0)) call abort
+  end if
+
+  if (ieee_selected_real_kind(0,0,3) /= -5) call abort
+  if (ieee_selected_real_kind(precision(0.d0)+1) /= -1) call abort
+  if (ieee_selected_real_kind(0,range(0.d0)+1) /= -2) call abort
+  if (ieee_selected_real_kind(precision(0.d0)+1,range(0.d0)+1) /= -3) call abort
+
+end
diff --git a/gcc/testsuite/gfortran.dg/ieee/ieee_rounding_1.f90 b/gcc/testsuite/gfortran.dg/ieee/ieee_rounding_1.f90
new file mode 100644
index 0000000..e6bf612
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/ieee/ieee_rounding_1.f90
@@ -0,0 +1,151 @@
+! { dg-do run }
+
+  use, intrinsic :: ieee_features, only : ieee_rounding
+  use, intrinsic :: ieee_arithmetic
+  implicit none
+
+  interface check_equal
+    procedure check_equal_float, check_equal_double
+  end interface
+
+  interface check_not_equal
+    procedure check_not_equal_float, check_not_equal_double
+  end interface
+
+  interface divide
+    procedure divide_float, divide_double
+  end interface
+
+  real :: sx1, sx2, sx3
+  double precision :: dx1, dx2, dx3
+  type(ieee_round_type) :: mode
+
+  ! We should support at least C float and C double types
+  if (ieee_support_rounding(ieee_nearest)) then
+    if (.not. ieee_support_rounding(ieee_nearest, 0.)) call abort
+    if (.not. ieee_support_rounding(ieee_nearest, 0.d0)) call abort
+  end if
+
+  ! The initial rounding mode should probably be NEAREST
+  ! (at least on the platforms we currently support)
+  if (ieee_support_rounding(ieee_nearest, 0.)) then
+    call ieee_get_rounding_mode (mode)
+    if (mode /= ieee_nearest) call abort
+  end if
+
+
+  if (ieee_support_rounding(ieee_up, sx1) .and. &
+      ieee_support_rounding(ieee_down, sx1) .and. &
+      ieee_support_rounding(ieee_nearest, sx1) .and. &
+      ieee_support_rounding(ieee_to_zero, sx1)) then
+
+    sx1 = 1
+    sx2 = 3
+    sx1 = divide(sx1, sx2, ieee_up)
+
+    sx3 = 1
+    sx2 = 3
+    sx3 = divide(sx3, sx2, ieee_down)
+    call check_not_equal(sx1, sx3)
+    call check_equal(sx3, nearest(sx1, -1.))
+    call check_equal(sx1, nearest(sx3,  1.))
+
+    call check_equal(1./3., divide(1., 3., ieee_nearest))
+    call check_equal(-1./3., divide(-1., 3., ieee_nearest))
+
+    call check_equal(divide(3., 7., ieee_to_zero), &
+                    divide(3., 7., ieee_down))
+    call check_equal(divide(-3., 7., ieee_to_zero), &
+                    divide(-3., 7., ieee_up))
+
+  end if
+
+  if (ieee_support_rounding(ieee_up, dx1) .and. &
+      ieee_support_rounding(ieee_down, dx1) .and. &
+      ieee_support_rounding(ieee_nearest, dx1) .and. &
+      ieee_support_rounding(ieee_to_zero, dx1)) then
+
+    dx1 = 1
+    dx2 = 3
+    dx1 = divide(dx1, dx2, ieee_up)
+
+    dx3 = 1
+    dx2 = 3
+    dx3 = divide(dx3, dx2, ieee_down)
+    call check_not_equal(dx1, dx3)
+    call check_equal(dx3, nearest(dx1, -1.d0))
+    call check_equal(dx1, nearest(dx3,  1.d0))
+
+    call check_equal(1.d0/3.d0, divide(1.d0, 3.d0, ieee_nearest))
+    call check_equal(-1.d0/3.d0, divide(-1.d0, 3.d0, ieee_nearest))
+
+    call check_equal(divide(3.d0, 7.d0, ieee_to_zero), &
+                    divide(3.d0, 7.d0, ieee_down))
+    call check_equal(divide(-3.d0, 7.d0, ieee_to_zero), &
+                    divide(-3.d0, 7.d0, ieee_up))
+
+  end if
+
+contains
+
+  real function divide_float (x, y, rounding) result(res)
+    use, intrinsic :: ieee_arithmetic
+    real, intent(in) :: x, y
+    type(ieee_round_type), intent(in) :: rounding
+    type(ieee_round_type) :: old
+
+    call ieee_get_rounding_mode (old)
+    call ieee_set_rounding_mode (rounding)
+
+    res = x / y
+
+    call ieee_set_rounding_mode (old)
+  end function
+
+  double precision function divide_double (x, y, rounding) result(res)
+    use, intrinsic :: ieee_arithmetic
+    double precision, intent(in) :: x, y
+    type(ieee_round_type), intent(in) :: rounding
+    type(ieee_round_type) :: old
+
+    call ieee_get_rounding_mode (old)
+    call ieee_set_rounding_mode (rounding)
+
+    res = x / y
+
+    call ieee_set_rounding_mode (old)
+  end function
+
+  subroutine check_equal_float (x, y)
+    real, intent(in) :: x, y
+    if (x /= y) then
+      print *, x, y
+      call abort
+    end if
+  end subroutine
+
+  subroutine check_equal_double (x, y)
+    double precision, intent(in) :: x, y
+    if (x /= y) then
+      print *, x, y
+      call abort
+    end if
+  end subroutine
+
+  subroutine check_not_equal_float (x, y)
+    real, intent(in) :: x, y
+    if (x == y) then
+      print *, x, y
+      call abort
+    end if
+  end subroutine
+
+  subroutine check_not_equal_double (x, y)
+    double precision, intent(in) :: x, y
+    if (x == y) then
+      print *, x, y
+      call abort
+    end if
+  end subroutine
+
+end
diff --git a/gcc/testsuite/lib/target-supports.exp b/gcc/testsuite/lib/target-supports.exp
index 83a8167..9b31a65 100644
--- a/gcc/testsuite/lib/target-supports.exp
+++ b/gcc/testsuite/lib/target-supports.exp
@@ -1110,6 +1110,20 @@ proc check_effective_target_fortran_real_16 { } {
 }
 
 
+# Return 1 if the target supports Fortran's IEEE modules,
+# 0 otherwise.
+#
+# When the target name changes, replace the cached result.
+
+proc check_effective_target_fortran_ieee { flags } {
+    return [check_no_compiler_messages fortran_ieee executable {
+	! Fortran
+	use, intrinsic :: ieee_features
+	end
+    } $flags ]
+}
+
+
 # Return 1 if the target supports SQRT for the largest floating-point
 # type. (Some targets lack the libm support for this FP type.)
 # On most targets, this check effectively checks either whether sqrtl is