obtained gcc-4.6.4.tar.bz2 from upstream website;upstream

verified gcc-4.6.4.tar.bz2.sig;
imported gcc-4.6.4 source tree from verified upstream tarball.

downloading a git-generated archive based on the 'upstream' tag
should provide you with a source tree that is binary identical
to the one extracted from the above tarball.

if you have obtained the source via the command 'git clone',
however, do note that line-endings of files in your working
directory might differ from line-endings of the respective
files in the upstream repository.

48 files changed, 6216 insertions, 0 deletions

diff --git a/libgfortran/m4/all.m4 b/libgfortran/m4/all.m4
new file mode 100644
index 000000000..cc1742078
--- /dev/null
+++ b/libgfortran/m4/all.m4
@@ -0,0 +1,44 @@
+`/* Implementation of the ALL intrinsic
+   Copyright 2002, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction_logical.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(1,
+`  /* Return true only if all the elements are set.  */
+  result = 1;',
+`  if (! *src)
+    {
+      result = 0;
+      break;
+    }', `')`
+
+#endif'
diff --git a/libgfortran/m4/any.m4 b/libgfortran/m4/any.m4
new file mode 100644
index 000000000..81ec31046
--- /dev/null
+++ b/libgfortran/m4/any.m4
@@ -0,0 +1,44 @@
+`/* Implementation of the ANY intrinsic
+   Copyright 2002, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction_logical.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(0,
+`  result = 0;',
+`  /* Return true if any of the elements are set.  */
+  if (*src)
+    {
+      result = 1;
+      break;
+    }', `')`
+
+#endif'
diff --git a/libgfortran/m4/bessel.m4 b/libgfortran/m4/bessel.m4
new file mode 100644
index 000000000..acef7eaa7
--- /dev/null
+++ b/libgfortran/m4/bessel.m4
@@ -0,0 +1,187 @@
+`/* Implementation of the BESSEL_JN and BESSEL_YN transformational
+   function using a recurrence algorithm.
+   Copyright 2010 Free Software Foundation, Inc.
+   Contributed by Tobias Burnus <burnus@net-b.de>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(`mtype.m4')dnl
+
+mathfunc_macro
+
+`#if defined (HAVE_'rtype_name`)
+
+
+
+#if 'hasmathfunc(jn)`
+extern void bessel_jn_r'rtype_kind` ('rtype` * const restrict ret, int n1,
+				     int n2, 'rtype_name` x);
+export_proto(bessel_jn_r'rtype_kind`);
+
+void
+bessel_jn_r'rtype_kind` ('rtype` * const restrict ret, int n1, int n2, 'rtype_name` x)
+{
+  int i;
+  index_type stride;
+
+  'rtype_name` last1, last2, x2rev;
+
+  stride = GFC_DESCRIPTOR_STRIDE(ret,0);
+
+  if (ret->data == NULL)
+    {
+      size_t size = n2 < n1 ? 0 : n2-n1+1; 
+      GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1);
+      ret->data = internal_malloc_size (sizeof ('rtype_name`) * size);
+      ret->offset = 0;
+    }
+
+  if (unlikely (n2 < n1))
+    return;
+
+  if (unlikely (compile_options.bounds_check)
+      && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1))
+    runtime_error("Incorrect extent in return value of BESSEL_JN "
+		  "(%ld vs. %ld)", (long int) n2-n1,
+		  (long int) GFC_DESCRIPTOR_EXTENT(ret,0));
+
+  stride = GFC_DESCRIPTOR_STRIDE(ret,0);
+
+  if (unlikely (x == 0))
+    {
+      ret->data[0] = 1;
+      for (i = 1; i <= n2-n1; i++)
+        ret->data[i*stride] = 0;
+      return;
+    }
+
+  ret->data = ret->data;
+  last1 = MATHFUNC(jn) (n2, x);
+  ret->data[(n2-n1)*stride] = last1;
+
+  if (n1 == n2)
+    return;
+
+  last2 = MATHFUNC(jn) (n2 - 1, x);
+  ret->data[(n2-n1-1)*stride] = last2;
+
+  if (n1 + 1 == n2)
+    return;
+
+  x2rev = GFC_REAL_'rtype_kind`_LITERAL(2.)/x;
+
+  for (i = n2-n1-2; i >= 0; i--)
+    {
+      ret->data[i*stride] = x2rev * (i+1+n1) * last2 - last1;
+      last1 = last2;
+      last2 = ret->data[i*stride];
+    }
+}
+
+#endif
+
+#if 'hasmathfunc(yn)`
+extern void bessel_yn_r'rtype_kind` ('rtype` * const restrict ret,
+				     int n1, int n2, 'rtype_name` x);
+export_proto(bessel_yn_r'rtype_kind`);
+
+void
+bessel_yn_r'rtype_kind` ('rtype` * const restrict ret, int n1, int n2,
+			 'rtype_name` x)
+{
+  int i;
+  index_type stride;
+
+  'rtype_name` last1, last2, x2rev;
+
+  stride = GFC_DESCRIPTOR_STRIDE(ret,0);
+
+  if (ret->data == NULL)
+    {
+      size_t size = n2 < n1 ? 0 : n2-n1+1; 
+      GFC_DIMENSION_SET(ret->dim[0], 0, size-1, 1);
+      ret->data = internal_malloc_size (sizeof ('rtype_name`) * size);
+      ret->offset = 0;
+    }
+
+  if (unlikely (n2 < n1))
+    return;
+
+  if (unlikely (compile_options.bounds_check)
+      && GFC_DESCRIPTOR_EXTENT(ret,0) != (n2-n1+1))
+    runtime_error("Incorrect extent in return value of BESSEL_JN "
+		  "(%ld vs. %ld)", (long int) n2-n1,
+		  (long int) GFC_DESCRIPTOR_EXTENT(ret,0));
+
+  stride = GFC_DESCRIPTOR_STRIDE(ret,0);
+
+  if (unlikely (x == 0))
+    {
+      for (i = 0; i <= n2-n1; i++)
+#if defined('rtype_name`_INFINITY)
+        ret->data[i*stride] = -'rtype_name`_INFINITY;
+#else
+        ret->data[i*stride] = -'rtype_name`_HUGE;
+#endif
+      return;
+    }
+
+  ret->data = ret->data;
+  last1 = MATHFUNC(yn) (n1, x);
+  ret->data[0] = last1;
+
+  if (n1 == n2)
+    return;
+
+  last2 = MATHFUNC(yn) (n1 + 1, x);
+  ret->data[1*stride] = last2;
+
+  if (n1 + 1 == n2)
+    return;
+
+  x2rev = GFC_REAL_'rtype_kind`_LITERAL(2.)/x;
+
+  for (i = 2; i <= n1+n2; i++)
+    {
+#if defined('rtype_name`_INFINITY)
+      if (unlikely (last2 == -'rtype_name`_INFINITY))
+	{
+	  ret->data[i*stride] = -'rtype_name`_INFINITY;
+	}
+      else
+#endif
+	{
+	  ret->data[i*stride] = x2rev * (i-1+n1) * last2 - last1;
+	  last1 = last2;
+	  last2 = ret->data[i*stride];
+	}
+    }
+}
+#endif
+
+#endif'
+
diff --git a/libgfortran/m4/count.m4 b/libgfortran/m4/count.m4
new file mode 100644
index 000000000..5998c20e6
--- /dev/null
+++ b/libgfortran/m4/count.m4
@@ -0,0 +1,40 @@
+`/* Implementation of the COUNT intrinsic
+   Copyright 2002, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction_logical.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(0,
+`  result = 0;',
+`  if (*src)
+    result++;', `')`
+
+#endif'
diff --git a/libgfortran/m4/cshift0.m4 b/libgfortran/m4/cshift0.m4
new file mode 100644
index 000000000..0c5e0158e
--- /dev/null
+++ b/libgfortran/m4/cshift0.m4
@@ -0,0 +1,172 @@
+`/* Helper function for cshift functions.
+   Copyright 2008, 2009 Free Software Foundation, Inc.
+   Contributed by Thomas Koenig <tkoenig@gcc.gnu.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+void
+cshift0_'rtype_code` ('rtype` *ret, const 'rtype` *array, ssize_t shift,
+		     int which)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type roffset;
+  'rtype_name` *rptr;
+
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type soffset;
+  const 'rtype_name` *sptr;
+
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type dim;
+  index_type len;
+  index_type n;
+
+  which = which - 1;
+  sstride[0] = 0;
+  rstride[0] = 0;
+
+  extent[0] = 1;
+  count[0] = 0;
+  n = 0;
+  /* Initialized for avoiding compiler warnings.  */
+  roffset = 1;
+  soffset = 1;
+  len = 0;
+
+  for (dim = 0; dim < GFC_DESCRIPTOR_RANK (array); dim++)
+    {
+      if (dim == which)
+        {
+          roffset = GFC_DESCRIPTOR_STRIDE(ret,dim);
+          if (roffset == 0)
+            roffset = 1;
+          soffset = GFC_DESCRIPTOR_STRIDE(array,dim);
+          if (soffset == 0)
+            soffset = 1;
+          len = GFC_DESCRIPTOR_EXTENT(array,dim);
+        }
+      else
+        {
+          count[n] = 0;
+          extent[n] = GFC_DESCRIPTOR_EXTENT(array,dim);
+          rstride[n] = GFC_DESCRIPTOR_STRIDE(ret,dim);
+          sstride[n] = GFC_DESCRIPTOR_STRIDE(array,dim);
+          n++;
+        }
+    }
+  if (sstride[0] == 0)
+    sstride[0] = 1;
+  if (rstride[0] == 0)
+    rstride[0] = 1;
+
+  dim = GFC_DESCRIPTOR_RANK (array);
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+  rptr = ret->data;
+  sptr = array->data;
+
+  shift = len == 0 ? 0 : shift % (ssize_t)len;
+  if (shift < 0)
+    shift += len;
+
+  while (rptr)
+    {
+      /* Do the shift for this dimension.  */
+
+      /* If elements are contiguous, perform the operation
+	 in two block moves.  */
+      if (soffset == 1 && roffset == 1)
+	{
+	  size_t len1 = shift * sizeof ('rtype_name`);
+	  size_t len2 = (len - shift) * sizeof ('rtype_name`);
+	  memcpy (rptr, sptr + shift, len2);
+	  memcpy (rptr + (len - shift), sptr, len1);
+	}
+      else
+	{
+	  /* Otherwise, we will have to perform the copy one element at
+	     a time.  */
+	  'rtype_name` *dest = rptr;
+	  const 'rtype_name` *src = &sptr[shift * soffset];
+
+	  for (n = 0; n < len - shift; n++)
+	    {
+	      *dest = *src;
+	      dest += roffset;
+	      src += soffset;
+	    }
+	  for (src = sptr, n = 0; n < shift; n++)
+	    {
+	      *dest = *src;
+	      dest += roffset;
+	      src += soffset;
+	    }
+	}
+
+      /* Advance to the next section.  */
+      rptr += rstride0;
+      sptr += sstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          rptr -= rstride[n] * extent[n];
+          sptr -= sstride[n] * extent[n];
+          n++;
+          if (n >= dim - 1)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              rptr += rstride[n];
+              sptr += sstride[n];
+            }
+        }
+    }
+
+  return;
+}
+
+#endif'
diff --git a/libgfortran/m4/cshift1.m4 b/libgfortran/m4/cshift1.m4
new file mode 100644
index 000000000..49a4f7340
--- /dev/null
+++ b/libgfortran/m4/cshift1.m4
@@ -0,0 +1,274 @@
+`/* Implementation of the CSHIFT intrinsic
+   Copyright 2003, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Feng Wang <wf_cs@yahoo.com>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Ligbfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'atype_name`)
+
+static void
+cshift1 (gfc_array_char * const restrict ret, 
+	const gfc_array_char * const restrict array,
+	const 'atype` * const restrict h, 
+	const 'atype_name` * const restrict pwhich)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type roffset;
+  char *rptr;
+  char *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type soffset;
+  const char *sptr;
+  const char *src;
+  /* h.* indicates the shift array.  */
+  index_type hstride[GFC_MAX_DIMENSIONS];
+  index_type hstride0;
+  const 'atype_name` *hptr;
+
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type dim;
+  index_type len;
+  index_type n;
+  int which;
+  'atype_name` sh;
+  index_type arraysize;
+  index_type size;
+
+  if (pwhich)
+    which = *pwhich - 1;
+  else
+    which = 0;
+
+  if (which < 0 || (which + 1) > GFC_DESCRIPTOR_RANK (array))
+    runtime_error ("Argument ''`DIM''` is out of range in call to ''`CSHIFT''`");
+
+  size = GFC_DESCRIPTOR_SIZE(array);
+
+  arraysize = size0 ((array_t *)array);
+
+  if (ret->data == NULL)
+    {
+      int i;
+
+      ret->data = internal_malloc_size (size * arraysize);
+      ret->offset = 0;
+      ret->dtype = array->dtype;
+      for (i = 0; i < GFC_DESCRIPTOR_RANK (array); i++)
+        {
+	  index_type ub, str;
+
+          ub = GFC_DESCRIPTOR_EXTENT(array,i) - 1;
+
+          if (i == 0)
+            str = 1;
+          else
+	    str = GFC_DESCRIPTOR_EXTENT(ret,i-1) *
+	      GFC_DESCRIPTOR_STRIDE(ret,i-1);
+
+	  GFC_DIMENSION_SET(ret->dim[i], 0, ub, str);
+        }
+    }
+  else if (unlikely (compile_options.bounds_check))
+    {
+      bounds_equal_extents ((array_t *) ret, (array_t *) array,
+				 "return value", "CSHIFT");
+    }
+
+  if (unlikely (compile_options.bounds_check))
+    {
+      bounds_reduced_extents ((array_t *) h, (array_t *) array, which,
+      			      "SHIFT argument", "CSHIFT");
+    }
+
+  if (arraysize == 0)
+    return;
+
+  extent[0] = 1;
+  count[0] = 0;
+  n = 0;
+
+  /* Initialized for avoiding compiler warnings.  */
+  roffset = size;
+  soffset = size;
+  len = 0;
+
+  for (dim = 0; dim < GFC_DESCRIPTOR_RANK (array); dim++)
+    {
+      if (dim == which)
+        {
+          roffset = GFC_DESCRIPTOR_STRIDE_BYTES(ret,dim);
+          if (roffset == 0)
+            roffset = size;
+          soffset = GFC_DESCRIPTOR_STRIDE_BYTES(array,dim);
+          if (soffset == 0)
+            soffset = size;
+          len = GFC_DESCRIPTOR_EXTENT(array,dim);
+        }
+      else
+        {
+          count[n] = 0;
+          extent[n] = GFC_DESCRIPTOR_EXTENT(array,dim);
+          rstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(ret,dim);
+          sstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(array,dim);
+
+          hstride[n] = GFC_DESCRIPTOR_STRIDE(h,n);
+          n++;
+        }
+    }
+  if (sstride[0] == 0)
+    sstride[0] = size;
+  if (rstride[0] == 0)
+    rstride[0] = size;
+  if (hstride[0] == 0)
+    hstride[0] = 1;
+
+  dim = GFC_DESCRIPTOR_RANK (array);
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+  hstride0 = hstride[0];
+  rptr = ret->data;
+  sptr = array->data;
+  hptr = h->data;
+
+  while (rptr)
+    {
+      /* Do the shift for this dimension.  */
+      sh = *hptr;
+      sh = (div (sh, len)).rem;
+      if (sh < 0)
+        sh += len;
+
+      src = &sptr[sh * soffset];
+      dest = rptr;
+
+      for (n = 0; n < len; n++)
+        {
+          memcpy (dest, src, size);
+          dest += roffset;
+          if (n == len - sh - 1)
+            src = sptr;
+          else
+            src += soffset;
+        }
+
+      /* Advance to the next section.  */
+      rptr += rstride0;
+      sptr += sstride0;
+      hptr += hstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          rptr -= rstride[n] * extent[n];
+          sptr -= sstride[n] * extent[n];
+	  hptr -= hstride[n] * extent[n];
+          n++;
+          if (n >= dim - 1)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              rptr += rstride[n];
+              sptr += sstride[n];
+	      hptr += hstride[n];
+            }
+        }
+    }
+}
+
+void cshift1_'atype_kind` (gfc_array_char * const restrict, 
+	const gfc_array_char * const restrict,
+	const 'atype` * const restrict, 
+	const 'atype_name` * const restrict);
+export_proto(cshift1_'atype_kind`);
+
+void
+cshift1_'atype_kind` (gfc_array_char * const restrict ret,
+	const gfc_array_char * const restrict array,
+	const 'atype` * const restrict h, 
+	const 'atype_name` * const restrict pwhich)
+{
+  cshift1 (ret, array, h, pwhich);
+}
+
+
+void cshift1_'atype_kind`_char (gfc_array_char * const restrict ret, 
+	GFC_INTEGER_4,
+	const gfc_array_char * const restrict array,
+	const 'atype` * const restrict h, 
+	const 'atype_name` * const restrict pwhich,
+	GFC_INTEGER_4);
+export_proto(cshift1_'atype_kind`_char);
+
+void
+cshift1_'atype_kind`_char (gfc_array_char * const restrict ret,
+	GFC_INTEGER_4 ret_length __attribute__((unused)),
+	const gfc_array_char * const restrict array,
+	const 'atype` * const restrict h, 
+	const 'atype_name` * const restrict pwhich,
+	GFC_INTEGER_4 array_length __attribute__((unused)))
+{
+  cshift1 (ret, array, h, pwhich);
+}
+
+
+void cshift1_'atype_kind`_char4 (gfc_array_char * const restrict ret, 
+	GFC_INTEGER_4,
+	const gfc_array_char * const restrict array,
+	const 'atype` * const restrict h, 
+	const 'atype_name` * const restrict pwhich,
+	GFC_INTEGER_4);
+export_proto(cshift1_'atype_kind`_char4);
+
+void
+cshift1_'atype_kind`_char4 (gfc_array_char * const restrict ret,
+	GFC_INTEGER_4 ret_length __attribute__((unused)),
+	const gfc_array_char * const restrict array,
+	const 'atype` * const restrict h, 
+	const 'atype_name` * const restrict pwhich,
+	GFC_INTEGER_4 array_length __attribute__((unused)))
+{
+  cshift1 (ret, array, h, pwhich);
+}
+
+#endif'
diff --git a/libgfortran/m4/eoshift1.m4 b/libgfortran/m4/eoshift1.m4
new file mode 100644
index 000000000..b99b3d1a0
--- /dev/null
+++ b/libgfortran/m4/eoshift1.m4
@@ -0,0 +1,319 @@
+`/* Implementation of the EOSHIFT intrinsic
+   Copyright 2002, 2005, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'atype_name`)
+
+static void
+eoshift1 (gfc_array_char * const restrict ret, 
+	const gfc_array_char * const restrict array, 
+	const 'atype` * const restrict h,
+	const char * const restrict pbound, 
+	const 'atype_name` * const restrict pwhich, 
+	const char * filler, index_type filler_len)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type roffset;
+  char *rptr;
+  char * restrict dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type soffset;
+  const char *sptr;
+  const char *src;
+  /* h.* indicates the shift array.  */
+  index_type hstride[GFC_MAX_DIMENSIONS];
+  index_type hstride0;
+  const 'atype_name` *hptr;
+
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type dim;
+  index_type len;
+  index_type n;
+  index_type size;
+  index_type arraysize;
+  int which;
+  'atype_name` sh;
+  'atype_name` delta;
+
+  /* The compiler cannot figure out that these are set, initialize
+     them to avoid warnings.  */
+  len = 0;
+  soffset = 0;
+  roffset = 0;
+
+  size = GFC_DESCRIPTOR_SIZE(array);
+
+  if (pwhich)
+    which = *pwhich - 1;
+  else
+    which = 0;
+
+  extent[0] = 1;
+  count[0] = 0;
+
+  arraysize = size0 ((array_t *) array);
+  if (ret->data == NULL)
+    {
+      int i;
+
+      ret->offset = 0;
+      ret->dtype = array->dtype;
+      for (i = 0; i < GFC_DESCRIPTOR_RANK (array); i++)
+        {
+	  index_type ub, str;
+
+	  ub = GFC_DESCRIPTOR_EXTENT(array,i) - 1;
+
+          if (i == 0)
+            str = 1;
+          else
+            str = GFC_DESCRIPTOR_EXTENT(ret,i-1)
+	      * GFC_DESCRIPTOR_STRIDE(ret,i-1);
+
+	  GFC_DIMENSION_SET(ret->dim[i], 0, ub, str);
+
+        }
+      if (arraysize > 0)
+	ret->data = internal_malloc_size (size * arraysize);
+      else
+	ret->data = internal_malloc_size (1);
+
+    }
+  else if (unlikely (compile_options.bounds_check))
+    {
+      bounds_equal_extents ((array_t *) ret, (array_t *) array,
+				 "return value", "EOSHIFT");
+    }
+
+  if (unlikely (compile_options.bounds_check))
+    {
+      bounds_reduced_extents ((array_t *) h, (array_t *) array, which,
+      			      "SHIFT argument", "EOSHIFT");
+    }
+
+  if (arraysize == 0)
+    return;
+
+  n = 0;
+  for (dim = 0; dim < GFC_DESCRIPTOR_RANK (array); dim++)
+    {
+      if (dim == which)
+        {
+          roffset = GFC_DESCRIPTOR_STRIDE_BYTES(ret,dim);
+          if (roffset == 0)
+            roffset = size;
+          soffset = GFC_DESCRIPTOR_STRIDE_BYTES(array,dim);
+          if (soffset == 0)
+            soffset = size;
+          len = GFC_DESCRIPTOR_EXTENT(array,dim);
+        }
+      else
+        {
+          count[n] = 0;
+          extent[n] = GFC_DESCRIPTOR_EXTENT(array,dim);
+          rstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(ret,dim);
+          sstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(array,dim);
+
+          hstride[n] = GFC_DESCRIPTOR_STRIDE(h,n);
+          n++;
+        }
+    }
+  if (sstride[0] == 0)
+    sstride[0] = size;
+  if (rstride[0] == 0)
+    rstride[0] = size;
+  if (hstride[0] == 0)
+    hstride[0] = 1;
+
+  dim = GFC_DESCRIPTOR_RANK (array);
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+  hstride0 = hstride[0];
+  rptr = ret->data;
+  sptr = array->data;
+  hptr = h->data;
+
+  while (rptr)
+    {
+      /* Do the shift for this dimension.  */
+      sh = *hptr;
+      if (( sh >= 0 ? sh : -sh ) > len)
+	{
+	  delta = len;
+	  sh = len;
+	}
+      else
+	delta = (sh >= 0) ? sh: -sh;
+
+      if (sh > 0)
+        {
+          src = &sptr[delta * soffset];
+          dest = rptr;
+        }
+      else
+        {
+          src = sptr;
+          dest = &rptr[delta * roffset];
+        }
+      for (n = 0; n < len - delta; n++)
+        {
+          memcpy (dest, src, size);
+          dest += roffset;
+          src += soffset;
+        }
+      if (sh < 0)
+        dest = rptr;
+      n = delta;
+
+      if (pbound)
+	while (n--)
+	  {
+	    memcpy (dest, pbound, size);
+	    dest += roffset;
+	  }
+      else
+	while (n--)
+	  {
+	    index_type i;
+
+	    if (filler_len == 1)
+	      memset (dest, filler[0], size);
+	    else
+	      for (i = 0; i < size; i += filler_len)
+		memcpy (&dest[i], filler, filler_len);
+
+	    dest += roffset;
+	  }
+
+      /* Advance to the next section.  */
+      rptr += rstride0;
+      sptr += sstride0;
+      hptr += hstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          rptr -= rstride[n] * extent[n];
+          sptr -= sstride[n] * extent[n];
+	  hptr -= hstride[n] * extent[n];
+          n++;
+          if (n >= dim - 1)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              rptr += rstride[n];
+              sptr += sstride[n];
+	      hptr += hstride[n];
+            }
+        }
+    }
+}
+
+void eoshift1_'atype_kind` (gfc_array_char * const restrict, 
+	const gfc_array_char * const restrict,
+	const 'atype` * const restrict, const char * const restrict, 
+	const 'atype_name` * const restrict);
+export_proto(eoshift1_'atype_kind`);
+
+void
+eoshift1_'atype_kind` (gfc_array_char * const restrict ret, 
+	const gfc_array_char * const restrict array,
+	const 'atype` * const restrict h, 
+	const char * const restrict pbound,
+	const 'atype_name` * const restrict pwhich)
+{
+  eoshift1 (ret, array, h, pbound, pwhich, "\0", 1);
+}
+
+
+void eoshift1_'atype_kind`_char (gfc_array_char * const restrict, 
+	GFC_INTEGER_4,
+	const gfc_array_char * const restrict, 
+	const 'atype` * const restrict,
+	const char * const restrict, 
+	const 'atype_name` * const restrict,
+	GFC_INTEGER_4, GFC_INTEGER_4);
+export_proto(eoshift1_'atype_kind`_char);
+
+void
+eoshift1_'atype_kind`_char (gfc_array_char * const restrict ret,
+	GFC_INTEGER_4 ret_length __attribute__((unused)),
+	const gfc_array_char * const restrict array, 
+	const 'atype` * const restrict h,
+	const char *  const restrict pbound, 
+	const 'atype_name` * const restrict pwhich,
+	GFC_INTEGER_4 array_length __attribute__((unused)),
+	GFC_INTEGER_4 bound_length __attribute__((unused)))
+{
+  eoshift1 (ret, array, h, pbound, pwhich, " ", 1);
+}
+
+
+void eoshift1_'atype_kind`_char4 (gfc_array_char * const restrict, 
+	GFC_INTEGER_4,
+	const gfc_array_char * const restrict, 
+	const 'atype` * const restrict,
+	const char * const restrict, 
+	const 'atype_name` * const restrict,
+	GFC_INTEGER_4, GFC_INTEGER_4);
+export_proto(eoshift1_'atype_kind`_char4);
+
+void
+eoshift1_'atype_kind`_char4 (gfc_array_char * const restrict ret,
+	GFC_INTEGER_4 ret_length __attribute__((unused)),
+	const gfc_array_char * const restrict array, 
+	const 'atype` * const restrict h,
+	const char *  const restrict pbound, 
+	const 'atype_name` * const restrict pwhich,
+	GFC_INTEGER_4 array_length __attribute__((unused)),
+	GFC_INTEGER_4 bound_length __attribute__((unused)))
+{
+  static const gfc_char4_t space = (unsigned char) ''` ''`;
+  eoshift1 (ret, array, h, pbound, pwhich,
+	    (const char *) &space, sizeof (gfc_char4_t));
+}
+
+#endif'
diff --git a/libgfortran/m4/eoshift3.m4 b/libgfortran/m4/eoshift3.m4
new file mode 100644
index 000000000..d1b10b56c
--- /dev/null
+++ b/libgfortran/m4/eoshift3.m4
@@ -0,0 +1,337 @@
+`/* Implementation of the EOSHIFT intrinsic
+   Copyright 2002, 2005, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'atype_name`)
+
+static void
+eoshift3 (gfc_array_char * const restrict ret, 
+	const gfc_array_char * const restrict array, 
+	const 'atype` * const restrict h,
+	const gfc_array_char * const restrict bound, 
+	const 'atype_name` * const restrict pwhich,
+	const char * filler, index_type filler_len)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type roffset;
+  char *rptr;
+  char * restrict dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type soffset;
+  const char *sptr;
+  const char *src;
+  /* h.* indicates the shift array.  */
+  index_type hstride[GFC_MAX_DIMENSIONS];
+  index_type hstride0;
+  const 'atype_name` *hptr;
+  /* b.* indicates the bound array.  */
+  index_type bstride[GFC_MAX_DIMENSIONS];
+  index_type bstride0;
+  const char *bptr;
+
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type dim;
+  index_type len;
+  index_type n;
+  index_type size;
+  index_type arraysize;
+  int which;
+  'atype_name` sh;
+  'atype_name` delta;
+
+  /* The compiler cannot figure out that these are set, initialize
+     them to avoid warnings.  */
+  len = 0;
+  soffset = 0;
+  roffset = 0;
+
+  arraysize = size0 ((array_t *) array);
+  size = GFC_DESCRIPTOR_SIZE(array);
+
+  if (pwhich)
+    which = *pwhich - 1;
+  else
+    which = 0;
+
+  if (ret->data == NULL)
+    {
+      int i;
+
+      ret->offset = 0;
+      ret->dtype = array->dtype;
+      for (i = 0; i < GFC_DESCRIPTOR_RANK (array); i++)
+        {
+	  index_type ub, str;
+
+	  ub = GFC_DESCRIPTOR_EXTENT(array,i) - 1;
+
+          if (i == 0)
+            str = 1;
+          else
+            str = GFC_DESCRIPTOR_EXTENT(ret,i-1)
+	      * GFC_DESCRIPTOR_STRIDE(ret,i-1);
+
+	  GFC_DIMENSION_SET(ret->dim[i], 0, ub, str);
+
+        }
+      if (arraysize > 0)
+	ret->data = internal_malloc_size (size * arraysize);
+      else
+	ret->data = internal_malloc_size (1);
+
+    }
+  else if (unlikely (compile_options.bounds_check))
+    {
+      bounds_equal_extents ((array_t *) ret, (array_t *) array,
+				 "return value", "EOSHIFT");
+    }
+
+  if (unlikely (compile_options.bounds_check))
+    {
+      bounds_reduced_extents ((array_t *) h, (array_t *) array, which,
+      			      "SHIFT argument", "EOSHIFT");
+    }
+
+  if (arraysize == 0)
+    return;
+
+  extent[0] = 1;
+  count[0] = 0;
+  n = 0;
+  for (dim = 0; dim < GFC_DESCRIPTOR_RANK (array); dim++)
+    {
+      if (dim == which)
+        {
+          roffset = GFC_DESCRIPTOR_STRIDE_BYTES(ret,dim);
+          if (roffset == 0)
+            roffset = size;
+          soffset = GFC_DESCRIPTOR_STRIDE_BYTES(array,dim);
+          if (soffset == 0)
+            soffset = size;
+          len = GFC_DESCRIPTOR_EXTENT(array,dim);
+        }
+      else
+        {
+          count[n] = 0;
+          extent[n] = GFC_DESCRIPTOR_EXTENT(array,dim);
+          rstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(ret,dim);
+          sstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(array,dim);
+
+          hstride[n] = GFC_DESCRIPTOR_STRIDE(h,n);
+          if (bound)
+            bstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(bound,n);
+          else
+            bstride[n] = 0;
+          n++;
+        }
+    }
+  if (sstride[0] == 0)
+    sstride[0] = size;
+  if (rstride[0] == 0)
+    rstride[0] = size;
+  if (hstride[0] == 0)
+    hstride[0] = 1;
+  if (bound && bstride[0] == 0)
+    bstride[0] = size;
+
+  dim = GFC_DESCRIPTOR_RANK (array);
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+  hstride0 = hstride[0];
+  bstride0 = bstride[0];
+  rptr = ret->data;
+  sptr = array->data;
+  hptr = h->data;
+  if (bound)
+    bptr = bound->data;
+  else
+    bptr = NULL;
+
+  while (rptr)
+    {
+      /* Do the shift for this dimension.  */
+      sh = *hptr;
+      if (( sh >= 0 ? sh : -sh ) > len)
+	{
+	  delta = len;
+	  sh = len;
+	}
+      else
+	delta = (sh >= 0) ? sh: -sh;
+
+      if (sh > 0)
+        {
+          src = &sptr[delta * soffset];
+          dest = rptr;
+        }
+      else
+        {
+          src = sptr;
+          dest = &rptr[delta * roffset];
+        }
+      for (n = 0; n < len - delta; n++)
+        {
+          memcpy (dest, src, size);
+          dest += roffset;
+          src += soffset;
+        }
+      if (sh < 0)
+        dest = rptr;
+      n = delta;
+
+      if (bptr)
+	while (n--)
+	  {
+	    memcpy (dest, bptr, size);
+	    dest += roffset;
+	  }
+      else
+	while (n--)
+	  {
+	    index_type i;
+
+	    if (filler_len == 1)
+	      memset (dest, filler[0], size);
+	    else
+	      for (i = 0; i < size; i += filler_len)
+		memcpy (&dest[i], filler, filler_len);
+
+	    dest += roffset;
+	  }
+
+      /* Advance to the next section.  */
+      rptr += rstride0;
+      sptr += sstride0;
+      hptr += hstride0;
+      bptr += bstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          rptr -= rstride[n] * extent[n];
+          sptr -= sstride[n] * extent[n];
+	  hptr -= hstride[n] * extent[n];
+          bptr -= bstride[n] * extent[n];
+          n++;
+          if (n >= dim - 1)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              rptr += rstride[n];
+              sptr += sstride[n];
+	      hptr += hstride[n];
+              bptr += bstride[n];
+            }
+        }
+    }
+}
+
+extern void eoshift3_'atype_kind` (gfc_array_char * const restrict, 
+	const gfc_array_char * const restrict,
+	const 'atype` * const restrict, 
+	const gfc_array_char * const restrict,
+	const 'atype_name` *);
+export_proto(eoshift3_'atype_kind`);
+
+void
+eoshift3_'atype_kind` (gfc_array_char * const restrict ret, 
+	const gfc_array_char * const restrict array,
+	const 'atype` * const restrict h, 
+	const gfc_array_char * const restrict bound,
+	const 'atype_name` * const restrict pwhich)
+{
+  eoshift3 (ret, array, h, bound, pwhich, "\0", 1);
+}
+
+
+extern void eoshift3_'atype_kind`_char (gfc_array_char * const restrict, 
+	GFC_INTEGER_4,
+	const gfc_array_char * const restrict,
+	const 'atype` * const restrict,
+	const gfc_array_char * const restrict,
+	const 'atype_name` * const restrict, 
+	GFC_INTEGER_4, GFC_INTEGER_4);
+export_proto(eoshift3_'atype_kind`_char);
+
+void
+eoshift3_'atype_kind`_char (gfc_array_char * const restrict ret,
+	GFC_INTEGER_4 ret_length __attribute__((unused)),
+	const gfc_array_char * const restrict array, 
+	const 'atype` *  const restrict h,
+	const gfc_array_char * const restrict bound,
+	const 'atype_name` * const restrict pwhich,
+	GFC_INTEGER_4 array_length __attribute__((unused)),
+	GFC_INTEGER_4 bound_length __attribute__((unused)))
+{
+  eoshift3 (ret, array, h, bound, pwhich, " ", 1);
+}
+
+
+extern void eoshift3_'atype_kind`_char4 (gfc_array_char * const restrict, 
+	GFC_INTEGER_4,
+	const gfc_array_char * const restrict,
+	const 'atype` * const restrict,
+	const gfc_array_char * const restrict,
+	const 'atype_name` * const restrict, 
+	GFC_INTEGER_4, GFC_INTEGER_4);
+export_proto(eoshift3_'atype_kind`_char4);
+
+void
+eoshift3_'atype_kind`_char4 (gfc_array_char * const restrict ret,
+	GFC_INTEGER_4 ret_length __attribute__((unused)),
+	const gfc_array_char * const restrict array, 
+	const 'atype` *  const restrict h,
+	const gfc_array_char * const restrict bound,
+	const 'atype_name` * const restrict pwhich,
+	GFC_INTEGER_4 array_length __attribute__((unused)),
+	GFC_INTEGER_4 bound_length __attribute__((unused)))
+{
+  static const gfc_char4_t space = (unsigned char) ''` ''`;
+  eoshift3 (ret, array, h, bound, pwhich,
+	    (const char *) &space, sizeof (gfc_char4_t));
+}
+
+#endif'
diff --git a/libgfortran/m4/exponent.m4 b/libgfortran/m4/exponent.m4
new file mode 100644
index 000000000..043f150fc
--- /dev/null
+++ b/libgfortran/m4/exponent.m4
@@ -0,0 +1,45 @@
+`/* Implementation of the EXPONENT intrinsic
+   Copyright 2003, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Richard Henderson <rth@redhat.com>.
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+mathfunc_macro
+
+`#if defined (HAVE_'real_type`) && 'hasmathfunc(frexp)`
+
+extern GFC_INTEGER_4 exponent_r'kind` ('real_type` s);
+export_proto(exponent_r'kind`);
+
+GFC_INTEGER_4
+exponent_r'kind` ('real_type` s)
+{
+  int ret;
+  MATHFUNC(frexp) (s, &ret);
+  return ret;
+}
+
+#endif'
diff --git a/libgfortran/m4/fraction.m4 b/libgfortran/m4/fraction.m4
new file mode 100644
index 000000000..cca2db980
--- /dev/null
+++ b/libgfortran/m4/fraction.m4
@@ -0,0 +1,44 @@
+`/* Implementation of the FRACTION intrinsic
+   Copyright 2003, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Richard Henderson <rth@redhat.com>.
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+mathfunc_macro
+
+`#if defined (HAVE_'real_type`) && 'hasmathfunc(frexp)`
+
+extern 'real_type` fraction_r'kind` ('real_type` s);
+export_proto(fraction_r'kind`);
+
+'real_type`
+fraction_r'kind` ('real_type` s)
+{
+  int dummy_exp;
+  return MATHFUNC(frexp) (s, &dummy_exp);
+}
+
+#endif'
diff --git a/libgfortran/m4/head.m4 b/libgfortran/m4/head.m4
new file mode 100644
index 000000000..30cdea892
--- /dev/null
+++ b/libgfortran/m4/head.m4
@@ -0,0 +1,25 @@
+`!   Copyright 2002, 2007, 2009 Free Software Foundation, Inc.
+!   Contributed by Paul Brook <paul@nowt.org>
+!
+!This file is part of the GNU Fortran 95 runtime library (libgfortran).
+!
+!GNU libgfortran 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 of the License, or (at your option) any later version.
+
+!GNU libgfortran 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.
+!
+!Under Section 7 of GPL version 3, you are granted additional
+!permissions described in the GCC Runtime Library Exception, version
+!3.1, as published by the Free Software Foundation.
+!
+!You should have received a copy of the GNU General Public License and
+!a copy of the GCC Runtime Library Exception along with this program;
+!see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+!<http://www.gnu.org/licenses/>.
+!
+!This file is machine generated.'
diff --git a/libgfortran/m4/iall.m4 b/libgfortran/m4/iall.m4
new file mode 100644
index 000000000..2e6667e26
--- /dev/null
+++ b/libgfortran/m4/iall.m4
@@ -0,0 +1,46 @@
+`/* Implementation of the IALL intrinsic
+   Copyright 2010 Free Software Foundation, Inc.
+   Contributed by Tobias Burnus <burnus@net-b.de>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(0,
+`  result = ('rtype_name`) -1;',
+`  result &= *src;')
+
+MASKED_ARRAY_FUNCTION(0,
+`  result = 0;',
+`  if (*msrc)
+    result &= *src;')
+
+SCALAR_ARRAY_FUNCTION(0)
+
+#endif
diff --git a/libgfortran/m4/iany.m4 b/libgfortran/m4/iany.m4
new file mode 100644
index 000000000..a17d951e5
--- /dev/null
+++ b/libgfortran/m4/iany.m4
@@ -0,0 +1,46 @@
+`/* Implementation of the IANY intrinsic
+   Copyright 2010 Free Software Foundation, Inc.
+   Contributed by Tobias Burnus <burnus@net-b.de>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(0,
+`  result = 0;',
+`  result |= *src;')
+
+MASKED_ARRAY_FUNCTION(0,
+`  result = 0;',
+`  if (*msrc)
+    result |= *src;')
+
+SCALAR_ARRAY_FUNCTION(0)
+
+#endif
diff --git a/libgfortran/m4/iforeach.m4 b/libgfortran/m4/iforeach.m4
new file mode 100644
index 000000000..14e501c67
--- /dev/null
+++ b/libgfortran/m4/iforeach.m4
@@ -0,0 +1,279 @@
+dnl Support macro file for intrinsic functions.
+dnl Contains the generic sections of the array functions.
+dnl This file is part of the GNU Fortran 95 Runtime Library (libgfortran)
+dnl Distributed under the GNU GPL with exception.  See COPYING for details.
+define(START_FOREACH_FUNCTION,
+`
+extern void name`'rtype_qual`_'atype_code (rtype * const restrict retarray, 
+	atype * const restrict array);
+export_proto(name`'rtype_qual`_'atype_code);
+
+void
+name`'rtype_qual`_'atype_code (rtype * const restrict retarray, 
+	atype * const restrict array)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type dstride;
+  const atype_name *base;
+  rtype_name * restrict dest;
+  index_type rank;
+  index_type n;
+
+  rank = GFC_DESCRIPTOR_RANK (array);
+  if (rank <= 0)
+    runtime_error ("Rank of array needs to be > 0");
+
+  if (retarray->data == NULL)
+    {
+      GFC_DIMENSION_SET(retarray->dim[0], 0, rank-1, 1);
+      retarray->dtype = (retarray->dtype & ~GFC_DTYPE_RANK_MASK) | 1;
+      retarray->offset = 0;
+      retarray->data = internal_malloc_size (sizeof (rtype_name) * rank);
+    }
+  else
+    {
+      if (unlikely (compile_options.bounds_check))
+	bounds_iforeach_return ((array_t *) retarray, (array_t *) array,
+				"u_name");
+    }
+
+  dstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
+  dest = retarray->data;
+  for (n = 0; n < rank; n++)
+    {
+      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
+      count[n] = 0;
+      if (extent[n] <= 0)
+	{
+	  /* Set the return value.  */
+	  for (n = 0; n < rank; n++)
+	    dest[n * dstride] = 0;
+	  return;
+	}
+    }
+
+  base = array->data;
+
+  /* Initialize the return value.  */
+  for (n = 0; n < rank; n++)
+    dest[n * dstride] = 1;
+  {
+')dnl
+define(START_FOREACH_BLOCK,
+`  while (base)
+    {
+      do
+	{
+	  /* Implementation start.  */
+')dnl
+define(FINISH_FOREACH_FUNCTION,
+`	  /* Implementation end.  */
+	  /* Advance to the next element.  */
+	  base += sstride[0];
+	}
+      while (++count[0] != extent[0]);
+      n = 0;
+      do
+	{
+	  /* When we get to the end of a dimension, reset it and increment
+	     the next dimension.  */
+	  count[n] = 0;
+	  /* We could precalculate these products, but this is a less
+	     frequently used path so probably not worth it.  */
+	  base -= sstride[n] * extent[n];
+	  n++;
+	  if (n == rank)
+	    {
+	      /* Break out of the loop.  */
+	      base = NULL;
+	      break;
+	    }
+	  else
+	    {
+	      count[n]++;
+	      base += sstride[n];
+	    }
+	}
+      while (count[n] == extent[n]);
+    }
+  }
+}')dnl
+define(START_MASKED_FOREACH_FUNCTION,
+`
+extern void `m'name`'rtype_qual`_'atype_code (rtype * const restrict, 
+	atype * const restrict, gfc_array_l1 * const restrict);
+export_proto(`m'name`'rtype_qual`_'atype_code);
+
+void
+`m'name`'rtype_qual`_'atype_code (rtype * const restrict retarray, 
+	atype * const restrict array,
+	gfc_array_l1 * const restrict mask)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type mstride[GFC_MAX_DIMENSIONS];
+  index_type dstride;
+  rtype_name *dest;
+  const atype_name *base;
+  GFC_LOGICAL_1 *mbase;
+  int rank;
+  index_type n;
+  int mask_kind;
+
+  rank = GFC_DESCRIPTOR_RANK (array);
+  if (rank <= 0)
+    runtime_error ("Rank of array needs to be > 0");
+
+  if (retarray->data == NULL)
+    {
+      GFC_DIMENSION_SET(retarray->dim[0], 0, rank - 1, 1);
+      retarray->dtype = (retarray->dtype & ~GFC_DTYPE_RANK_MASK) | 1;
+      retarray->offset = 0;
+      retarray->data = internal_malloc_size (sizeof (rtype_name) * rank);
+    }
+  else
+    {
+      if (unlikely (compile_options.bounds_check))
+	{
+
+	  bounds_iforeach_return ((array_t *) retarray, (array_t *) array,
+				  "u_name");
+	  bounds_equal_extents ((array_t *) mask, (array_t *) array,
+				  "MASK argument", "u_name");
+	}
+    }
+
+  mask_kind = GFC_DESCRIPTOR_SIZE (mask);
+
+  mbase = mask->data;
+
+  if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
+#ifdef HAVE_GFC_LOGICAL_16
+      || mask_kind == 16
+#endif
+      )
+    mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
+  else
+    runtime_error ("Funny sized logical array");
+
+  dstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
+  dest = retarray->data;
+  for (n = 0; n < rank; n++)
+    {
+      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
+      mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
+      count[n] = 0;
+      if (extent[n] <= 0)
+	{
+	  /* Set the return value.  */
+	  for (n = 0; n < rank; n++)
+	    dest[n * dstride] = 0;
+	  return;
+	}
+    }
+
+  base = array->data;
+
+  /* Initialize the return value.  */
+  for (n = 0; n < rank; n++)
+    dest[n * dstride] = 0;
+  {
+')dnl
+define(START_MASKED_FOREACH_BLOCK, `START_FOREACH_BLOCK')dnl
+define(FINISH_MASKED_FOREACH_FUNCTION,
+`	  /* Implementation end.  */
+	  /* Advance to the next element.  */
+	  base += sstride[0];
+	  mbase += mstride[0];
+	}
+      while (++count[0] != extent[0]);
+      n = 0;
+      do
+	{
+	  /* When we get to the end of a dimension, reset it and increment
+	     the next dimension.  */
+	  count[n] = 0;
+	  /* We could precalculate these products, but this is a less
+	     frequently used path so probably not worth it.  */
+	  base -= sstride[n] * extent[n];
+	  mbase -= mstride[n] * extent[n];
+	  n++;
+	  if (n == rank)
+	    {
+	      /* Break out of the loop.  */
+	      base = NULL;
+	      break;
+	    }
+	  else
+	    {
+	      count[n]++;
+	      base += sstride[n];
+	      mbase += mstride[n];
+	    }
+	}
+      while (count[n] == extent[n]);
+    }
+  }
+}')dnl
+define(FOREACH_FUNCTION,
+`START_FOREACH_FUNCTION
+$1
+START_FOREACH_BLOCK
+$2
+FINISH_FOREACH_FUNCTION')dnl
+define(MASKED_FOREACH_FUNCTION,
+`START_MASKED_FOREACH_FUNCTION
+$1
+START_MASKED_FOREACH_BLOCK
+$2
+FINISH_MASKED_FOREACH_FUNCTION')dnl
+define(SCALAR_FOREACH_FUNCTION,
+`
+extern void `s'name`'rtype_qual`_'atype_code (rtype * const restrict, 
+	atype * const restrict, GFC_LOGICAL_4 *);
+export_proto(`s'name`'rtype_qual`_'atype_code);
+
+void
+`s'name`'rtype_qual`_'atype_code (rtype * const restrict retarray, 
+	atype * const restrict array,
+	GFC_LOGICAL_4 * mask)
+{
+  index_type rank;
+  index_type dstride;
+  index_type n;
+  rtype_name *dest;
+
+  if (*mask)
+    {
+      name`'rtype_qual`_'atype_code (retarray, array);
+      return;
+    }
+
+  rank = GFC_DESCRIPTOR_RANK (array);
+
+  if (rank <= 0)
+    runtime_error ("Rank of array needs to be > 0");
+
+  if (retarray->data == NULL)
+    {
+      GFC_DIMENSION_SET(retarray->dim[0], 0, rank-1, 1);
+      retarray->dtype = (retarray->dtype & ~GFC_DTYPE_RANK_MASK) | 1;
+      retarray->offset = 0;
+      retarray->data = internal_malloc_size (sizeof (rtype_name) * rank);
+    }
+  else if (unlikely (compile_options.bounds_check))
+    {
+       bounds_iforeach_return ((array_t *) retarray, (array_t *) array,
+			       "u_name");
+    }
+
+  dstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
+  dest = retarray->data;
+  for (n = 0; n<rank; n++)
+    dest[n * dstride] = $1 ;
+}')dnl
diff --git a/libgfortran/m4/ifunction.m4 b/libgfortran/m4/ifunction.m4
new file mode 100644
index 000000000..68e1c5e18
--- /dev/null
+++ b/libgfortran/m4/ifunction.m4
@@ -0,0 +1,505 @@
+dnl Support macro file for intrinsic functions.
+dnl Contains the generic sections of the array functions.
+dnl This file is part of the GNU Fortran Runtime Library (libgfortran)
+dnl Distributed under the GNU GPL with exception.  See COPYING for details.
+dnl
+dnl Pass the implementation for a single section as the parameter to
+dnl {MASK_}ARRAY_FUNCTION.
+dnl The variables base, delta, and len describe the input section.
+dnl For masked section the mask is described by mbase and mdelta.
+dnl These should not be modified. The result should be stored in *dest.
+dnl The names count, extent, sstride, dstride, base, dest, rank, dim
+dnl retarray, array, pdim and mstride should not be used.
+dnl The variable n is declared as index_type and may be used.
+dnl Other variable declarations may be placed at the start of the code,
+dnl The types of the array parameter and the return value are
+dnl atype_name and rtype_name respectively.
+dnl Execution should be allowed to continue to the end of the block.
+dnl You should not return or break from the inner loop of the implementation.
+dnl Care should also be taken to avoid using the names defined in iparm.m4
+define(START_ARRAY_FUNCTION,
+`
+extern void name`'rtype_qual`_'atype_code (rtype * const restrict, 
+	atype * const restrict, const index_type * const restrict);
+export_proto(name`'rtype_qual`_'atype_code);
+
+void
+name`'rtype_qual`_'atype_code (rtype * const restrict retarray, 
+	atype * const restrict array, 
+	const index_type * const restrict pdim)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type dstride[GFC_MAX_DIMENSIONS];
+  const atype_name * restrict base;
+  rtype_name * restrict dest;
+  index_type rank;
+  index_type n;
+  index_type len;
+  index_type delta;
+  index_type dim;
+  int continue_loop;
+
+  /* Make dim zero based to avoid confusion.  */
+  dim = (*pdim) - 1;
+  rank = GFC_DESCRIPTOR_RANK (array) - 1;
+
+  len = GFC_DESCRIPTOR_EXTENT(array,dim);
+  if (len < 0)
+    len = 0;
+  delta = GFC_DESCRIPTOR_STRIDE(array,dim);
+
+  for (n = 0; n < dim; n++)
+    {
+      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
+
+      if (extent[n] < 0)
+	extent[n] = 0;
+    }
+  for (n = dim; n < rank; n++)
+    {
+      sstride[n] = GFC_DESCRIPTOR_STRIDE(array, n + 1);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
+
+      if (extent[n] < 0)
+	extent[n] = 0;
+    }
+
+  if (retarray->data == NULL)
+    {
+      size_t alloc_size, str;
+
+      for (n = 0; n < rank; n++)
+	{
+	  if (n == 0)
+	    str = 1;
+	  else
+	    str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
+
+	  GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
+
+	}
+
+      retarray->offset = 0;
+      retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
+
+      alloc_size = sizeof (rtype_name) * GFC_DESCRIPTOR_STRIDE(retarray,rank-1)
+    		   * extent[rank-1];
+
+      retarray->data = internal_malloc_size (alloc_size);
+      if (alloc_size == 0)
+	{
+	  /* Make sure we have a zero-sized array.  */
+	  GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
+	  return;
+
+	}
+    }
+  else
+    {
+      if (rank != GFC_DESCRIPTOR_RANK (retarray))
+	runtime_error ("rank of return array incorrect in"
+		       " u_name intrinsic: is %ld, should be %ld",
+		       (long int) (GFC_DESCRIPTOR_RANK (retarray)),
+		       (long int) rank);
+
+      if (unlikely (compile_options.bounds_check))
+	bounds_ifunction_return ((array_t *) retarray, extent,
+				 "return value", "u_name");
+    }
+
+  for (n = 0; n < rank; n++)
+    {
+      count[n] = 0;
+      dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
+      if (extent[n] <= 0)
+	return;
+    }
+
+  base = array->data;
+  dest = retarray->data;
+
+  continue_loop = 1;
+  while (continue_loop)
+    {
+      const atype_name * restrict src;
+      rtype_name result;
+      src = base;
+      {
+')dnl
+define(START_ARRAY_BLOCK,
+`	if (len <= 0)
+	  *dest = '$1`;
+	else
+	  {
+	    for (n = 0; n < len; n++, src += delta)
+	      {
+')dnl
+define(FINISH_ARRAY_FUNCTION,
+`	      }
+	    '$1`
+	    *dest = result;
+	  }
+      }
+      /* Advance to the next element.  */
+      count[0]++;
+      base += sstride[0];
+      dest += dstride[0];
+      n = 0;
+      while (count[n] == extent[n])
+	{
+	  /* When we get to the end of a dimension, reset it and increment
+	     the next dimension.  */
+	  count[n] = 0;
+	  /* We could precalculate these products, but this is a less
+	     frequently used path so probably not worth it.  */
+	  base -= sstride[n] * extent[n];
+	  dest -= dstride[n] * extent[n];
+	  n++;
+	  if (n == rank)
+	    {
+	      /* Break out of the look.  */
+	      continue_loop = 0;
+	      break;
+	    }
+	  else
+	    {
+	      count[n]++;
+	      base += sstride[n];
+	      dest += dstride[n];
+	    }
+	}
+    }
+}')dnl
+define(START_MASKED_ARRAY_FUNCTION,
+`
+extern void `m'name`'rtype_qual`_'atype_code (rtype * const restrict, 
+	atype * const restrict, const index_type * const restrict,
+	gfc_array_l1 * const restrict);
+export_proto(`m'name`'rtype_qual`_'atype_code);
+
+void
+`m'name`'rtype_qual`_'atype_code (rtype * const restrict retarray, 
+	atype * const restrict array, 
+	const index_type * const restrict pdim, 
+	gfc_array_l1 * const restrict mask)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type dstride[GFC_MAX_DIMENSIONS];
+  index_type mstride[GFC_MAX_DIMENSIONS];
+  rtype_name * restrict dest;
+  const atype_name * restrict base;
+  const GFC_LOGICAL_1 * restrict mbase;
+  int rank;
+  int dim;
+  index_type n;
+  index_type len;
+  index_type delta;
+  index_type mdelta;
+  int mask_kind;
+
+  dim = (*pdim) - 1;
+  rank = GFC_DESCRIPTOR_RANK (array) - 1;
+
+  len = GFC_DESCRIPTOR_EXTENT(array,dim);
+  if (len <= 0)
+    return;
+
+  mbase = mask->data;
+
+  mask_kind = GFC_DESCRIPTOR_SIZE (mask);
+
+  if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
+#ifdef HAVE_GFC_LOGICAL_16
+      || mask_kind == 16
+#endif
+      )
+    mbase = GFOR_POINTER_TO_L1 (mbase, mask_kind);
+  else
+    runtime_error ("Funny sized logical array");
+
+  delta = GFC_DESCRIPTOR_STRIDE(array,dim);
+  mdelta = GFC_DESCRIPTOR_STRIDE_BYTES(mask,dim);
+
+  for (n = 0; n < dim; n++)
+    {
+      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
+      mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
+
+      if (extent[n] < 0)
+	extent[n] = 0;
+
+    }
+  for (n = dim; n < rank; n++)
+    {
+      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n + 1);
+      mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask, n + 1);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array, n + 1);
+
+      if (extent[n] < 0)
+	extent[n] = 0;
+    }
+
+  if (retarray->data == NULL)
+    {
+      size_t alloc_size, str;
+
+      for (n = 0; n < rank; n++)
+	{
+	  if (n == 0)
+	    str = 1;
+	  else
+	    str= GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
+
+	  GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
+
+	}
+
+      alloc_size = sizeof (rtype_name) * GFC_DESCRIPTOR_STRIDE(retarray,rank-1)
+    		   * extent[rank-1];
+
+      retarray->offset = 0;
+      retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
+
+      if (alloc_size == 0)
+	{
+	  /* Make sure we have a zero-sized array.  */
+	  GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
+	  return;
+	}
+      else
+	retarray->data = internal_malloc_size (alloc_size);
+
+    }
+  else
+    {
+      if (rank != GFC_DESCRIPTOR_RANK (retarray))
+	runtime_error ("rank of return array incorrect in u_name intrinsic");
+
+      if (unlikely (compile_options.bounds_check))
+	{
+	  bounds_ifunction_return ((array_t *) retarray, extent,
+				   "return value", "u_name");
+	  bounds_equal_extents ((array_t *) mask, (array_t *) array,
+	  			"MASK argument", "u_name");
+	}
+    }
+
+  for (n = 0; n < rank; n++)
+    {
+      count[n] = 0;
+      dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
+      if (extent[n] <= 0)
+	return;
+    }
+
+  dest = retarray->data;
+  base = array->data;
+
+  while (base)
+    {
+      const atype_name * restrict src;
+      const GFC_LOGICAL_1 * restrict msrc;
+      rtype_name result;
+      src = base;
+      msrc = mbase;
+      {
+')dnl
+define(START_MASKED_ARRAY_BLOCK,
+`	if (len <= 0)
+	  *dest = '$1`;
+	else
+	  {
+	    for (n = 0; n < len; n++, src += delta, msrc += mdelta)
+	      {
+')dnl
+define(FINISH_MASKED_ARRAY_FUNCTION,
+`	      }
+	    *dest = result;
+	  }
+      }
+      /* Advance to the next element.  */
+      count[0]++;
+      base += sstride[0];
+      mbase += mstride[0];
+      dest += dstride[0];
+      n = 0;
+      while (count[n] == extent[n])
+	{
+	  /* When we get to the end of a dimension, reset it and increment
+	     the next dimension.  */
+	  count[n] = 0;
+	  /* We could precalculate these products, but this is a less
+	     frequently used path so probably not worth it.  */
+	  base -= sstride[n] * extent[n];
+	  mbase -= mstride[n] * extent[n];
+	  dest -= dstride[n] * extent[n];
+	  n++;
+	  if (n == rank)
+	    {
+	      /* Break out of the look.  */
+	      base = NULL;
+	      break;
+	    }
+	  else
+	    {
+	      count[n]++;
+	      base += sstride[n];
+	      mbase += mstride[n];
+	      dest += dstride[n];
+	    }
+	}
+    }
+}')dnl
+define(SCALAR_ARRAY_FUNCTION,
+`
+extern void `s'name`'rtype_qual`_'atype_code (rtype * const restrict, 
+	atype * const restrict, const index_type * const restrict,
+	GFC_LOGICAL_4 *);
+export_proto(`s'name`'rtype_qual`_'atype_code);
+
+void
+`s'name`'rtype_qual`_'atype_code (rtype * const restrict retarray, 
+	atype * const restrict array, 
+	const index_type * const restrict pdim, 
+	GFC_LOGICAL_4 * mask)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type dstride[GFC_MAX_DIMENSIONS];
+  rtype_name * restrict dest;
+  index_type rank;
+  index_type n;
+  index_type dim;
+
+
+  if (*mask)
+    {
+      name`'rtype_qual`_'atype_code (retarray, array, pdim);
+      return;
+    }
+  /* Make dim zero based to avoid confusion.  */
+  dim = (*pdim) - 1;
+  rank = GFC_DESCRIPTOR_RANK (array) - 1;
+
+  for (n = 0; n < dim; n++)
+    {
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
+
+      if (extent[n] <= 0)
+	extent[n] = 0;
+    }
+
+  for (n = dim; n < rank; n++)
+    {
+      extent[n] =
+	GFC_DESCRIPTOR_EXTENT(array,n + 1);
+
+      if (extent[n] <= 0)
+	extent[n] = 0;
+    }
+
+  if (retarray->data == NULL)
+    {
+      size_t alloc_size, str;
+
+      for (n = 0; n < rank; n++)
+	{
+	  if (n == 0)
+	    str = 1;
+	  else
+	    str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
+
+	  GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
+
+	}
+
+      retarray->offset = 0;
+      retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
+
+      alloc_size = sizeof (rtype_name) * GFC_DESCRIPTOR_STRIDE(retarray,rank-1)
+    		   * extent[rank-1];
+
+      if (alloc_size == 0)
+	{
+	  /* Make sure we have a zero-sized array.  */
+	  GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
+	  return;
+	}
+      else
+	retarray->data = internal_malloc_size (alloc_size);
+    }
+  else
+    {
+      if (rank != GFC_DESCRIPTOR_RANK (retarray))
+	runtime_error ("rank of return array incorrect in"
+		       " u_name intrinsic: is %ld, should be %ld",
+		       (long int) (GFC_DESCRIPTOR_RANK (retarray)),
+		       (long int) rank);
+
+      if (unlikely (compile_options.bounds_check))
+	{
+	  for (n=0; n < rank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n);
+	      if (extent[n] != ret_extent)
+		runtime_error ("Incorrect extent in return value of"
+			       " u_name intrinsic in dimension %ld:"
+			       " is %ld, should be %ld", (long int) n + 1,
+			       (long int) ret_extent, (long int) extent[n]);
+	    }
+	}
+    }
+
+  for (n = 0; n < rank; n++)
+    {
+      count[n] = 0;
+      dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
+    }
+
+  dest = retarray->data;
+
+  while(1)
+    {
+      *dest = '$1`;
+      count[0]++;
+      dest += dstride[0];
+      n = 0;
+      while (count[n] == extent[n])
+	{
+	  /* When we get to the end of a dimension, reset it and increment
+	     the next dimension.  */
+	  count[n] = 0;
+	  /* We could precalculate these products, but this is a less
+	     frequently used path so probably not worth it.  */
+	  dest -= dstride[n] * extent[n];
+	  n++;
+	  if (n == rank)
+	    return;
+	  else
+	    {
+	      count[n]++;
+	      dest += dstride[n];
+	    }
+      	}
+    }
+}')dnl
+define(ARRAY_FUNCTION,
+`START_ARRAY_FUNCTION
+$2
+START_ARRAY_BLOCK($1)
+$3
+FINISH_ARRAY_FUNCTION($4)')dnl
+define(MASKED_ARRAY_FUNCTION,
+`START_MASKED_ARRAY_FUNCTION
+$2
+START_MASKED_ARRAY_BLOCK($1)
+$3
+FINISH_MASKED_ARRAY_FUNCTION')dnl
diff --git a/libgfortran/m4/ifunction_logical.m4 b/libgfortran/m4/ifunction_logical.m4
new file mode 100644
index 000000000..e72f1d8a8
--- /dev/null
+++ b/libgfortran/m4/ifunction_logical.m4
@@ -0,0 +1,208 @@
+dnl Support macro file for intrinsic functions.
+dnl Contains the generic sections of the array functions.
+dnl This file is part of the GNU Fortran 95 Runtime Library (libgfortran)
+dnl Distributed under the GNU GPL with exception.  See COPYING for details.
+dnl
+dnl Pass the implementation for a single section as the parameter to
+dnl {MASK_}ARRAY_FUNCTION.
+dnl The variables base, delta, and len describe the input section.
+dnl For masked section the mask is described by mbase and mdelta.
+dnl These should not be modified. The result should be stored in *dest.
+dnl The names count, extent, sstride, dstride, base, dest, rank, dim
+dnl retarray, array, pdim and mstride should not be used.
+dnl The variable n is declared as index_type and may be used.
+dnl Other variable declarations may be placed at the start of the code,
+dnl The types of the array parameter and the return value are
+dnl atype_name and rtype_name respectively.
+dnl Execution should be allowed to continue to the end of the block.
+dnl You should not return or break from the inner loop of the implementation.
+dnl Care should also be taken to avoid using the names defined in iparm.m4
+define(START_ARRAY_FUNCTION,
+`
+extern void name`'rtype_qual`_'atype_code (rtype * const restrict, 
+	gfc_array_l1 * const restrict, const index_type * const restrict);
+export_proto(name`'rtype_qual`_'atype_code);
+
+void
+name`'rtype_qual`_'atype_code (rtype * const restrict retarray, 
+	gfc_array_l1 * const restrict array, 
+	const index_type * const restrict pdim)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type dstride[GFC_MAX_DIMENSIONS];
+  const GFC_LOGICAL_1 * restrict base;
+  rtype_name * restrict dest;
+  index_type rank;
+  index_type n;
+  index_type len;
+  index_type delta;
+  index_type dim;
+  int src_kind;
+  int continue_loop;
+
+  /* Make dim zero based to avoid confusion.  */
+  dim = (*pdim) - 1;
+  rank = GFC_DESCRIPTOR_RANK (array) - 1;
+
+  src_kind = GFC_DESCRIPTOR_SIZE (array);
+
+  len = GFC_DESCRIPTOR_EXTENT(array,dim);
+  if (len < 0)
+    len = 0;
+
+  delta = GFC_DESCRIPTOR_STRIDE_BYTES(array,dim);
+
+  for (n = 0; n < dim; n++)
+    {
+      sstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(array,n);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
+
+      if (extent[n] < 0)
+	extent[n] = 0;
+    }
+  for (n = dim; n < rank; n++)
+    {
+      sstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(array,n + 1);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n + 1);
+
+      if (extent[n] < 0)
+	extent[n] = 0;
+    }
+
+  if (retarray->data == NULL)
+    {
+      size_t alloc_size, str;
+
+      for (n = 0; n < rank; n++)
+        {
+          if (n == 0)
+            str = 1;
+          else
+            str = GFC_DESCRIPTOR_STRIDE(retarray,n-1) * extent[n-1];
+
+	  GFC_DIMENSION_SET(retarray->dim[n], 0, extent[n] - 1, str);
+
+        }
+
+      retarray->offset = 0;
+      retarray->dtype = (array->dtype & ~GFC_DTYPE_RANK_MASK) | rank;
+
+      alloc_size = sizeof (rtype_name) * GFC_DESCRIPTOR_STRIDE(retarray,rank-1)
+    		   * extent[rank-1];
+
+      if (alloc_size == 0)
+	{
+	  /* Make sure we have a zero-sized array.  */
+	  GFC_DIMENSION_SET(retarray->dim[0], 0, -1, 1);
+	  return;
+	}
+      else
+	retarray->data = internal_malloc_size (alloc_size);
+    }
+  else
+    {
+      if (rank != GFC_DESCRIPTOR_RANK (retarray))
+	runtime_error ("rank of return array incorrect in"
+		       " u_name intrinsic: is %ld, should be %ld",
+		       (long int) GFC_DESCRIPTOR_RANK (retarray),
+		       (long int) rank);
+
+      if (unlikely (compile_options.bounds_check))
+	{
+	  for (n=0; n < rank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,n);
+	      if (extent[n] != ret_extent)
+		runtime_error ("Incorrect extent in return value of"
+			       " u_name intrinsic in dimension %d:"
+			       " is %ld, should be %ld", (int) n + 1,
+			       (long int) ret_extent, (long int) extent[n]);
+	    }
+	}
+    }
+
+  for (n = 0; n < rank; n++)
+    {
+      count[n] = 0;
+      dstride[n] = GFC_DESCRIPTOR_STRIDE(retarray,n);
+      if (extent[n] <= 0)
+	return;
+    }
+
+  base = array->data;
+
+  if (src_kind == 1 || src_kind == 2 || src_kind == 4 || src_kind == 8
+#ifdef HAVE_GFC_LOGICAL_16
+      || src_kind == 16
+#endif
+    )
+    {
+      if (base)
+	base = GFOR_POINTER_TO_L1 (base, src_kind);
+    }
+  else
+    internal_error (NULL, "Funny sized logical array in u_name intrinsic");
+
+  dest = retarray->data;
+
+  continue_loop = 1;
+  while (continue_loop)
+    {
+      const GFC_LOGICAL_1 * restrict src;
+      rtype_name result;
+      src = base;
+      {
+')dnl
+define(START_ARRAY_BLOCK,
+`        if (len <= 0)
+	  *dest = '$1`;
+	else
+	  {
+	    for (n = 0; n < len; n++, src += delta)
+	      {
+')dnl
+define(FINISH_ARRAY_FUNCTION,
+    `          }
+	    *dest = result;
+	  }
+      }
+      /* Advance to the next element.  */
+      count[0]++;
+      base += sstride[0];
+      dest += dstride[0];
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          base -= sstride[n] * extent[n];
+          dest -= dstride[n] * extent[n];
+          n++;
+          if (n == rank)
+            {
+              /* Break out of the look.  */
+              continue_loop = 0;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              base += sstride[n];
+              dest += dstride[n];
+            }
+        }
+    }
+}')dnl
+define(ARRAY_FUNCTION,
+`START_ARRAY_FUNCTION
+$2
+START_ARRAY_BLOCK($1)
+$3
+FINISH_ARRAY_FUNCTION')dnl
diff --git a/libgfortran/m4/in_pack.m4 b/libgfortran/m4/in_pack.m4
new file mode 100644
index 000000000..a4337aad8
--- /dev/null
+++ b/libgfortran/m4/in_pack.m4
@@ -0,0 +1,122 @@
+`/* Helper function for repacking arrays.
+   Copyright 2003, 2006, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+/* Allocates a block of memory with internal_malloc if the array needs
+   repacking.  */
+'
+dnl The kind (ie size) is used to name the function for logicals, integers
+dnl and reals.  For complex, it's c4 or c8.
+rtype_name` *
+internal_pack_'rtype_ccode` ('rtype` * source)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type stride[GFC_MAX_DIMENSIONS];
+  index_type stride0;
+  index_type dim;
+  index_type ssize;
+  const 'rtype_name` *src;
+  'rtype_name` * restrict dest;
+  'rtype_name` *destptr;
+  int n;
+  int packed;
+
+  /* TODO: Investigate how we can figure out if this is a temporary
+     since the stride=0 thing has been removed from the frontend.  */
+
+  dim = GFC_DESCRIPTOR_RANK (source);
+  ssize = 1;
+  packed = 1;
+  for (n = 0; n < dim; n++)
+    {
+      count[n] = 0;
+      stride[n] = GFC_DESCRIPTOR_STRIDE(source,n);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(source,n);
+      if (extent[n] <= 0)
+        {
+          /* Do nothing.  */
+          packed = 1;
+          break;
+        }
+
+      if (ssize != stride[n])
+        packed = 0;
+
+      ssize *= extent[n];
+    }
+
+  if (packed)
+    return source->data;
+
+  /* Allocate storage for the destination.  */
+  destptr = ('rtype_name` *)internal_malloc_size (ssize * sizeof ('rtype_name`));
+  dest = destptr;
+  src = source->data;
+  stride0 = stride[0];
+
+
+  while (src)
+    {
+      /* Copy the data.  */
+      *(dest++) = *src;
+      /* Advance to the next element.  */
+      src += stride0;
+      count[0]++;
+      /* Advance to the next source element.  */
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          src -= stride[n] * extent[n];
+          n++;
+          if (n == dim)
+            {
+              src = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              src += stride[n];
+            }
+        }
+    }
+  return destptr;
+}
+
+#endif
+'
diff --git a/libgfortran/m4/in_unpack.m4 b/libgfortran/m4/in_unpack.m4
new file mode 100644
index 000000000..661c54e1d
--- /dev/null
+++ b/libgfortran/m4/in_unpack.m4
@@ -0,0 +1,110 @@
+`/* Helper function for repacking arrays.
+   Copyright 2003, 2006, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)'
+
+dnl Only the kind (ie size) is used to name the function for integers,
+dnl reals and logicals.  For complex, it's c4 and c8.
+`void
+internal_unpack_'rtype_ccode` ('rtype` * d, const 'rtype_name` * src)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type stride[GFC_MAX_DIMENSIONS];
+  index_type stride0;
+  index_type dim;
+  index_type dsize;
+  'rtype_name` * restrict dest;
+  int n;
+
+  dest = d->data;
+  if (src == dest || !src)
+    return;
+
+  dim = GFC_DESCRIPTOR_RANK (d);
+  dsize = 1;
+  for (n = 0; n < dim; n++)
+    {
+      count[n] = 0;
+      stride[n] = GFC_DESCRIPTOR_STRIDE(d,n);
+      extent[n] = GFC_DESCRIPTOR_EXTENT(d,n);
+      if (extent[n] <= 0)
+	return;
+
+      if (dsize == stride[n])
+	dsize *= extent[n];
+      else
+	dsize = 0;
+    }
+
+  if (dsize != 0)
+    {
+      memcpy (dest, src, dsize * sizeof ('rtype_name`));
+      return;
+    }
+
+  stride0 = stride[0];
+
+  while (dest)
+    {
+      /* Copy the data.  */
+      *dest = *(src++);
+      /* Advance to the next element.  */
+      dest += stride0;
+      count[0]++;
+      /* Advance to the next source element.  */
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          dest -= stride[n] * extent[n];
+          n++;
+          if (n == dim)
+            {
+              dest = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              dest += stride[n];
+            }
+        }
+    }
+}
+
+#endif
+'
diff --git a/libgfortran/m4/iparity.m4 b/libgfortran/m4/iparity.m4
new file mode 100644
index 000000000..78dbc3dd4
--- /dev/null
+++ b/libgfortran/m4/iparity.m4
@@ -0,0 +1,46 @@
+`/* Implementation of the IPARITY intrinsic
+   Copyright 2010 Free Software Foundation, Inc.
+   Contributed by Tobias Burnus <burnus@net-b.de>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(0,
+`  result = 0;',
+`  result ^= *src;')
+
+MASKED_ARRAY_FUNCTION(0,
+`  result = 0;',
+`  if (*msrc)
+    result ^= *src;')
+
+SCALAR_ARRAY_FUNCTION(0)
+
+#endif
diff --git a/libgfortran/m4/iparm.m4 b/libgfortran/m4/iparm.m4
new file mode 100644
index 000000000..4ff624741
--- /dev/null
+++ b/libgfortran/m4/iparm.m4
@@ -0,0 +1,37 @@
+dnl Support macro file for intrinsic functions.
+dnl Works out all the function types from the filename.
+dnl This file is part of the GNU Fortran 95 Runtime Library (libgfortran)
+dnl Distributed under the GNU GPL with exception.  See COPYING for details.
+dnl M4 macro file to get type names from filenames
+define(get_typename2, `GFC_$1_$2')dnl
+define(get_typename, `get_typename2(ifelse($1,i,INTEGER,ifelse($1,r,REAL,ifelse($1,l,LOGICAL,ifelse($1,c,COMPLEX,unknown)))),`$2')')dnl
+define(get_arraytype, `gfc_array_$1$2')dnl
+define(define_type, `dnl
+ifelse(regexp($2,`^[0-9]'),-1,`dnl
+define($1_letter, substr($2, 0, 1))dnl
+define($1_kind, substr($2, 1))dnl
+',`dnl
+define($1_letter,i)dnl
+define($1_kind,$2)dnl
+')dnl
+define($1_code,$1_letter`'$1_kind)dnl
+define($1,get_arraytype($1_letter,$1_kind))dnl
+define($1_name, get_typename($1_letter, $1_kind))')dnl
+dnl
+define_type(atype, regexp(file, `_\(.?[0-9]*\)\.c$', `\1'))dnl
+define(rtype_tmp, regexp(file, `_\(.?[0-9]*\)_[^_]*\.c$', `\1'))dnl
+ifelse(rtype_tmp,,`dnl
+define_type(rtype, atype_code)dnl
+define(rtype_qual,`')dnl
+',`dnl
+define_type(rtype, rtype_tmp)dnl
+define(rtype_qual,`_'rtype_kind)dnl
+')dnl
+define(atype_max, atype_name`_HUGE')dnl
+define(atype_min,ifelse(regexp(file, `_\(.\)[0-9]*\.c$', `\1'),`i',`(-'atype_max`-1)',`-'atype_max))dnl
+define(atype_inf, atype_name`_INFINITY')dnl
+define(atype_nan, atype_name`_QUIET_NAN')dnl
+define(name, regexp(regexp(file, `[^/]*$', `\&'), `^\([^_]*\)_', `\1'))dnl
+define(`upcase', `translit(`$*', `a-z', `A-Z')')dnl
+define(`u_name',`regexp(upcase(name),`\([A-Z]*\)',`\1')')dnl
+define(rtype_ccode,ifelse(rtype_letter,`i',rtype_kind,rtype_code))dnl
diff --git a/libgfortran/m4/matmul.m4 b/libgfortran/m4/matmul.m4
new file mode 100644
index 000000000..bb42f2a6c
--- /dev/null
+++ b/libgfortran/m4/matmul.m4
@@ -0,0 +1,378 @@
+`/* Implementation of the MATMUL intrinsic
+   Copyright 2002, 2005, 2006, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+/* Prototype for the BLAS ?gemm subroutine, a pointer to which can be
+   passed to us by the front-end, in which case we''`ll call it for large
+   matrices.  */
+
+typedef void (*blas_call)(const char *, const char *, const int *, const int *,
+                          const int *, const 'rtype_name` *, const 'rtype_name` *,
+                          const int *, const 'rtype_name` *, const int *,
+                          const 'rtype_name` *, 'rtype_name` *, const int *,
+                          int, int);
+
+/* The order of loops is different in the case of plain matrix
+   multiplication C=MATMUL(A,B), and in the frequent special case where
+   the argument A is the temporary result of a TRANSPOSE intrinsic:
+   C=MATMUL(TRANSPOSE(A),B).  Transposed temporaries are detected by
+   looking at their strides.
+
+   The equivalent Fortran pseudo-code is:
+
+   DIMENSION A(M,COUNT), B(COUNT,N), C(M,N)
+   IF (.NOT.IS_TRANSPOSED(A)) THEN
+     C = 0
+     DO J=1,N
+       DO K=1,COUNT
+         DO I=1,M
+           C(I,J) = C(I,J)+A(I,K)*B(K,J)
+   ELSE
+     DO J=1,N
+       DO I=1,M
+         S = 0
+         DO K=1,COUNT
+           S = S+A(I,K)*B(K,J)
+         C(I,J) = S
+   ENDIF
+*/
+
+/* If try_blas is set to a nonzero value, then the matmul function will
+   see if there is a way to perform the matrix multiplication by a call
+   to the BLAS gemm function.  */
+
+extern void matmul_'rtype_code` ('rtype` * const restrict retarray, 
+	'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm);
+export_proto(matmul_'rtype_code`);
+
+void
+matmul_'rtype_code` ('rtype` * const restrict retarray, 
+	'rtype` * const restrict a, 'rtype` * const restrict b, int try_blas,
+	int blas_limit, blas_call gemm)
+{
+  const 'rtype_name` * restrict abase;
+  const 'rtype_name` * restrict bbase;
+  'rtype_name` * restrict dest;
+
+  index_type rxstride, rystride, axstride, aystride, bxstride, bystride;
+  index_type x, y, n, count, xcount, ycount;
+
+  assert (GFC_DESCRIPTOR_RANK (a) == 2
+          || GFC_DESCRIPTOR_RANK (b) == 2);
+
+/* C[xcount,ycount] = A[xcount, count] * B[count,ycount]
+
+   Either A or B (but not both) can be rank 1:
+
+   o One-dimensional argument A is implicitly treated as a row matrix
+     dimensioned [1,count], so xcount=1.
+
+   o One-dimensional argument B is implicitly treated as a column matrix
+     dimensioned [count, 1], so ycount=1.
+  */
+
+  if (retarray->data == NULL)
+    {
+      if (GFC_DESCRIPTOR_RANK (a) == 1)
+        {
+	  GFC_DIMENSION_SET(retarray->dim[0], 0,
+	                    GFC_DESCRIPTOR_EXTENT(b,1) - 1, 1);
+        }
+      else if (GFC_DESCRIPTOR_RANK (b) == 1)
+        {
+	  GFC_DIMENSION_SET(retarray->dim[0], 0,
+	                    GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
+        }
+      else
+        {
+	  GFC_DIMENSION_SET(retarray->dim[0], 0,
+	                    GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
+
+          GFC_DIMENSION_SET(retarray->dim[1], 0,
+	                    GFC_DESCRIPTOR_EXTENT(b,1) - 1,
+			    GFC_DESCRIPTOR_EXTENT(retarray,0));
+        }
+
+      retarray->data
+	= internal_malloc_size (sizeof ('rtype_name`) * size0 ((array_t *) retarray));
+      retarray->offset = 0;
+    }
+    else if (unlikely (compile_options.bounds_check))
+      {
+	index_type ret_extent, arg_extent;
+
+	if (GFC_DESCRIPTOR_RANK (a) == 1)
+	  {
+	    arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
+	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
+	    if (arg_extent != ret_extent)
+	      runtime_error ("Incorrect extent in return array in"
+			     " MATMUL intrinsic: is %ld, should be %ld",
+			     (long int) ret_extent, (long int) arg_extent);
+	  }
+	else if (GFC_DESCRIPTOR_RANK (b) == 1)
+	  {
+	    arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
+	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
+	    if (arg_extent != ret_extent)
+	      runtime_error ("Incorrect extent in return array in"
+			     " MATMUL intrinsic: is %ld, should be %ld",
+			     (long int) ret_extent, (long int) arg_extent);	    
+	  }
+	else
+	  {
+	    arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
+	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
+	    if (arg_extent != ret_extent)
+	      runtime_error ("Incorrect extent in return array in"
+			     " MATMUL intrinsic for dimension 1:"
+			     " is %ld, should be %ld",
+			     (long int) ret_extent, (long int) arg_extent);
+
+	    arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
+	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
+	    if (arg_extent != ret_extent)
+	      runtime_error ("Incorrect extent in return array in"
+			     " MATMUL intrinsic for dimension 2:"
+			     " is %ld, should be %ld",
+			     (long int) ret_extent, (long int) arg_extent);
+	  }
+      }
+'
+sinclude(`matmul_asm_'rtype_code`.m4')dnl
+`
+  if (GFC_DESCRIPTOR_RANK (retarray) == 1)
+    {
+      /* One-dimensional result may be addressed in the code below
+	 either as a row or a column matrix. We want both cases to
+	 work. */
+      rxstride = rystride = GFC_DESCRIPTOR_STRIDE(retarray,0);
+    }
+  else
+    {
+      rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
+      rystride = GFC_DESCRIPTOR_STRIDE(retarray,1);
+    }
+
+
+  if (GFC_DESCRIPTOR_RANK (a) == 1)
+    {
+      /* Treat it as a a row matrix A[1,count]. */
+      axstride = GFC_DESCRIPTOR_STRIDE(a,0);
+      aystride = 1;
+
+      xcount = 1;
+      count = GFC_DESCRIPTOR_EXTENT(a,0);
+    }
+  else
+    {
+      axstride = GFC_DESCRIPTOR_STRIDE(a,0);
+      aystride = GFC_DESCRIPTOR_STRIDE(a,1);
+
+      count = GFC_DESCRIPTOR_EXTENT(a,1);
+      xcount = GFC_DESCRIPTOR_EXTENT(a,0);
+    }
+
+  if (count != GFC_DESCRIPTOR_EXTENT(b,0))
+    {
+      if (count > 0 || GFC_DESCRIPTOR_EXTENT(b,0) > 0)
+	runtime_error ("dimension of array B incorrect in MATMUL intrinsic");
+    }
+
+  if (GFC_DESCRIPTOR_RANK (b) == 1)
+    {
+      /* Treat it as a column matrix B[count,1] */
+      bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
+
+      /* bystride should never be used for 1-dimensional b.
+	 in case it is we want it to cause a segfault, rather than
+	 an incorrect result. */
+      bystride = 0xDEADBEEF;
+      ycount = 1;
+    }
+  else
+    {
+      bxstride = GFC_DESCRIPTOR_STRIDE(b,0);
+      bystride = GFC_DESCRIPTOR_STRIDE(b,1);
+      ycount = GFC_DESCRIPTOR_EXTENT(b,1);
+    }
+
+  abase = a->data;
+  bbase = b->data;
+  dest = retarray->data;
+
+
+  /* Now that everything is set up, we''`re performing the multiplication
+     itself.  */
+
+#define POW3(x) (((float) (x)) * ((float) (x)) * ((float) (x)))
+
+  if (try_blas && rxstride == 1 && (axstride == 1 || aystride == 1)
+      && (bxstride == 1 || bystride == 1)
+      && (((float) xcount) * ((float) ycount) * ((float) count)
+          > POW3(blas_limit)))
+  {
+    const int m = xcount, n = ycount, k = count, ldc = rystride;
+    const 'rtype_name` one = 1, zero = 0;
+    const int lda = (axstride == 1) ? aystride : axstride,
+              ldb = (bxstride == 1) ? bystride : bxstride;
+
+    if (lda > 0 && ldb > 0 && ldc > 0 && m > 1 && n > 1 && k > 1)
+      {
+        assert (gemm != NULL);
+        gemm (axstride == 1 ? "N" : "T", bxstride == 1 ? "N" : "T", &m, &n, &k,
+              &one, abase, &lda, bbase, &ldb, &zero, dest, &ldc, 1, 1);
+        return;
+      }
+  }
+
+  if (rxstride == 1 && axstride == 1 && bxstride == 1)
+    {
+      const 'rtype_name` * restrict bbase_y;
+      'rtype_name` * restrict dest_y;
+      const 'rtype_name` * restrict abase_n;
+      'rtype_name` bbase_yn;
+
+      if (rystride == xcount)
+	memset (dest, 0, (sizeof ('rtype_name`) * xcount * ycount));
+      else
+	{
+	  for (y = 0; y < ycount; y++)
+	    for (x = 0; x < xcount; x++)
+	      dest[x + y*rystride] = ('rtype_name`)0;
+	}
+
+      for (y = 0; y < ycount; y++)
+	{
+	  bbase_y = bbase + y*bystride;
+	  dest_y = dest + y*rystride;
+	  for (n = 0; n < count; n++)
+	    {
+	      abase_n = abase + n*aystride;
+	      bbase_yn = bbase_y[n];
+	      for (x = 0; x < xcount; x++)
+		{
+		  dest_y[x] += abase_n[x] * bbase_yn;
+		}
+	    }
+	}
+    }
+  else if (rxstride == 1 && aystride == 1 && bxstride == 1)
+    {
+      if (GFC_DESCRIPTOR_RANK (a) != 1)
+	{
+	  const 'rtype_name` *restrict abase_x;
+	  const 'rtype_name` *restrict bbase_y;
+	  'rtype_name` *restrict dest_y;
+	  'rtype_name` s;
+
+	  for (y = 0; y < ycount; y++)
+	    {
+	      bbase_y = &bbase[y*bystride];
+	      dest_y = &dest[y*rystride];
+	      for (x = 0; x < xcount; x++)
+		{
+		  abase_x = &abase[x*axstride];
+		  s = ('rtype_name`) 0;
+		  for (n = 0; n < count; n++)
+		    s += abase_x[n] * bbase_y[n];
+		  dest_y[x] = s;
+		}
+	    }
+	}
+      else
+	{
+	  const 'rtype_name` *restrict bbase_y;
+	  'rtype_name` s;
+
+	  for (y = 0; y < ycount; y++)
+	    {
+	      bbase_y = &bbase[y*bystride];
+	      s = ('rtype_name`) 0;
+	      for (n = 0; n < count; n++)
+		s += abase[n*axstride] * bbase_y[n];
+	      dest[y*rystride] = s;
+	    }
+	}
+    }
+  else if (axstride < aystride)
+    {
+      for (y = 0; y < ycount; y++)
+	for (x = 0; x < xcount; x++)
+	  dest[x*rxstride + y*rystride] = ('rtype_name`)0;
+
+      for (y = 0; y < ycount; y++)
+	for (n = 0; n < count; n++)
+	  for (x = 0; x < xcount; x++)
+	    /* dest[x,y] += a[x,n] * b[n,y] */
+	    dest[x*rxstride + y*rystride] += abase[x*axstride + n*aystride] * bbase[n*bxstride + y*bystride];
+    }
+  else if (GFC_DESCRIPTOR_RANK (a) == 1)
+    {
+      const 'rtype_name` *restrict bbase_y;
+      'rtype_name` s;
+
+      for (y = 0; y < ycount; y++)
+	{
+	  bbase_y = &bbase[y*bystride];
+	  s = ('rtype_name`) 0;
+	  for (n = 0; n < count; n++)
+	    s += abase[n*axstride] * bbase_y[n*bxstride];
+	  dest[y*rxstride] = s;
+	}
+    }
+  else
+    {
+      const 'rtype_name` *restrict abase_x;
+      const 'rtype_name` *restrict bbase_y;
+      'rtype_name` *restrict dest_y;
+      'rtype_name` s;
+
+      for (y = 0; y < ycount; y++)
+	{
+	  bbase_y = &bbase[y*bystride];
+	  dest_y = &dest[y*rystride];
+	  for (x = 0; x < xcount; x++)
+	    {
+	      abase_x = &abase[x*axstride];
+	      s = ('rtype_name`) 0;
+	      for (n = 0; n < count; n++)
+		s += abase_x[n*aystride] * bbase_y[n*bxstride];
+	      dest_y[x*rxstride] = s;
+	    }
+	}
+    }
+}
+
+#endif'
diff --git a/libgfortran/m4/matmull.m4 b/libgfortran/m4/matmull.m4
new file mode 100644
index 000000000..c5bad25f7
--- /dev/null
+++ b/libgfortran/m4/matmull.m4
@@ -0,0 +1,241 @@
+`/* Implementation of the MATMUL intrinsic
+   Copyright 2002, 2005, 2006, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+/* Dimensions: retarray(x,y) a(x, count) b(count,y).
+   Either a or b can be rank 1.  In this case x or y is 1.  */
+
+extern void matmul_'rtype_code` ('rtype` * const restrict, 
+	gfc_array_l1 * const restrict, gfc_array_l1 * const restrict);
+export_proto(matmul_'rtype_code`);
+
+void
+matmul_'rtype_code` ('rtype` * const restrict retarray, 
+	gfc_array_l1 * const restrict a, gfc_array_l1 * const restrict b)
+{
+  const GFC_LOGICAL_1 * restrict abase;
+  const GFC_LOGICAL_1 * restrict bbase;
+  'rtype_name` * restrict dest;
+  index_type rxstride;
+  index_type rystride;
+  index_type xcount;
+  index_type ycount;
+  index_type xstride;
+  index_type ystride;
+  index_type x;
+  index_type y;
+  int a_kind;
+  int b_kind;
+
+  const GFC_LOGICAL_1 * restrict pa;
+  const GFC_LOGICAL_1 * restrict pb;
+  index_type astride;
+  index_type bstride;
+  index_type count;
+  index_type n;
+
+  assert (GFC_DESCRIPTOR_RANK (a) == 2
+          || GFC_DESCRIPTOR_RANK (b) == 2);
+
+  if (retarray->data == NULL)
+    {
+      if (GFC_DESCRIPTOR_RANK (a) == 1)
+        {
+	  GFC_DIMENSION_SET(retarray->dim[0], 0,
+	                    GFC_DESCRIPTOR_EXTENT(b,1) - 1, 1);
+        }
+      else if (GFC_DESCRIPTOR_RANK (b) == 1)
+        {
+	  GFC_DIMENSION_SET(retarray->dim[0], 0,
+	                    GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
+        }
+      else
+        {
+	  GFC_DIMENSION_SET(retarray->dim[0], 0,
+	                    GFC_DESCRIPTOR_EXTENT(a,0) - 1, 1);
+
+          GFC_DIMENSION_SET(retarray->dim[1], 0,
+	                    GFC_DESCRIPTOR_EXTENT(b,1) - 1,
+			    GFC_DESCRIPTOR_EXTENT(retarray,0));
+        }
+          
+      retarray->data
+	= internal_malloc_size (sizeof ('rtype_name`) * size0 ((array_t *) retarray));
+      retarray->offset = 0;
+    }
+    else if (unlikely (compile_options.bounds_check))
+      {
+	index_type ret_extent, arg_extent;
+
+	if (GFC_DESCRIPTOR_RANK (a) == 1)
+	  {
+	    arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
+	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
+	    if (arg_extent != ret_extent)
+	      runtime_error ("Incorrect extent in return array in"
+			     " MATMUL intrinsic: is %ld, should be %ld",
+			     (long int) ret_extent, (long int) arg_extent);
+	  }
+	else if (GFC_DESCRIPTOR_RANK (b) == 1)
+	  {
+	    arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
+	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
+	    if (arg_extent != ret_extent)
+	      runtime_error ("Incorrect extent in return array in"
+			     " MATMUL intrinsic: is %ld, should be %ld",
+			     (long int) ret_extent, (long int) arg_extent);	    
+	  }
+	else
+	  {
+	    arg_extent = GFC_DESCRIPTOR_EXTENT(a,0);
+	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,0);
+	    if (arg_extent != ret_extent)
+	      runtime_error ("Incorrect extent in return array in"
+			     " MATMUL intrinsic for dimension 1:"
+			     " is %ld, should be %ld",
+			     (long int) ret_extent, (long int) arg_extent);
+
+	    arg_extent = GFC_DESCRIPTOR_EXTENT(b,1);
+	    ret_extent = GFC_DESCRIPTOR_EXTENT(retarray,1);
+	    if (arg_extent != ret_extent)
+	      runtime_error ("Incorrect extent in return array in"
+			     " MATMUL intrinsic for dimension 2:"
+			     " is %ld, should be %ld",
+			     (long int) ret_extent, (long int) arg_extent);
+	  }
+      }
+
+  abase = a->data;
+  a_kind = GFC_DESCRIPTOR_SIZE (a);
+
+  if (a_kind == 1 || a_kind == 2 || a_kind == 4 || a_kind == 8
+#ifdef HAVE_GFC_LOGICAL_16
+     || a_kind == 16
+#endif
+     )
+    abase = GFOR_POINTER_TO_L1 (abase, a_kind);
+  else
+    internal_error (NULL, "Funny sized logical array");
+
+  bbase = b->data;
+  b_kind = GFC_DESCRIPTOR_SIZE (b);
+
+  if (b_kind == 1 || b_kind == 2 || b_kind == 4 || b_kind == 8
+#ifdef HAVE_GFC_LOGICAL_16
+     || b_kind == 16
+#endif
+     )
+    bbase = GFOR_POINTER_TO_L1 (bbase, b_kind);
+  else
+    internal_error (NULL, "Funny sized logical array");
+
+  dest = retarray->data;
+'
+sinclude(`matmul_asm_'rtype_code`.m4')dnl
+`
+  if (GFC_DESCRIPTOR_RANK (retarray) == 1)
+    {
+      rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
+      rystride = rxstride;
+    }
+  else
+    {
+      rxstride = GFC_DESCRIPTOR_STRIDE(retarray,0);
+      rystride = GFC_DESCRIPTOR_STRIDE(retarray,1);
+    }
+
+  /* If we have rank 1 parameters, zero the absent stride, and set the size to
+     one.  */
+  if (GFC_DESCRIPTOR_RANK (a) == 1)
+    {
+      astride = GFC_DESCRIPTOR_STRIDE_BYTES(a,0);
+      count = GFC_DESCRIPTOR_EXTENT(a,0);
+      xstride = 0;
+      rxstride = 0;
+      xcount = 1;
+    }
+  else
+    {
+      astride = GFC_DESCRIPTOR_STRIDE_BYTES(a,1);
+      count = GFC_DESCRIPTOR_EXTENT(a,1);
+      xstride = GFC_DESCRIPTOR_STRIDE_BYTES(a,0);
+      xcount = GFC_DESCRIPTOR_EXTENT(a,0);
+    }
+  if (GFC_DESCRIPTOR_RANK (b) == 1)
+    {
+      bstride = GFC_DESCRIPTOR_STRIDE_BYTES(b,0);
+      assert(count == GFC_DESCRIPTOR_EXTENT(b,0));
+      ystride = 0;
+      rystride = 0;
+      ycount = 1;
+    }
+  else
+    {
+      bstride = GFC_DESCRIPTOR_STRIDE_BYTES(b,0);
+      assert(count == GFC_DESCRIPTOR_EXTENT(b,0));
+      ystride = GFC_DESCRIPTOR_STRIDE_BYTES(b,1);
+      ycount = GFC_DESCRIPTOR_EXTENT(b,1);
+    }
+
+  for (y = 0; y < ycount; y++)
+    {
+      for (x = 0; x < xcount; x++)
+        {
+          /* Do the summation for this element.  For real and integer types
+             this is the same as DOT_PRODUCT.  For complex types we use do
+             a*b, not conjg(a)*b.  */
+          pa = abase;
+          pb = bbase;
+          *dest = 0;
+
+          for (n = 0; n < count; n++)
+            {
+              if (*pa && *pb)
+                {
+                  *dest = 1;
+                  break;
+                }
+              pa += astride;
+              pb += bstride;
+            }
+
+          dest += rxstride;
+          abase += xstride;
+        }
+      abase -= xstride * xcount;
+      bbase += ystride;
+      dest += rystride - (rxstride * xcount);
+    }
+}
+
+#endif
+'
diff --git a/libgfortran/m4/maxloc0.m4 b/libgfortran/m4/maxloc0.m4
new file mode 100644
index 000000000..aef7b9b8d
--- /dev/null
+++ b/libgfortran/m4/maxloc0.m4
@@ -0,0 +1,126 @@
+`/* Implementation of the MAXLOC intrinsic
+   Copyright 2002, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <limits.h>'
+
+include(iparm.m4)dnl
+include(iforeach.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+FOREACH_FUNCTION(
+`    atype_name maxval;
+#if defined('atype_nan`)
+    int fast = 0;
+#endif
+
+#if defined('atype_inf`)
+    maxval = -atype_inf;
+#else
+    maxval = atype_min;
+#endif',
+`#if defined('atype_nan`)
+	}
+      while (0);
+      if (unlikely (!fast))
+	{
+	  do
+	    {
+	      if (*base >= maxval)
+		{
+		  fast = 1;
+		  maxval = *base;
+		  for (n = 0; n < rank; n++)
+		    dest[n * dstride] = count[n] + 1;
+		  break;
+		}
+	      base += sstride[0];
+	    }
+	  while (++count[0] != extent[0]);
+	  if (likely (fast))
+	    continue;
+	}
+      else do
+	{
+#endif
+	  if (*base > maxval)
+	    {
+	      maxval = *base;
+	      for (n = 0; n < rank; n++)
+		dest[n * dstride] = count[n] + 1;
+	    }')
+
+MASKED_FOREACH_FUNCTION(
+`  atype_name maxval;
+   int fast = 0;
+
+#if defined('atype_inf`)
+    maxval = -atype_inf;
+#else
+    maxval = atype_min;
+#endif',
+`	}
+      while (0);
+      if (unlikely (!fast))
+	{
+	  do
+	    {
+	      if (*mbase)
+		{
+#if defined('atype_nan`)
+		  if (unlikely (dest[0] == 0))
+		    for (n = 0; n < rank; n++)
+		      dest[n * dstride] = count[n] + 1;
+		  if (*base >= maxval)
+#endif
+		    {
+		      fast = 1;
+		      maxval = *base;
+		      for (n = 0; n < rank; n++)
+			dest[n * dstride] = count[n] + 1;
+		      break;
+		    }
+		}
+	      base += sstride[0];
+	      mbase += mstride[0];
+	    }
+	  while (++count[0] != extent[0]);
+	  if (likely (fast))
+	    continue;
+	}
+      else do
+	{
+	  if (*mbase && *base > maxval)
+	    {
+	      maxval = *base;
+	      for (n = 0; n < rank; n++)
+		dest[n * dstride] = count[n] + 1;
+	    }')
+
+SCALAR_FOREACH_FUNCTION(`0')
+#endif
diff --git a/libgfortran/m4/maxloc1.m4 b/libgfortran/m4/maxloc1.m4
new file mode 100644
index 000000000..44872ee61
--- /dev/null
+++ b/libgfortran/m4/maxloc1.m4
@@ -0,0 +1,101 @@
+`/* Implementation of the MAXLOC intrinsic
+   Copyright 2002, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <limits.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(0,
+`	atype_name maxval;
+#if defined ('atype_inf`)
+	maxval = -atype_inf;
+#else
+	maxval = atype_min;
+#endif
+	result = 1;',
+`#if defined ('atype_nan`)
+		if (*src >= maxval)
+		  {
+		    maxval = *src;
+		    result = (rtype_name)n + 1;
+		    break;
+		  }
+	      }
+	    for (; n < len; n++, src += delta)
+	      {
+#endif
+		if (*src > maxval)
+		  {
+		    maxval = *src;
+		    result = (rtype_name)n + 1;
+		  }', `')
+
+MASKED_ARRAY_FUNCTION(0,
+`	atype_name maxval;
+#if defined ('atype_inf`)
+	maxval = -atype_inf;
+#else
+	maxval = atype_min;
+#endif
+#if defined ('atype_nan`)
+	rtype_name result2 = 0;
+#endif
+	result = 0;',
+`		if (*msrc)
+		  {
+#if defined ('atype_nan`)
+		    if (!result2)
+		      result2 = (rtype_name)n + 1;
+		    if (*src >= maxval)
+#endif
+		      {
+			maxval = *src;
+			result = (rtype_name)n + 1;
+			break;
+		      }
+		  }
+	      }
+#if defined ('atype_nan`)
+	    if (unlikely (n >= len))
+	      result = result2;
+	    else
+#endif
+	    for (; n < len; n++, src += delta, msrc += mdelta)
+	      {
+		if (*msrc && *src > maxval)
+		  {
+		    maxval = *src;
+		    result = (rtype_name)n + 1;
+		  }')
+
+SCALAR_ARRAY_FUNCTION(0)
+
+#endif
diff --git a/libgfortran/m4/maxval.m4 b/libgfortran/m4/maxval.m4
new file mode 100644
index 000000000..153906880
--- /dev/null
+++ b/libgfortran/m4/maxval.m4
@@ -0,0 +1,88 @@
+`/* Implementation of the MAXVAL intrinsic
+   Copyright 2002, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(atype_min,
+`#if defined ('atype_inf`)
+	result = -atype_inf;
+#else
+	result = atype_min;
+#endif',
+`#if defined ('atype_nan`)
+		if (*src >= result)
+		  break;
+	      }
+	    if (unlikely (n >= len))
+	      result = atype_nan;
+	    else for (; n < len; n++, src += delta)
+	      {
+#endif
+		if (*src > result)
+		  result = *src;', `')
+
+MASKED_ARRAY_FUNCTION(atype_min,
+`#if defined ('atype_inf`)
+	result = -atype_inf;
+#else
+	result = atype_min;
+#endif
+#if defined ('atype_nan`)
+	int non_empty_p = 0;
+#endif',
+`#if defined ('atype_inf`) || defined ('atype_nan`)
+		if (*msrc)
+		  {
+#if defined ('atype_nan`)
+		    non_empty_p = 1;
+		    if (*src >= result)
+#endif
+		      break;
+		  }
+	      }
+	    if (unlikely (n >= len))
+	      {
+#if defined ('atype_nan`)
+		result = non_empty_p ? atype_nan : atype_min;
+#else
+		result = atype_min;
+#endif
+	      }
+	    else for (; n < len; n++, src += delta, msrc += mdelta)
+	      {
+#endif
+		if (*msrc && *src > result)
+		  result = *src;')
+
+SCALAR_ARRAY_FUNCTION(atype_min)
+
+#endif
diff --git a/libgfortran/m4/minloc0.m4 b/libgfortran/m4/minloc0.m4
new file mode 100644
index 000000000..1345f2df4
--- /dev/null
+++ b/libgfortran/m4/minloc0.m4
@@ -0,0 +1,126 @@
+`/* Implementation of the MINLOC intrinsic
+   Copyright 2002, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <limits.h>'
+
+include(iparm.m4)dnl
+include(iforeach.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+FOREACH_FUNCTION(
+`    atype_name minval;
+#if defined('atype_nan`)
+    int fast = 0;
+#endif
+
+#if defined('atype_inf`)
+    minval = atype_inf;
+#else
+    minval = atype_max;
+#endif',
+`#if defined('atype_nan`)
+	}
+      while (0);
+      if (unlikely (!fast))
+	{
+	  do
+	    {
+	      if (*base <= minval)
+		{
+		  fast = 1;
+		  minval = *base;
+		  for (n = 0; n < rank; n++)
+		    dest[n * dstride] = count[n] + 1;
+		  break;
+		}
+	      base += sstride[0];
+	    }
+	  while (++count[0] != extent[0]);
+	  if (likely (fast))
+	    continue;
+	}
+      else do
+	{
+#endif
+	  if (*base < minval)
+	    {
+	      minval = *base;
+	      for (n = 0; n < rank; n++)
+		dest[n * dstride] = count[n] + 1;
+	    }')
+
+MASKED_FOREACH_FUNCTION(
+`  atype_name minval;
+   int fast = 0;
+
+#if defined('atype_inf`)
+    minval = atype_inf;
+#else
+    minval = atype_max;
+#endif',
+`	}
+      while (0);
+      if (unlikely (!fast))
+	{
+	  do
+	    {
+	      if (*mbase)
+		{
+#if defined('atype_nan`)
+		  if (unlikely (dest[0] == 0))
+		    for (n = 0; n < rank; n++)
+		      dest[n * dstride] = count[n] + 1;
+		  if (*base <= minval)
+#endif
+		    {
+		      fast = 1;
+		      minval = *base;
+		      for (n = 0; n < rank; n++)
+			dest[n * dstride] = count[n] + 1;
+		      break;
+		    }
+		}
+	      base += sstride[0];
+	      mbase += mstride[0];
+	    }
+	  while (++count[0] != extent[0]);
+	  if (likely (fast))
+	    continue;
+	}
+      else do
+	{
+	  if (*mbase && *base < minval)
+	    {
+	      minval = *base;
+	      for (n = 0; n < rank; n++)
+		dest[n * dstride] = count[n] + 1;
+	    }')
+
+SCALAR_FOREACH_FUNCTION(`0')
+#endif
diff --git a/libgfortran/m4/minloc1.m4 b/libgfortran/m4/minloc1.m4
new file mode 100644
index 000000000..e0bd3cd88
--- /dev/null
+++ b/libgfortran/m4/minloc1.m4
@@ -0,0 +1,101 @@
+`/* Implementation of the MINLOC intrinsic
+   Copyright 2002, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <limits.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(0,
+`	atype_name minval;
+#if defined ('atype_inf`)
+	minval = atype_inf;
+#else
+	minval = atype_max;
+#endif
+	result = 1;',
+`#if defined ('atype_nan`)
+		if (*src <= minval)
+		  {
+		    minval = *src;
+		    result = (rtype_name)n + 1;
+		    break;
+		  }
+	      }
+	    for (; n < len; n++, src += delta)
+	      {
+#endif
+		if (*src < minval)
+		  {
+		    minval = *src;
+		    result = (rtype_name)n + 1;
+		  }')
+
+MASKED_ARRAY_FUNCTION(0,
+`	atype_name minval;
+#if defined ('atype_inf`)
+	minval = atype_inf;
+#else
+	minval = atype_max;
+#endif
+#if defined ('atype_nan`)
+	rtype_name result2 = 0;
+#endif
+	result = 0;',
+`		if (*msrc)
+		  {
+#if defined ('atype_nan`)
+		    if (!result2)
+		      result2 = (rtype_name)n + 1;
+		    if (*src <= minval)
+#endif
+		      {
+			minval = *src;
+			result = (rtype_name)n + 1;
+			break;
+		      }
+		  }
+	      }
+#if defined ('atype_nan`)
+	    if (unlikely (n >= len))
+	      result = result2;
+	    else
+#endif
+	    for (; n < len; n++, src += delta, msrc += mdelta)
+	      {
+		if (*msrc && *src < minval)
+		  {
+		    minval = *src;
+		    result = (rtype_name)n + 1;
+		  }', `')
+
+SCALAR_ARRAY_FUNCTION(0)
+
+#endif
diff --git a/libgfortran/m4/minval.m4 b/libgfortran/m4/minval.m4
new file mode 100644
index 000000000..037dd6725
--- /dev/null
+++ b/libgfortran/m4/minval.m4
@@ -0,0 +1,88 @@
+`/* Implementation of the MINVAL intrinsic
+   Copyright 2002, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(atype_max,
+`#if defined ('atype_inf`)
+	result = atype_inf;
+#else
+	result = atype_max;
+#endif',
+`#if defined ('atype_nan`)
+		if (*src <= result)
+		  break;
+	      }
+	    if (unlikely (n >= len))
+	      result = atype_nan;
+	    else for (; n < len; n++, src += delta)
+	      {
+#endif
+		if (*src < result)
+		  result = *src;')
+
+MASKED_ARRAY_FUNCTION(atype_max,
+`#if defined ('atype_inf`)
+	result = atype_inf;
+#else
+	result = atype_max;
+#endif
+#if defined ('atype_nan`)
+	int non_empty_p = 0;
+#endif',
+`#if defined ('atype_inf`) || defined ('atype_nan`)
+		if (*msrc)
+		  {
+#if defined ('atype_nan`)
+		    non_empty_p = 1;
+		    if (*src <= result)
+#endif
+		      break;
+		  }
+	      }
+	    if (unlikely (n >= len))
+	      {
+#if defined ('atype_nan`)
+		result = non_empty_p ? atype_nan : atype_max;
+#else
+		result = atype_max;
+#endif
+	      }
+	    else for (; n < len; n++, src += delta, msrc += mdelta)
+	      {
+#endif
+		if (*msrc && *src < result)
+		  result = *src;', `')
+
+SCALAR_ARRAY_FUNCTION(atype_max)
+
+#endif
diff --git a/libgfortran/m4/misc_specifics.m4 b/libgfortran/m4/misc_specifics.m4
new file mode 100644
index 000000000..3e40bf017
--- /dev/null
+++ b/libgfortran/m4/misc_specifics.m4
@@ -0,0 +1,64 @@
+include(head.m4)dnl
+dnl
+dnl This file contains the specific functions that are not handled in the
+dnl m4/specific.m4 file.
+
+#include "config.h"
+#include "kinds.inc"
+
+dnl This is from GNU m4 examples file foreach.m4:
+divert(-1)
+# foreach(x, (item_1, item_2, ..., item_n), stmt)
+define(`foreach', `pushdef(`$1', `')_foreach(`$1', `$2',
+`$3')popdef(`$1')')
+define(`_arg1', `$1')
+define(`_foreach',
+        `ifelse(`$2', `()', ,
+                `define(`$1', _arg1$2)$3`'_foreach(`$1', (shift$2),
+`$3')')')
+# traceon(`define', `foreach', `_foreach', `ifelse')
+divert
+
+dnl   NINT specifics
+foreach(`ikind', `(4, 8, 16)', `foreach(`rkind', `(4, 8, 10, 16)', `
+`#if defined (HAVE_GFC_REAL_'rkind`) && defined (HAVE_GFC_INTEGER_'ikind`)'
+elemental function _gfortran_specific__nint_`'ikind`_'rkind (parm)
+   real (kind=rkind) , intent (in) :: parm
+   integer (kind=ikind) :: _gfortran_specific__nint_`'ikind`_'rkind
+   _gfortran_specific__nint_`'ikind`_'rkind = nint (parm)
+end function
+#endif
+')')
+
+dnl   CHAR specifics
+foreach(`ckind', `(1)', `foreach(`ikind', `(4, 8, 16)', `
+`#if defined (HAVE_GFC_INTEGER_'ikind`)'
+elemental function _gfortran_specific__char_`'ckind`_i'ikind (parm)
+   integer (kind=ikind) , intent (in) :: parm
+   character (kind=ckind,len=1) :: _gfortran_specific__char_`'ckind`_i'ikind
+   _gfortran_specific__char_`'ckind`_i'ikind` = char (parm, kind='ckind`)'
+end function
+#endif
+')')
+
+dnl   LEN specifics
+foreach(`ckind', `(1)', `foreach(`ikind', `(4, 8, 16)', `
+`#if defined (HAVE_GFC_INTEGER_'ikind`)'
+elemental function _gfortran_specific__len_`'ckind`_i'ikind (parm)
+   character (kind=ckind,len=*) , intent (in) :: parm
+   integer (kind=ikind) :: _gfortran_specific__len_`'ckind`_i'ikind
+   _gfortran_specific__len_`'ckind`_i'ikind` = len (parm)'
+end function
+#endif
+')')
+
+dnl   INDEX specifics
+foreach(`ckind', `(1)', `foreach(`ikind', `(4, 8, 16)', `
+`#if defined (HAVE_GFC_INTEGER_'ikind`)'
+elemental function _gfortran_specific__index_`'ckind`_i'ikind (parm1, parm2)
+   character (kind=ckind,len=*) , intent (in) :: parm1, parm2
+   integer (kind=ikind) :: _gfortran_specific__index_`'ckind`_i'ikind
+   _gfortran_specific__index_`'ckind`_i'ikind` = index (parm1, parm2)'
+end function
+#endif
+')')
diff --git a/libgfortran/m4/mtype.m4 b/libgfortran/m4/mtype.m4
new file mode 100644
index 000000000..b133e578b
--- /dev/null
+++ b/libgfortran/m4/mtype.m4
@@ -0,0 +1,13 @@
+dnl Get type kind from filename.
+define(kind,regexp(file, `_.\([0-9]+\).c$', `\1'))dnl
+define(complex_type, `GFC_COMPLEX_'kind)dnl
+define(real_type, `GFC_REAL_'kind)dnl
+define(`upcase', `translit(`$*', `a-z', `A-Z')')dnl
+define(q,ifelse(kind,4,f,ifelse(kind,8,`',ifelse(kind,10,l,ifelse(kind,16,l,`_'kind)))))dnl
+define(Q,translit(q,`a-z',`A-Z'))dnl
+define(hasmathfunc,`ifelse(kind,4,`defined (HAVE_'upcase($1)`F)',ifelse(kind,8,`defined (HAVE_'upcase($1)`)',ifelse(kind,10,`defined (HAVE_'upcase($1)`L)',ifelse(kind,16,`(defined(GFC_REAL_16_IS_FLOAT128) || defined(HAVE_'upcase($1)`L))',`error out'))))')
+define(mathfunc_macro,`ifelse(kind,16,`#if defined(GFC_REAL_16_IS_FLOAT128)
+#define MATHFUNC(funcname) funcname ## q
+#else
+#define MATHFUNC(funcname) funcname ## l
+#endif',ifelse(kind,8,``#''`define MATHFUNC(funcname) funcname',```#'''`define MATHFUNC(funcname) funcname '```#'''```#'''` 'q))')dnl
diff --git a/libgfortran/m4/nearest.m4 b/libgfortran/m4/nearest.m4
new file mode 100644
index 000000000..643443048
--- /dev/null
+++ b/libgfortran/m4/nearest.m4
@@ -0,0 +1,51 @@
+`/* Implementation of the NEAREST intrinsic
+   Copyright 2003, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Richard Henderson <rth@redhat.com>.
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+mathfunc_macro
+
+`#if defined (HAVE_'real_type`) && 'hasmathfunc(copysign) && hasmathfunc(nextafter)`
+
+extern 'real_type` nearest_r'kind` ('real_type` s, 'real_type` dir);
+export_proto(nearest_r'kind`);
+
+'real_type`
+nearest_r'kind` ('real_type` s, 'real_type` dir)
+{
+  dir = MATHFUNC(copysign) (MATHFUNC(__builtin_inf) (), dir);
+  if (FLT_EVAL_METHOD != 0)
+    {
+      /* ??? Work around glibc bug on x86.  */
+      volatile 'real_type` r = MATHFUNC(nextafter) (s, dir);
+      return r;
+    }
+  else
+    return MATHFUNC(nextafter) (s, dir);
+}
+
+#endif'
diff --git a/libgfortran/m4/norm2.m4 b/libgfortran/m4/norm2.m4
new file mode 100644
index 000000000..b2162ac23
--- /dev/null
+++ b/libgfortran/m4/norm2.m4
@@ -0,0 +1,61 @@
+`/* Implementation of the NORM2 intrinsic
+   Copyright 2010 Free Software Foundation, Inc.
+   Contributed by Tobias Burnus  <burnus@net-b.de>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <math.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+include(`mtype.m4')dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`) && 'hasmathfunc(sqrt) && hasmathfunc(fabs)
+
+mathfunc_macro
+
+ARRAY_FUNCTION(0,
+`	'rtype_name` scale;
+	result = 0;
+	scale = 1;',
+`	  if (*src != 0)
+	    {
+	      'rtype_name` absX, val;
+	      absX = MATHFUNC(fabs) (*src);
+	      if (scale < absX)
+		{
+		  val = scale / absX;
+		  result = 1 + result * val * val;
+		  scale = absX;
+		}
+	      else
+		{
+		  val = absX / scale;
+		  result += val * val;
+		}
+	    }',
+`   result = scale * MATHFUNC(sqrt) (result);')
+
+#endif
diff --git a/libgfortran/m4/pack.m4 b/libgfortran/m4/pack.m4
new file mode 100644
index 000000000..c5fd2fd81
--- /dev/null
+++ b/libgfortran/m4/pack.m4
@@ -0,0 +1,262 @@
+`/* Specific implementation of the PACK intrinsic
+   Copyright (C) 2002, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Ligbfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+/* PACK is specified as follows:
+
+   13.14.80 PACK (ARRAY, MASK, [VECTOR])
+
+   Description: Pack an array into an array of rank one under the
+   control of a mask.
+
+   Class: Transformational function.
+
+   Arguments:
+      ARRAY   may be of any type. It shall not be scalar.
+      MASK    shall be of type LOGICAL. It shall be conformable with ARRAY.
+      VECTOR  (optional) shall be of the same type and type parameters
+              as ARRAY. VECTOR shall have at least as many elements as
+              there are true elements in MASK. If MASK is a scalar
+              with the value true, VECTOR shall have at least as many
+              elements as there are in ARRAY.
+
+   Result Characteristics: The result is an array of rank one with the
+   same type and type parameters as ARRAY. If VECTOR is present, the
+   result size is that of VECTOR; otherwise, the result size is the
+   number /t/ of true elements in MASK unless MASK is scalar with the
+   value true, in which case the result size is the size of ARRAY.
+
+   Result Value: Element /i/ of the result is the element of ARRAY
+   that corresponds to the /i/th true element of MASK, taking elements
+   in array element order, for /i/ = 1, 2, ..., /t/. If VECTOR is
+   present and has size /n/ > /t/, element /i/ of the result has the
+   value VECTOR(/i/), for /i/ = /t/ + 1, ..., /n/.
+
+   Examples: The nonzero elements of an array M with the value
+   | 0 0 0 |
+   | 9 0 0 | may be "gathered" by the function PACK. The result of
+   | 0 0 7 |
+   PACK (M, MASK = M.NE.0) is [9,7] and the result of PACK (M, M.NE.0,
+   VECTOR = (/ 2,4,6,8,10,12 /)) is [9,7,6,8,10,12].
+
+There are two variants of the PACK intrinsic: one, where MASK is
+array valued, and the other one where MASK is scalar.  */
+
+void
+pack_'rtype_code` ('rtype` *ret, const 'rtype` *array,
+	       const gfc_array_l1 *mask, const 'rtype` *vector)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride0;
+  'rtype_name` * restrict rptr;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  const 'rtype_name` *sptr;
+  /* m.* indicates the mask array.  */
+  index_type mstride[GFC_MAX_DIMENSIONS];
+  index_type mstride0;
+  const GFC_LOGICAL_1 *mptr;
+
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  int zero_sized;
+  index_type n;
+  index_type dim;
+  index_type nelem;
+  index_type total;
+  int mask_kind;
+
+  dim = GFC_DESCRIPTOR_RANK (array);
+
+  mptr = mask->data;
+
+  /* Use the same loop for all logical types, by using GFC_LOGICAL_1
+     and using shifting to address size and endian issues.  */
+
+  mask_kind = GFC_DESCRIPTOR_SIZE (mask);
+
+  if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
+#ifdef HAVE_GFC_LOGICAL_16
+      || mask_kind == 16
+#endif
+      )
+    {
+      /*  Do not convert a NULL pointer as we use test for NULL below.  */
+      if (mptr)
+	mptr = GFOR_POINTER_TO_L1 (mptr, mask_kind);
+    }
+  else
+    runtime_error ("Funny sized logical array");
+
+  zero_sized = 0;
+  for (n = 0; n < dim; n++)
+    {
+      count[n] = 0;
+      extent[n] = GFC_DESCRIPTOR_EXTENT(array,n);
+      if (extent[n] <= 0)
+       zero_sized = 1;
+      sstride[n] = GFC_DESCRIPTOR_STRIDE(array,n);
+      mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
+    }
+  if (sstride[0] == 0)
+    sstride[0] = 1;
+  if (mstride[0] == 0)
+    mstride[0] = mask_kind;
+
+  if (zero_sized)
+    sptr = NULL;
+  else
+    sptr = array->data;
+
+  if (ret->data == NULL || unlikely (compile_options.bounds_check))
+    {
+      /* Count the elements, either for allocating memory or
+	 for bounds checking.  */
+
+      if (vector != NULL)
+	{
+	  /* The return array will have as many
+	     elements as there are in VECTOR.  */
+	  total = GFC_DESCRIPTOR_EXTENT(vector,0);
+	  if (total < 0)
+	    {
+	      total = 0;
+	      vector = NULL;
+	    }
+	}
+      else
+        {
+      	  /* We have to count the true elements in MASK.  */
+	  total = count_0 (mask);
+        }
+
+      if (ret->data == NULL)
+	{
+	  /* Setup the array descriptor.  */
+	  GFC_DIMENSION_SET(ret->dim[0], 0, total-1, 1);
+
+	  ret->offset = 0;
+	  if (total == 0)
+	    {
+	      /* In this case, nothing remains to be done.  */
+	      ret->data = internal_malloc_size (1);
+	      return;
+	    }
+	  else
+	    ret->data = internal_malloc_size (sizeof ('rtype_name`) * total);
+	}
+      else 
+	{
+	  /* We come here because of range checking.  */
+	  index_type ret_extent;
+
+	  ret_extent = GFC_DESCRIPTOR_EXTENT(ret,0);
+	  if (total != ret_extent)
+	    runtime_error ("Incorrect extent in return value of PACK intrinsic;"
+			   " is %ld, should be %ld", (long int) total,
+			   (long int) ret_extent);
+	}
+    }
+
+  rstride0 = GFC_DESCRIPTOR_STRIDE(ret,0);
+  if (rstride0 == 0)
+    rstride0 = 1;
+  sstride0 = sstride[0];
+  mstride0 = mstride[0];
+  rptr = ret->data;
+
+  while (sptr && mptr)
+    {
+      /* Test this element.  */
+      if (*mptr)
+        {
+          /* Add it.  */
+	  *rptr = *sptr;
+          rptr += rstride0;
+        }
+      /* Advance to the next element.  */
+      sptr += sstride0;
+      mptr += mstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          sptr -= sstride[n] * extent[n];
+          mptr -= mstride[n] * extent[n];
+          n++;
+          if (n >= dim)
+            {
+              /* Break out of the loop.  */
+              sptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              sptr += sstride[n];
+              mptr += mstride[n];
+            }
+        }
+    }
+
+  /* Add any remaining elements from VECTOR.  */
+  if (vector)
+    {
+      n = GFC_DESCRIPTOR_EXTENT(vector,0);
+      nelem = ((rptr - ret->data) / rstride0);
+      if (n > nelem)
+        {
+          sstride0 = GFC_DESCRIPTOR_STRIDE(vector,0);
+          if (sstride0 == 0)
+            sstride0 = 1;
+
+          sptr = vector->data + sstride0 * nelem;
+          n -= nelem;
+          while (n--)
+            {
+	      *rptr = *sptr;
+              rptr += rstride0;
+              sptr += sstride0;
+            }
+        }
+    }
+}
+
+#endif
+'
diff --git a/libgfortran/m4/parity.m4 b/libgfortran/m4/parity.m4
new file mode 100644
index 000000000..037e96db9
--- /dev/null
+++ b/libgfortran/m4/parity.m4
@@ -0,0 +1,40 @@
+`/* Implementation of the NORM2 intrinsic
+   Copyright 2010 Free Software Foundation, Inc.
+   Contributed by Tobias Burnus  <burnus@net-b.de>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <math.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(0,
+`  result = 0;',
+`  result = result != *src;', `')
+
+#endif
diff --git a/libgfortran/m4/pow.m4 b/libgfortran/m4/pow.m4
new file mode 100644
index 000000000..3814f87fe
--- /dev/null
+++ b/libgfortran/m4/pow.m4
@@ -0,0 +1,83 @@
+`/* Support routines for the intrinsic power (**) operator.
+   Copyright 2004, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+
+/* Use Binary Method to calculate the powi. This is not an optimal but
+   a simple and reasonable arithmetic. See section 4.6.3, "Evaluation of
+   Powers" of Donald E. Knuth, "Seminumerical Algorithms", Vol. 2, "The Art
+   of Computer Programming", 3rd Edition, 1998.  */
+
+`#if defined (HAVE_'rtype_name`) && defined (HAVE_'atype_name`)'
+
+rtype_name `pow_'rtype_code`_'atype_code` ('rtype_name` a, 'atype_name` b);
+export_proto(pow_'rtype_code`_'atype_code`);
+
+'rtype_name`
+pow_'rtype_code`_'atype_code` ('rtype_name` a, 'atype_name` b)
+{
+  'rtype_name` pow, x;
+  'atype_name` n;
+  GFC_UINTEGER_'atype_kind` u;
+  
+  n = b;
+  x = a;
+  pow = 1;
+  if (n != 0)
+    {
+      if (n < 0)
+	{
+'ifelse(rtype_letter,i,`dnl
+	  if (x == 1)
+	    return 1;
+	  if (x == -1)
+	    return (n & 1) ? -1 : 1;
+	  return (x == 0) ? 1 / x : 0;
+',`
+	  u = -n;
+	  x = pow / x;
+')dnl
+`	}
+      else
+	{
+	   u = n;
+	}
+      for (;;)
+	{
+	  if (u & 1)
+	    pow *= x;
+	  u >>= 1;
+	  if (u)
+	    x *= x;
+	  else
+	    break;
+	}
+    }
+  return pow;
+}
+
+#endif'
diff --git a/libgfortran/m4/product.m4 b/libgfortran/m4/product.m4
new file mode 100644
index 000000000..6f6f2c011
--- /dev/null
+++ b/libgfortran/m4/product.m4
@@ -0,0 +1,46 @@
+`/* Implementation of the PRODUCT intrinsic
+   Copyright 2002, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(1,
+`  result = 1;',
+`  result *= *src;', `')
+
+MASKED_ARRAY_FUNCTION(1,
+`  result = 1;',
+`  if (*msrc)
+    result *= *src;')
+
+SCALAR_ARRAY_FUNCTION(1)
+
+`#endif'
diff --git a/libgfortran/m4/reshape.m4 b/libgfortran/m4/reshape.m4
new file mode 100644
index 000000000..4052a5ecc
--- /dev/null
+++ b/libgfortran/m4/reshape.m4
@@ -0,0 +1,356 @@
+`/* Implementation of the RESHAPE intrinsic
+   Copyright 2002, 2006, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+typedef GFC_ARRAY_DESCRIPTOR(1, 'index_type`) 'shape_type`;'
+
+dnl For integer routines, only the kind (ie size) is used to name the
+dnl function.  The same function will be used for integer and logical
+dnl arrays of the same kind.
+
+`extern void reshape_'rtype_ccode` ('rtype` * const restrict, 
+	'rtype` * const restrict, 
+	'shape_type` * const restrict,
+	'rtype` * const restrict, 
+	'shape_type` * const restrict);
+export_proto(reshape_'rtype_ccode`);
+
+void
+reshape_'rtype_ccode` ('rtype` * const restrict ret, 
+	'rtype` * const restrict source, 
+	'shape_type` * const restrict shape,
+	'rtype` * const restrict pad, 
+	'shape_type` * const restrict order)
+{
+  /* r.* indicates the return array.  */
+  index_type rcount[GFC_MAX_DIMENSIONS];
+  index_type rextent[GFC_MAX_DIMENSIONS];
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type rdim;
+  index_type rsize;
+  index_type rs;
+  index_type rex;
+  'rtype_name` *rptr;
+  /* s.* indicates the source array.  */
+  index_type scount[GFC_MAX_DIMENSIONS];
+  index_type sextent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type sdim;
+  index_type ssize;
+  const 'rtype_name` *sptr;
+  /* p.* indicates the pad array.  */
+  index_type pcount[GFC_MAX_DIMENSIONS];
+  index_type pextent[GFC_MAX_DIMENSIONS];
+  index_type pstride[GFC_MAX_DIMENSIONS];
+  index_type pdim;
+  index_type psize;
+  const 'rtype_name` *pptr;
+
+  const 'rtype_name` *src;
+  int n;
+  int dim;
+  int sempty, pempty, shape_empty;
+  index_type shape_data[GFC_MAX_DIMENSIONS];
+
+  rdim = GFC_DESCRIPTOR_EXTENT(shape,0);
+  if (rdim != GFC_DESCRIPTOR_RANK(ret))
+    runtime_error("rank of return array incorrect in RESHAPE intrinsic");
+
+  shape_empty = 0;
+
+  for (n = 0; n < rdim; n++)
+    {
+      shape_data[n] = shape->data[n * GFC_DESCRIPTOR_STRIDE(shape,0)];
+      if (shape_data[n] <= 0)
+      {
+        shape_data[n] = 0;
+	shape_empty = 1;
+      }
+    }
+
+  if (ret->data == NULL)
+    {
+      rs = 1;
+      for (n = 0; n < rdim; n++)
+	{
+	  rex = shape_data[n];
+
+	  GFC_DIMENSION_SET(ret->dim[n], 0, rex - 1, rs);
+
+	  rs *= rex;
+	}
+      ret->offset = 0;
+      ret->data = internal_malloc_size ( rs * sizeof ('rtype_name`));
+      ret->dtype = (source->dtype & ~GFC_DTYPE_RANK_MASK) | rdim;
+    }
+
+  if (shape_empty)
+    return;
+
+  if (pad)
+    {
+      pdim = GFC_DESCRIPTOR_RANK (pad);
+      psize = 1;
+      pempty = 0;
+      for (n = 0; n < pdim; n++)
+        {
+          pcount[n] = 0;
+          pstride[n] = GFC_DESCRIPTOR_STRIDE(pad,n);
+          pextent[n] = GFC_DESCRIPTOR_EXTENT(pad,n);
+          if (pextent[n] <= 0)
+	    {
+	      pempty = 1;
+	      pextent[n] = 0;
+	    }
+
+          if (psize == pstride[n])
+            psize *= pextent[n];
+          else
+            psize = 0;
+        }
+      pptr = pad->data;
+    }
+  else
+    {
+      pdim = 0;
+      psize = 1;
+      pempty = 1;
+      pptr = NULL;
+    }
+
+  if (unlikely (compile_options.bounds_check))
+    {
+      index_type ret_extent, source_extent;
+
+      rs = 1;
+      for (n = 0; n < rdim; n++)
+	{
+	  rs *= shape_data[n];
+	  ret_extent = GFC_DESCRIPTOR_EXTENT(ret,n);
+	  if (ret_extent != shape_data[n])
+	    runtime_error("Incorrect extent in return value of RESHAPE"
+			  " intrinsic in dimension %ld: is %ld,"
+			  " should be %ld", (long int) n+1,
+			  (long int) ret_extent, (long int) shape_data[n]);
+	}
+
+      source_extent = 1;
+      sdim = GFC_DESCRIPTOR_RANK (source);
+      for (n = 0; n < sdim; n++)
+	{
+	  index_type se;
+	  se = GFC_DESCRIPTOR_EXTENT(source,n);
+	  source_extent *= se > 0 ? se : 0;
+	}
+
+      if (rs > source_extent && (!pad || pempty))
+	runtime_error("Incorrect size in SOURCE argument to RESHAPE"
+		      " intrinsic: is %ld, should be %ld",
+		      (long int) source_extent, (long int) rs);
+
+      if (order)
+	{
+	  int seen[GFC_MAX_DIMENSIONS];
+	  index_type v;
+
+	  for (n = 0; n < rdim; n++)
+	    seen[n] = 0;
+
+	  for (n = 0; n < rdim; n++)
+	    {
+	      v = order->data[n * GFC_DESCRIPTOR_STRIDE(order,0)] - 1;
+
+	      if (v < 0 || v >= rdim)
+		runtime_error("Value %ld out of range in ORDER argument"
+			      " to RESHAPE intrinsic", (long int) v + 1);
+
+	      if (seen[v] != 0)
+		runtime_error("Duplicate value %ld in ORDER argument to"
+			      " RESHAPE intrinsic", (long int) v + 1);
+		
+	      seen[v] = 1;
+	    }
+	}
+    }
+
+  rsize = 1;
+  for (n = 0; n < rdim; n++)
+    {
+      if (order)
+        dim = order->data[n * GFC_DESCRIPTOR_STRIDE(order,0)] - 1;
+      else
+        dim = n;
+
+      rcount[n] = 0;
+      rstride[n] = GFC_DESCRIPTOR_STRIDE(ret,dim);
+      rextent[n] = GFC_DESCRIPTOR_EXTENT(ret,dim);
+      if (rextent[n] < 0)
+        rextent[n] = 0;
+
+      if (rextent[n] != shape_data[dim])
+        runtime_error ("shape and target do not conform");
+
+      if (rsize == rstride[n])
+        rsize *= rextent[n];
+      else
+        rsize = 0;
+      if (rextent[n] <= 0)
+        return;
+    }
+
+  sdim = GFC_DESCRIPTOR_RANK (source);
+  ssize = 1;
+  sempty = 0;
+  for (n = 0; n < sdim; n++)
+    {
+      scount[n] = 0;
+      sstride[n] = GFC_DESCRIPTOR_STRIDE(source,n);
+      sextent[n] = GFC_DESCRIPTOR_EXTENT(source,n);
+      if (sextent[n] <= 0)
+	{
+	  sempty = 1;
+	  sextent[n] = 0;
+	}
+
+      if (ssize == sstride[n])
+        ssize *= sextent[n];
+      else
+        ssize = 0;
+    }
+
+  if (rsize != 0 && ssize != 0 && psize != 0)
+    {
+      rsize *= sizeof ('rtype_name`);
+      ssize *= sizeof ('rtype_name`);
+      psize *= sizeof ('rtype_name`);
+      reshape_packed ((char *)ret->data, rsize, (char *)source->data,
+		      ssize, pad ? (char *)pad->data : NULL, psize);
+      return;
+    }
+  rptr = ret->data;
+  src = sptr = source->data;
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+
+  if (sempty && pempty)
+    abort ();
+
+  if (sempty)
+    {
+      /* Pretend we are using the pad array the first time around, too.  */
+      src = pptr;
+      sptr = pptr;
+      sdim = pdim;
+      for (dim = 0; dim < pdim; dim++)
+	{
+	  scount[dim] = pcount[dim];
+	  sextent[dim] = pextent[dim];
+	  sstride[dim] = pstride[dim];
+	  sstride0 = pstride[0];
+	}
+    }
+
+  while (rptr)
+    {
+      /* Select between the source and pad arrays.  */
+      *rptr = *src;
+      /* Advance to the next element.  */
+      rptr += rstride0;
+      src += sstride0;
+      rcount[0]++;
+      scount[0]++;
+
+      /* Advance to the next destination element.  */
+      n = 0;
+      while (rcount[n] == rextent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          rcount[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          rptr -= rstride[n] * rextent[n];
+          n++;
+          if (n == rdim)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              rcount[n]++;
+              rptr += rstride[n];
+            }
+        }
+      /* Advance to the next source element.  */
+      n = 0;
+      while (scount[n] == sextent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          scount[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          src -= sstride[n] * sextent[n];
+          n++;
+          if (n == sdim)
+            {
+              if (sptr && pad)
+                {
+                  /* Switch to the pad array.  */
+                  sptr = NULL;
+                  sdim = pdim;
+                  for (dim = 0; dim < pdim; dim++)
+                    {
+                      scount[dim] = pcount[dim];
+                      sextent[dim] = pextent[dim];
+                      sstride[dim] = pstride[dim];
+                      sstride0 = sstride[0];
+                    }
+                }
+              /* We now start again from the beginning of the pad array.  */
+              src = pptr;
+              break;
+            }
+          else
+            {
+              scount[n]++;
+              src += sstride[n];
+            }
+        }
+    }
+}
+
+#endif'
diff --git a/libgfortran/m4/rrspacing.m4 b/libgfortran/m4/rrspacing.m4
new file mode 100644
index 000000000..5f11d5b6b
--- /dev/null
+++ b/libgfortran/m4/rrspacing.m4
@@ -0,0 +1,54 @@
+`/* Implementation of the RRSPACING intrinsic
+   Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Steven G. Kargl <kargl@gcc.gnu.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+mathfunc_macro
+
+`#if defined (HAVE_'real_type`) && 'hasmathfunc(fabs) && hasmathfunc(frexp)`
+
+extern 'real_type` rrspacing_r'kind` ('real_type` s, int p);
+export_proto(rrspacing_r'kind`);
+
+'real_type`
+rrspacing_r'kind` ('real_type` s, int p)
+{
+  int e;
+  'real_type` x;
+  x = MATHFUNC(fabs) (s);
+  if (x == 0.)
+    return 0.;
+  MATHFUNC(frexp) (s, &e);
+#if 'hasmathfunc(ldexp)`
+  return MATHFUNC(ldexp) (x, p - e);
+#else
+  return MATHFUNC(scalbn) (x, p - e);
+#endif
+
+}
+
+#endif'
diff --git a/libgfortran/m4/set_exponent.m4 b/libgfortran/m4/set_exponent.m4
new file mode 100644
index 000000000..5f31ab681
--- /dev/null
+++ b/libgfortran/m4/set_exponent.m4
@@ -0,0 +1,44 @@
+`/* Implementation of the SET_EXPONENT intrinsic
+   Copyright 2003, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Richard Henderson <rth@redhat.com>.
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+mathfunc_macro
+
+`#if defined (HAVE_'real_type`) && 'hasmathfunc(scalbn) && hasmathfunc(frexp)`
+
+extern 'real_type` set_exponent_r'kind` ('real_type` s, GFC_INTEGER_4 i);
+export_proto(set_exponent_r'kind`);
+
+'real_type`
+set_exponent_r'kind` ('real_type` s, GFC_INTEGER_4 i)
+{
+  int dummy_exp;
+  return MATHFUNC(scalbn) (MATHFUNC(frexp) (s, &dummy_exp), i);
+}
+
+#endif'
diff --git a/libgfortran/m4/shape.m4 b/libgfortran/m4/shape.m4
new file mode 100644
index 000000000..c079cca99
--- /dev/null
+++ b/libgfortran/m4/shape.m4
@@ -0,0 +1,68 @@
+`/* Implementation of the SHAPE intrinsic
+   Copyright 2002, 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+extern void shape_'rtype_kind` ('rtype` * const restrict ret, 
+	const 'rtype` * const restrict array);
+export_proto(shape_'rtype_kind`);
+
+void
+shape_'rtype_kind` ('rtype` * const restrict ret, 
+	const 'rtype` * const restrict array)
+{
+  int n;
+  index_type stride;
+  index_type extent;
+  int rank;
+
+  rank = GFC_DESCRIPTOR_RANK (array);
+
+  if (ret->data == NULL)
+    {
+      GFC_DIMENSION_SET(ret->dim[0], 0, rank - 1, 1);
+      ret->offset = 0;
+      ret->data = internal_malloc_size (sizeof ('rtype_name`) * rank);
+    }
+
+  stride = GFC_DESCRIPTOR_STRIDE(ret,0);
+
+  if (GFC_DESCRIPTOR_EXTENT(ret,0) < 1)
+    return;
+
+  for (n = 0; n < rank; n++)
+    {
+      extent = GFC_DESCRIPTOR_EXTENT(array,n);
+      ret->data[n * stride] = extent > 0 ? extent : 0 ;
+    }
+}
+
+#endif'
diff --git a/libgfortran/m4/spacing.m4 b/libgfortran/m4/spacing.m4
new file mode 100644
index 000000000..faae6d124
--- /dev/null
+++ b/libgfortran/m4/spacing.m4
@@ -0,0 +1,53 @@
+`/* Implementation of the SPACING intrinsic
+   Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Contributed by Steven G. Kargl <kargl@gcc.gnu.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+mathfunc_macro
+
+`#if defined (HAVE_'real_type`) && 'hasmathfunc(frexp)`
+
+extern 'real_type` spacing_r'kind` ('real_type` s, int p, int emin, 'real_type` tiny);
+export_proto(spacing_r'kind`);
+
+'real_type`
+spacing_r'kind` ('real_type` s, int p, int emin, 'real_type` tiny)
+{
+  int e;
+  if (s == 0.)
+    return tiny;
+  MATHFUNC(frexp) (s, &e);
+  e = e - p;
+  e = e > emin ? e : emin;
+#if 'hasmathfunc(ldexp)`
+  return MATHFUNC(ldexp) (1., e);
+#else
+  return MATHFUNC(scalbn) (1., e);
+#endif
+}
+
+#endif'
diff --git a/libgfortran/m4/specific.m4 b/libgfortran/m4/specific.m4
new file mode 100644
index 000000000..ebc89839a
--- /dev/null
+++ b/libgfortran/m4/specific.m4
@@ -0,0 +1,43 @@
+include(head.m4)
+define(atype_code,regexp(file,`_\([ircl][0-9]+\).[fF]90',`\1'))dnl
+define(atype_letter,substr(atype_code, 0, 1))dnl
+define(atype_kind,substr(atype_code, 1))dnl
+define(get_typename2, `$1 (kind=$2)')dnl
+define(get_typename, `get_typename2(ifelse($1,i,integer,ifelse($1,r,real,ifelse($1,l,logical,ifelse($1,c,complex,unknown)))),`$2')')dnl
+define(atype_name, get_typename(atype_letter,atype_kind))dnl
+define(name, regexp(regexp(file, `[^/]*$', `\&'), `^_\([^_]*\)_', `\1'))dnl
+define(rtype_name,get_typename(ifelse(name,abs,ifelse(atype_letter,c,r,atype_letter),ifelse(name,aimag,ifelse(atype_letter,c,r,atype_letter),atype_letter)),atype_kind))dnl
+define(function_name,ifelse(name,conjg,`_gfortran_specific__conjg_'atype_kind,`_gfortran_specific__'name`_'atype_code))dnl
+
+define(type,ifelse(atype_letter,l,LOGICAL,ifelse(atype_letter,i,INTEGER,ifelse(atype_letter,r,REAL,ifelse(atype_letter,c,COMPLEX,UNKNOW)))))dnl
+define(Q,ifelse(atype_kind,4,F,ifelse(atype_kind,8,`',ifelse(atype_kind,10,L,ifelse(atype_kind,16,L,`_'atype_kind)))))dnl
+
+dnl A few specifics require a function other than their name, or
+dnl nothing. The list is currently:
+dnl    - integer and logical specifics require no libm function
+dnl    - AINT requires the trunc() family functions
+dnl    - ANINT requires round()
+dnl    - AIMAG, CONJG, DIM, SIGN require no libm function
+define(needed,ifelse(atype_letter,i,`none',ifelse(atype_letter,l,`none',ifelse(name,aint,trunc,ifelse(name,anint,round,ifelse(name,aimag,none,ifelse(name,conjg,none,ifelse(name,dim,none,ifelse(name,sign,none,ifelse(name,abs,fabs,name))))))))))dnl
+define(prefix,ifelse(atype_letter,c,C,`'))dnl
+
+dnl Special case for fabs, for which the corresponding complex function
+dnl is not cfabs but cabs.
+define(NEEDED,translit(ifelse(prefix`'needed,`Cfabs',`abs',needed),`a-z',`A-Z'))dnl
+
+#include "config.h"
+#include "kinds.inc"
+#include "c99_protos.inc"
+
+`#if defined (HAVE_GFC_'type`_'atype_kind`)'
+ifelse(NEEDED,NONE,`',`#ifdef HAVE_'prefix`'NEEDED`'Q)
+
+elemental function function_name (parm)
+   atype_name, intent (in) :: parm
+   rtype_name :: function_name
+
+   function_name = name (parm)
+end function
+
+ifelse(NEEDED,NONE,`',`#endif')
+#endif
diff --git a/libgfortran/m4/specific2.m4 b/libgfortran/m4/specific2.m4
new file mode 100644
index 000000000..d05e8db14
--- /dev/null
+++ b/libgfortran/m4/specific2.m4
@@ -0,0 +1,30 @@
+include(head.m4)
+define(atype_code,regexp(file,`_\([ircl][0-9]+\).[fF]90',`\1'))dnl
+define(atype_letter,substr(atype_code, 0, 1))dnl
+define(atype_kind,substr(atype_code, 1))dnl
+define(get_typename2, `$1 (kind=$2)')dnl
+define(get_typename, `get_typename2(ifelse($1,i,integer,ifelse($1,r,real,ifelse($1,l,logical,ifelse($1,c,complex,unknown)))),`$2')')dnl
+define(atype_name, get_typename(atype_letter,atype_kind))dnl
+define(name, regexp(regexp(file, `[^/]*$', `\&'), `^_\([^_]*\)_', `\1'))dnl
+define(function_name,`_gfortran_specific__'name`_'atype_code)dnl
+
+define(Q,ifelse(atype_kind,4,F,ifelse(atype_kind,8,`',ifelse(atype_kind,10,L,ifelse(atype_kind,16,L,`_'atype_kind)))))dnl
+
+#include "config.h"
+#include "kinds.inc"
+#include "c99_protos.inc"
+
+`#if defined (HAVE_GFC_'ifelse(atype_letter,l,LOGICAL,ifelse(atype_letter,i,INTEGER,ifelse(atype_letter,r,REAL,ifelse(atype_letter,c,COMPLEX,UNKNOW))))`_'atype_kind`)'
+
+ifelse(name,atan2,`#ifdef HAVE_ATAN2'Q,)
+
+elemental function function_name (p1, p2)
+   atype_name, intent (in) :: p1, p2
+   atype_name :: function_name
+
+   function_name = name (p1, p2)
+end function
+
+ifelse(name,atan2,`#endif',)
+
+#endif
diff --git a/libgfortran/m4/spread.m4 b/libgfortran/m4/spread.m4
new file mode 100644
index 000000000..5e73d9742
--- /dev/null
+++ b/libgfortran/m4/spread.m4
@@ -0,0 +1,272 @@
+`/* Special implementation of the SPREAD intrinsic
+   Copyright 2008, 2009 Free Software Foundation, Inc.
+   Contributed by Thomas Koenig <tkoenig@gcc.gnu.org>, based on
+   spread_generic.c written by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Ligbfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+void
+spread_'rtype_code` ('rtype` *ret, const 'rtype` *source,
+		 const index_type along, const index_type pncopies)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type rdelta = 0;
+  index_type rrank;
+  index_type rs;
+  'rtype_name` *rptr;
+  'rtype_name` * restrict dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const 'rtype_name` *sptr;
+
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type n;
+  index_type dim;
+  index_type ncopies;
+
+  srank = GFC_DESCRIPTOR_RANK(source);
+
+  rrank = srank + 1;
+  if (rrank > GFC_MAX_DIMENSIONS)
+    runtime_error ("return rank too large in spread()");
+
+  if (along > rrank)
+      runtime_error ("dim outside of rank in spread()");
+
+  ncopies = pncopies;
+
+  if (ret->data == NULL)
+    {
+
+      size_t ub, stride;
+
+      /* The front end has signalled that we need to populate the
+	 return array descriptor.  */
+      ret->dtype = (source->dtype & ~GFC_DTYPE_RANK_MASK) | rrank;
+      dim = 0;
+      rs = 1;
+      for (n = 0; n < rrank; n++)
+	{
+	  stride = rs;
+	  if (n == along - 1)
+	    {
+	      ub = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = GFC_DESCRIPTOR_EXTENT(source,dim);
+	      sstride[dim] = GFC_DESCRIPTOR_STRIDE(source,dim);
+	      rstride[dim] = rs;
+
+	      ub = extent[dim] - 1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	  GFC_DIMENSION_SET(ret->dim[n], 0, ub, stride);
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof('rtype_name`));
+      else
+	{
+	  ret->data = internal_malloc_size (1);
+	  return;
+	}
+    }
+  else
+    {
+      int zero_sized;
+
+      zero_sized = 0;
+
+      dim = 0;
+      if (GFC_DESCRIPTOR_RANK(ret) != rrank)
+	runtime_error ("rank mismatch in spread()");
+
+      if (unlikely (compile_options.bounds_check))
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = GFC_DESCRIPTOR_EXTENT(ret,n);
+	      if (n == along - 1)
+		{
+		  rdelta = GFC_DESCRIPTOR_STRIDE(ret,n);
+
+		  if (ret_extent != ncopies)
+		    runtime_error("Incorrect extent in return value of SPREAD"
+				  " intrinsic in dimension %ld: is %ld,"
+				  " should be %ld", (long int) n+1,
+				  (long int) ret_extent, (long int) ncopies);
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = GFC_DESCRIPTOR_EXTENT(source,dim);
+		  if (ret_extent != extent[dim])
+		    runtime_error("Incorrect extent in return value of SPREAD"
+				  " intrinsic in dimension %ld: is %ld,"
+				  " should be %ld", (long int) n+1,
+				  (long int) ret_extent,
+				  (long int) extent[dim]);
+		    
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = GFC_DESCRIPTOR_STRIDE(source,dim);
+		  rstride[dim] = GFC_DESCRIPTOR_STRIDE(ret,n);
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = GFC_DESCRIPTOR_STRIDE(ret,n);
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = GFC_DESCRIPTOR_EXTENT(source,dim);
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = GFC_DESCRIPTOR_STRIDE(source,dim);
+		  rstride[dim] = GFC_DESCRIPTOR_STRIDE(ret,n);
+		  dim++;
+		}
+	    }
+	}
+
+      if (zero_sized)
+	return;
+
+      if (sstride[0] == 0)
+	sstride[0] = 1;
+    }
+  sstride0 = sstride[0];
+  rstride0 = rstride[0];
+  rptr = ret->data;
+  sptr = source->data;
+
+  while (sptr)
+    {
+      /* Spread this element.  */
+      dest = rptr;
+      for (n = 0; n < ncopies; n++)
+        {
+	  *dest = *sptr;
+          dest += rdelta;
+        }
+      /* Advance to the next element.  */
+      sptr += sstride0;
+      rptr += rstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          sptr -= sstride[n] * extent[n];
+          rptr -= rstride[n] * extent[n];
+          n++;
+          if (n >= srank)
+            {
+              /* Break out of the loop.  */
+              sptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              sptr += sstride[n];
+              rptr += rstride[n];
+            }
+        }
+    }
+}
+
+/* This version of spread_internal treats the special case of a scalar
+   source.  This is much simpler than the more general case above.  */
+
+void
+spread_scalar_'rtype_code` ('rtype` *ret, const 'rtype_name` *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  'rtype_name` * restrict dest;
+  index_type stride;
+
+  if (GFC_DESCRIPTOR_RANK (ret) != 1)
+    runtime_error ("incorrect destination rank in spread()");
+
+  if (along > 1)
+    runtime_error ("dim outside of rank in spread()");
+
+  if (ret->data == NULL)
+    {
+      ret->data = internal_malloc_size (ncopies * sizeof ('rtype_name`));
+      ret->offset = 0;
+      GFC_DIMENSION_SET(ret->dim[0], 0, ncopies - 1, 1);
+    }
+  else
+    {
+      if (ncopies - 1 > (GFC_DESCRIPTOR_EXTENT(ret,0) - 1)
+			   / GFC_DESCRIPTOR_STRIDE(ret,0))
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = GFC_DESCRIPTOR_STRIDE(ret,0);
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
+'
diff --git a/libgfortran/m4/sum.m4 b/libgfortran/m4/sum.m4
new file mode 100644
index 000000000..b502c6e7d
--- /dev/null
+++ b/libgfortran/m4/sum.m4
@@ -0,0 +1,46 @@
+`/* Implementation of the SUM intrinsic
+   Copyright 2002, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+
+`#if defined (HAVE_'atype_name`) && defined (HAVE_'rtype_name`)'
+
+ARRAY_FUNCTION(0,
+`  result = 0;',
+`  result += *src;')
+
+MASKED_ARRAY_FUNCTION(0,
+`  result = 0;',
+`  if (*msrc)
+    result += *src;')
+
+SCALAR_ARRAY_FUNCTION(0)
+
+#endif
diff --git a/libgfortran/m4/transpose.m4 b/libgfortran/m4/transpose.m4
new file mode 100644
index 000000000..34c2d6c06
--- /dev/null
+++ b/libgfortran/m4/transpose.m4
@@ -0,0 +1,115 @@
+`/* Implementation of the TRANSPOSE intrinsic
+   Copyright 2003, 2005, 2006, 2007, 2009 Free Software Foundation, Inc.
+   Contributed by Tobias Schlüter
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Libgfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <assert.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+extern void transpose_'rtype_code` ('rtype` * const restrict ret, 
+	'rtype` * const restrict source);
+export_proto(transpose_'rtype_code`);
+
+void
+transpose_'rtype_code` ('rtype` * const restrict ret, 
+	'rtype` * const restrict source)
+{
+  /* r.* indicates the return array.  */
+  index_type rxstride, rystride;
+  'rtype_name` * restrict rptr;
+  /* s.* indicates the source array.  */
+  index_type sxstride, systride;
+  const 'rtype_name` *sptr;
+
+  index_type xcount, ycount;
+  index_type x, y;
+
+  assert (GFC_DESCRIPTOR_RANK (source) == 2);
+
+  if (ret->data == NULL)
+    {
+      assert (GFC_DESCRIPTOR_RANK (ret) == 2);
+      assert (ret->dtype == source->dtype);
+
+      GFC_DIMENSION_SET(ret->dim[0], 0, GFC_DESCRIPTOR_EXTENT(source,1) - 1,
+			1);
+
+      GFC_DIMENSION_SET(ret->dim[1], 0, GFC_DESCRIPTOR_EXTENT(source,0) - 1,
+			GFC_DESCRIPTOR_EXTENT(source, 1));
+
+      ret->data = internal_malloc_size (sizeof ('rtype_name`) * size0 ((array_t *) ret));
+      ret->offset = 0;
+    } else if (unlikely (compile_options.bounds_check))
+    {
+      index_type ret_extent, src_extent;
+
+      ret_extent = GFC_DESCRIPTOR_EXTENT(ret,0);
+      src_extent = GFC_DESCRIPTOR_EXTENT(source,1);
+
+      if (src_extent != ret_extent)
+	runtime_error ("Incorrect extent in return value of TRANSPOSE"
+		       " intrinsic in dimension 1: is %ld,"
+		       " should be %ld", (long int) src_extent,
+		       (long int) ret_extent);
+
+      ret_extent = GFC_DESCRIPTOR_EXTENT(ret,1);
+      src_extent = GFC_DESCRIPTOR_EXTENT(source,0);
+
+      if (src_extent != ret_extent)
+	runtime_error ("Incorrect extent in return value of TRANSPOSE"
+		       " intrinsic in dimension 2: is %ld,"
+		       " should be %ld", (long int) src_extent,
+		       (long int) ret_extent);
+
+    }
+
+  sxstride = GFC_DESCRIPTOR_STRIDE(source,0);
+  systride = GFC_DESCRIPTOR_STRIDE(source,1);
+  xcount = GFC_DESCRIPTOR_EXTENT(source,0);
+  ycount = GFC_DESCRIPTOR_EXTENT(source,1);
+
+  rxstride = GFC_DESCRIPTOR_STRIDE(ret,0);
+  rystride = GFC_DESCRIPTOR_STRIDE(ret,1);
+
+  rptr = ret->data;
+  sptr = source->data;
+
+  for (y=0; y < ycount; y++)
+    {
+      for (x=0; x < xcount; x++)
+        {
+          *rptr = *sptr;
+
+          sptr += sxstride;
+          rptr += rystride;
+        }
+        sptr += systride - (sxstride * xcount);
+        rptr += rxstride - (rystride * xcount);
+    }
+}
+
+#endif'
diff --git a/libgfortran/m4/types.m4 b/libgfortran/m4/types.m4
new file mode 100644
index 000000000..cb808290c
--- /dev/null
+++ b/libgfortran/m4/types.m4
@@ -0,0 +1,4 @@
+define(get_typename2, `GFC_$1_$2')dnl
+define(get_typename, `get_typename2(ifelse($1,i,INTEGER,ifelse($1,r,REAL,ifelse($1,l,LOGICAL,ifelse($1,c,COMPLEX,unknown)))),`$2')')dnl
+define(get_arraytype, `gfc_array_$1$2')dnl
+define(name, regexp(regexp(file, `[^/]*$', `\&'), `^\([^_]*\)_', `\1'))dnl
diff --git a/libgfortran/m4/unpack.m4 b/libgfortran/m4/unpack.m4
new file mode 100644
index 000000000..bf348aebe
--- /dev/null
+++ b/libgfortran/m4/unpack.m4
@@ -0,0 +1,332 @@
+`/* Specific implementation of the UNPACK intrinsic
+   Copyright 2008, 2009 Free Software Foundation, Inc.
+   Contributed by Thomas Koenig <tkoenig@gcc.gnu.org>, based on
+   unpack_generic.c by Paul Brook <paul@nowt.org>.
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran 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 of the License, or (at your option) any later version.
+
+Ligbfortran 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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>'
+
+include(iparm.m4)dnl
+
+`#if defined (HAVE_'rtype_name`)
+
+void
+unpack0_'rtype_code` ('rtype` *ret, const 'rtype` *vector,
+		 const gfc_array_l1 *mask, const 'rtype_name` *fptr)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type rs;
+  'rtype_name` * restrict rptr;
+  /* v.* indicates the vector array.  */
+  index_type vstride0;
+  'rtype_name` *vptr;
+  /* Value for field, this is constant.  */
+  const 'rtype_name` fval = *fptr;
+  /* m.* indicates the mask array.  */
+  index_type mstride[GFC_MAX_DIMENSIONS];
+  index_type mstride0;
+  const GFC_LOGICAL_1 *mptr;
+
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type n;
+  index_type dim;
+
+  int empty;
+  int mask_kind;
+
+  empty = 0;
+
+  mptr = mask->data;
+
+  /* Use the same loop for all logical types, by using GFC_LOGICAL_1
+     and using shifting to address size and endian issues.  */
+
+  mask_kind = GFC_DESCRIPTOR_SIZE (mask);
+
+  if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
+#ifdef HAVE_GFC_LOGICAL_16
+      || mask_kind == 16
+#endif
+      )
+    {
+      /*  Do not convert a NULL pointer as we use test for NULL below.  */
+      if (mptr)
+	mptr = GFOR_POINTER_TO_L1 (mptr, mask_kind);
+    }
+  else
+    runtime_error ("Funny sized logical array");
+
+  if (ret->data == NULL)
+    {
+      /* The front end has signalled that we need to populate the
+	 return array descriptor.  */
+      dim = GFC_DESCRIPTOR_RANK (mask);
+      rs = 1;
+      for (n = 0; n < dim; n++)
+	{
+	  count[n] = 0;
+	  GFC_DIMENSION_SET(ret->dim[n], 0,
+			    GFC_DESCRIPTOR_EXTENT(mask,n) - 1, rs);
+	  extent[n] = GFC_DESCRIPTOR_EXTENT(ret,n);
+	  empty = empty || extent[n] <= 0;
+	  rstride[n] = GFC_DESCRIPTOR_STRIDE(ret,n);
+	  mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
+	  rs *= extent[n];
+	}
+      ret->offset = 0;
+      ret->data = internal_malloc_size (rs * sizeof ('rtype_name`));
+    }
+  else
+    {
+      dim = GFC_DESCRIPTOR_RANK (ret);
+      for (n = 0; n < dim; n++)
+	{
+	  count[n] = 0;
+	  extent[n] = GFC_DESCRIPTOR_EXTENT(ret,n);
+	  empty = empty || extent[n] <= 0;
+	  rstride[n] = GFC_DESCRIPTOR_STRIDE(ret,n);
+	  mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
+	}
+      if (rstride[0] == 0)
+	rstride[0] = 1;
+    }
+
+  if (empty)
+    return;
+
+  if (mstride[0] == 0)
+    mstride[0] = 1;
+
+  vstride0 = GFC_DESCRIPTOR_STRIDE(vector,0);
+  if (vstride0 == 0)
+    vstride0 = 1;
+  rstride0 = rstride[0];
+  mstride0 = mstride[0];
+  rptr = ret->data;
+  vptr = vector->data;
+
+  while (rptr)
+    {
+      if (*mptr)
+        {
+	  /* From vector.  */
+	  *rptr = *vptr;
+	  vptr += vstride0;
+        }
+      else
+        {
+	  /* From field.  */
+	  *rptr = fval;
+        }
+      /* Advance to the next element.  */
+      rptr += rstride0;
+      mptr += mstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          rptr -= rstride[n] * extent[n];
+          mptr -= mstride[n] * extent[n];
+          n++;
+          if (n >= dim)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              rptr += rstride[n];
+              mptr += mstride[n];
+            }
+        }
+    }
+}
+
+void
+unpack1_'rtype_code` ('rtype` *ret, const 'rtype` *vector,
+		 const gfc_array_l1 *mask, const 'rtype` *field)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type rs;
+  'rtype_name` * restrict rptr;
+  /* v.* indicates the vector array.  */
+  index_type vstride0;
+  'rtype_name` *vptr;
+  /* f.* indicates the field array.  */
+  index_type fstride[GFC_MAX_DIMENSIONS];
+  index_type fstride0;
+  const 'rtype_name` *fptr;
+  /* m.* indicates the mask array.  */
+  index_type mstride[GFC_MAX_DIMENSIONS];
+  index_type mstride0;
+  const GFC_LOGICAL_1 *mptr;
+
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type n;
+  index_type dim;
+
+  int empty;
+  int mask_kind;
+
+  empty = 0;
+
+  mptr = mask->data;
+
+  /* Use the same loop for all logical types, by using GFC_LOGICAL_1
+     and using shifting to address size and endian issues.  */
+
+  mask_kind = GFC_DESCRIPTOR_SIZE (mask);
+
+  if (mask_kind == 1 || mask_kind == 2 || mask_kind == 4 || mask_kind == 8
+#ifdef HAVE_GFC_LOGICAL_16
+      || mask_kind == 16
+#endif
+      )
+    {
+      /*  Do not convert a NULL pointer as we use test for NULL below.  */
+      if (mptr)
+	mptr = GFOR_POINTER_TO_L1 (mptr, mask_kind);
+    }
+  else
+    runtime_error ("Funny sized logical array");
+
+  if (ret->data == NULL)
+    {
+      /* The front end has signalled that we need to populate the
+	 return array descriptor.  */
+      dim = GFC_DESCRIPTOR_RANK (mask);
+      rs = 1;
+      for (n = 0; n < dim; n++)
+	{
+	  count[n] = 0;
+	  GFC_DIMENSION_SET(ret->dim[n], 0,
+			    GFC_DESCRIPTOR_EXTENT(mask,n) - 1, rs);
+	  extent[n] = GFC_DESCRIPTOR_EXTENT(ret,n);
+	  empty = empty || extent[n] <= 0;
+	  rstride[n] = GFC_DESCRIPTOR_STRIDE(ret,n);
+	  fstride[n] = GFC_DESCRIPTOR_STRIDE(field,n);
+	  mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
+	  rs *= extent[n];
+	}
+      ret->offset = 0;
+      ret->data = internal_malloc_size (rs * sizeof ('rtype_name`));
+    }
+  else
+    {
+      dim = GFC_DESCRIPTOR_RANK (ret);
+      for (n = 0; n < dim; n++)
+	{
+	  count[n] = 0;
+	  extent[n] = GFC_DESCRIPTOR_EXTENT(ret,n);
+	  empty = empty || extent[n] <= 0;
+	  rstride[n] = GFC_DESCRIPTOR_STRIDE(ret,n);
+	  fstride[n] = GFC_DESCRIPTOR_STRIDE(field,n);
+	  mstride[n] = GFC_DESCRIPTOR_STRIDE_BYTES(mask,n);
+	}
+      if (rstride[0] == 0)
+	rstride[0] = 1;
+    }
+
+  if (empty)
+    return;
+
+  if (fstride[0] == 0)
+    fstride[0] = 1;
+  if (mstride[0] == 0)
+    mstride[0] = 1;
+
+  vstride0 = GFC_DESCRIPTOR_STRIDE(vector,0);
+  if (vstride0 == 0)
+    vstride0 = 1;
+  rstride0 = rstride[0];
+  fstride0 = fstride[0];
+  mstride0 = mstride[0];
+  rptr = ret->data;
+  fptr = field->data;
+  vptr = vector->data;
+
+  while (rptr)
+    {
+      if (*mptr)
+        {
+          /* From vector.  */
+	  *rptr = *vptr;
+          vptr += vstride0;
+        }
+      else
+        {
+          /* From field.  */
+	  *rptr = *fptr;
+        }
+      /* Advance to the next element.  */
+      rptr += rstride0;
+      fptr += fstride0;
+      mptr += mstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so probably not worth it.  */
+          rptr -= rstride[n] * extent[n];
+          fptr -= fstride[n] * extent[n];
+          mptr -= mstride[n] * extent[n];
+          n++;
+          if (n >= dim)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              rptr += rstride[n];
+              fptr += fstride[n];
+              mptr += mstride[n];
+            }
+        }
+    }
+}
+
+#endif
+'