2008-03-28 Thomas Koenig <tkoenig@gcc.gnu.org>

	PR libfortran/32972
	PR libfortran/32512
	* Makefile.am:  Add new variable, i_spread_c, containing
	pack_i1.c, pack_i2.c, pack_i4.c, pack_i8.c, spread_i16.c,
	spread_r4.c, spread_r8.c, spread_r10.c, spread_r16.c,
	spread_c4.c, spread_c8.c, spread_c10.c, spread_c16.c.
	* Makefile.in:  Regenerated.
	* libgfortran.h:  Add prototypes for spread_i1, spread_i2,
	spread_i4, spread_i8, spread_i16, spread_r4, spread_r8,
	spread_c4, spread_c8, spread_c10, spread_c16,
	spread_scalar_i1, spread_scalar_i2, spread_scalar_i4,
	spread_scalar_i8, spread_scalar_i16, spread_scalar_r4
	spread_scalar_r8, spread_scalar_c4, spread_scalar_c8,
	spread_scalar_c10 and spread_scalar_c16.
	Add macros to isolate both type and size information
	from array descriptors with a single mask operation.
	* intrinsics/spread_generic.c:  Add calls to specific
	spread functions.
	* m4/spread.m4:  New file.
	* generated/spread_i1.c:  New file.
	* generated/spread_i2.c:  New file.
	* generated/spread_i4.c:  New file.
	* generated/spread_i8.c:  New file.
	* generated/spread_i16.c:  New file.
	* generated/spread_r4.c:  New file.
	* generated/spread_r8.c:  New file.
	* generated/spread_r10.c:  New file.
	* generated/spread_r16.c:  New file.
	* generated/spread_c4.c:  New file.
	* generated/spread_c8.c:  New file.
	* generated/spread_c10.c:  New file.
	* generated/spread_c16.c:  New file.

2008-03-28  Thomas Koenig  <tkoenig@gcc.gnu.org>

	PR libfortran/32972
	PR libfortran/32512
	* intrinsic_spread_1.f90:  New file.
	* intrinsic_spread_2.f90:  New file.
	* intrinsic_spread_3.f90:  New file.



git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@133702 138bc75d-0d04-0410-961f-82ee72b054a4

13 files changed, 3601 insertions, 0 deletions

diff --git a/libgfortran/generated/spread_c10.c b/libgfortran/generated/spread_c10.c
new file mode 100644
index 00000000000..76a361406c1
--- /dev/null
+++ b/libgfortran/generated/spread_c10.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_COMPLEX_10)
+
+void
+spread_c10 (gfc_array_c10 *ret, const gfc_array_c10 *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;
+  GFC_COMPLEX_10 *rptr;
+  GFC_COMPLEX_10 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_COMPLEX_10 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_COMPLEX_10));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_c10 (gfc_array_c10 *ret, const GFC_COMPLEX_10 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_COMPLEX_10 * 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 (GFC_COMPLEX_10));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_c16.c b/libgfortran/generated/spread_c16.c
new file mode 100644
index 00000000000..0ea57561849
--- /dev/null
+++ b/libgfortran/generated/spread_c16.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_COMPLEX_16)
+
+void
+spread_c16 (gfc_array_c16 *ret, const gfc_array_c16 *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;
+  GFC_COMPLEX_16 *rptr;
+  GFC_COMPLEX_16 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_COMPLEX_16 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_COMPLEX_16));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_c16 (gfc_array_c16 *ret, const GFC_COMPLEX_16 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_COMPLEX_16 * 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 (GFC_COMPLEX_16));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_c4.c b/libgfortran/generated/spread_c4.c
new file mode 100644
index 00000000000..f86da84a58f
--- /dev/null
+++ b/libgfortran/generated/spread_c4.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_COMPLEX_4)
+
+void
+spread_c4 (gfc_array_c4 *ret, const gfc_array_c4 *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;
+  GFC_COMPLEX_4 *rptr;
+  GFC_COMPLEX_4 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_COMPLEX_4 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_COMPLEX_4));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_c4 (gfc_array_c4 *ret, const GFC_COMPLEX_4 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_COMPLEX_4 * 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 (GFC_COMPLEX_4));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_c8.c b/libgfortran/generated/spread_c8.c
new file mode 100644
index 00000000000..7a3f4dfd210
--- /dev/null
+++ b/libgfortran/generated/spread_c8.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_COMPLEX_8)
+
+void
+spread_c8 (gfc_array_c8 *ret, const gfc_array_c8 *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;
+  GFC_COMPLEX_8 *rptr;
+  GFC_COMPLEX_8 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_COMPLEX_8 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_COMPLEX_8));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_c8 (gfc_array_c8 *ret, const GFC_COMPLEX_8 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_COMPLEX_8 * 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 (GFC_COMPLEX_8));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_i1.c b/libgfortran/generated/spread_i1.c
new file mode 100644
index 00000000000..396a521eab8
--- /dev/null
+++ b/libgfortran/generated/spread_i1.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_INTEGER_1)
+
+void
+spread_i1 (gfc_array_i1 *ret, const gfc_array_i1 *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;
+  GFC_INTEGER_1 *rptr;
+  GFC_INTEGER_1 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_INTEGER_1 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_INTEGER_1));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_i1 (gfc_array_i1 *ret, const GFC_INTEGER_1 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_INTEGER_1 * 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 (GFC_INTEGER_1));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_i16.c b/libgfortran/generated/spread_i16.c
new file mode 100644
index 00000000000..55993424054
--- /dev/null
+++ b/libgfortran/generated/spread_i16.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_INTEGER_16)
+
+void
+spread_i16 (gfc_array_i16 *ret, const gfc_array_i16 *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;
+  GFC_INTEGER_16 *rptr;
+  GFC_INTEGER_16 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_INTEGER_16 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_INTEGER_16));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_i16 (gfc_array_i16 *ret, const GFC_INTEGER_16 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_INTEGER_16 * 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 (GFC_INTEGER_16));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_i2.c b/libgfortran/generated/spread_i2.c
new file mode 100644
index 00000000000..d8ac9dc9af1
--- /dev/null
+++ b/libgfortran/generated/spread_i2.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_INTEGER_2)
+
+void
+spread_i2 (gfc_array_i2 *ret, const gfc_array_i2 *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;
+  GFC_INTEGER_2 *rptr;
+  GFC_INTEGER_2 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_INTEGER_2 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_INTEGER_2));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_i2 (gfc_array_i2 *ret, const GFC_INTEGER_2 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_INTEGER_2 * 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 (GFC_INTEGER_2));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_i4.c b/libgfortran/generated/spread_i4.c
new file mode 100644
index 00000000000..c0890b666a1
--- /dev/null
+++ b/libgfortran/generated/spread_i4.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_INTEGER_4)
+
+void
+spread_i4 (gfc_array_i4 *ret, const gfc_array_i4 *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;
+  GFC_INTEGER_4 *rptr;
+  GFC_INTEGER_4 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_INTEGER_4 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_INTEGER_4));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_i4 (gfc_array_i4 *ret, const GFC_INTEGER_4 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_INTEGER_4 * 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 (GFC_INTEGER_4));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_i8.c b/libgfortran/generated/spread_i8.c
new file mode 100644
index 00000000000..b0032bf64dd
--- /dev/null
+++ b/libgfortran/generated/spread_i8.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_INTEGER_8)
+
+void
+spread_i8 (gfc_array_i8 *ret, const gfc_array_i8 *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;
+  GFC_INTEGER_8 *rptr;
+  GFC_INTEGER_8 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_INTEGER_8 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_INTEGER_8));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_i8 (gfc_array_i8 *ret, const GFC_INTEGER_8 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_INTEGER_8 * 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 (GFC_INTEGER_8));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_r10.c b/libgfortran/generated/spread_r10.c
new file mode 100644
index 00000000000..404aaa4654c
--- /dev/null
+++ b/libgfortran/generated/spread_r10.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_REAL_10)
+
+void
+spread_r10 (gfc_array_r10 *ret, const gfc_array_r10 *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;
+  GFC_REAL_10 *rptr;
+  GFC_REAL_10 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_REAL_10 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_REAL_10));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_r10 (gfc_array_r10 *ret, const GFC_REAL_10 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_REAL_10 * 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 (GFC_REAL_10));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_r16.c b/libgfortran/generated/spread_r16.c
new file mode 100644
index 00000000000..122673305e8
--- /dev/null
+++ b/libgfortran/generated/spread_r16.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_REAL_16)
+
+void
+spread_r16 (gfc_array_r16 *ret, const gfc_array_r16 *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;
+  GFC_REAL_16 *rptr;
+  GFC_REAL_16 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_REAL_16 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_REAL_16));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_r16 (gfc_array_r16 *ret, const GFC_REAL_16 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_REAL_16 * 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 (GFC_REAL_16));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_r4.c b/libgfortran/generated/spread_r4.c
new file mode 100644
index 00000000000..1569dbc09f3
--- /dev/null
+++ b/libgfortran/generated/spread_r4.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_REAL_4)
+
+void
+spread_r4 (gfc_array_r4 *ret, const gfc_array_r4 *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;
+  GFC_REAL_4 *rptr;
+  GFC_REAL_4 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_REAL_4 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_REAL_4));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_r4 (gfc_array_r4 *ret, const GFC_REAL_4 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_REAL_4 * 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 (GFC_REAL_4));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/spread_r8.c b/libgfortran/generated/spread_r8.c
new file mode 100644
index 00000000000..c028f804079
--- /dev/null
+++ b/libgfortran/generated/spread_r8.c
@@ -0,0 +1,277 @@
+/* Special implementation of the SPREAD intrinsic
+   Copyright 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include "libgfortran.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+
+
+#if defined (HAVE_GFC_REAL_8)
+
+void
+spread_r8 (gfc_array_r8 *ret, const gfc_array_r8 *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;
+  GFC_REAL_8 *rptr;
+  GFC_REAL_8 *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type srank;
+  const GFC_REAL_8 *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)
+    {
+      /* 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++)
+	{
+	  ret->dim[n].stride = rs;
+	  ret->dim[n].lbound = 0;
+	  if (n == along - 1)
+	    {
+	      ret->dim[n].ubound = ncopies - 1;
+	      rdelta = rs;
+	      rs *= ncopies;
+	    }
+	  else
+	    {
+	      count[dim] = 0;
+	      extent[dim] = source->dim[dim].ubound + 1
+		- source->dim[dim].lbound;
+	      sstride[dim] = source->dim[dim].stride;
+	      rstride[dim] = rs;
+
+	      ret->dim[n].ubound = extent[dim]-1;
+	      rs *= extent[dim];
+	      dim++;
+	    }
+	}
+      ret->offset = 0;
+      if (rs > 0)
+        ret->data = internal_malloc_size (rs * sizeof(GFC_REAL_8));
+      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 (compile_options.bounds_check)
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      index_type ret_extent;
+
+	      ret_extent = ret->dim[n].ubound + 1 - ret->dim[n].lbound;
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+
+		  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] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  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] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  dim++;
+		}
+	    }
+	}
+      else
+	{
+	  for (n = 0; n < rrank; n++)
+	    {
+	      if (n == along - 1)
+		{
+		  rdelta = ret->dim[n].stride;
+		}
+	      else
+		{
+		  count[dim] = 0;
+		  extent[dim] = source->dim[dim].ubound + 1
+		    - source->dim[dim].lbound;
+		  if (extent[dim] <= 0)
+		    zero_sized = 1;
+		  sstride[dim] = source->dim[dim].stride;
+		  rstride[dim] = ret->dim[n].stride;
+		  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_r8 (gfc_array_r8 *ret, const GFC_REAL_8 *source,
+			const index_type along, const index_type pncopies)
+{
+  int n;
+  int ncopies = pncopies;
+  GFC_REAL_8 * 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 (GFC_REAL_8));
+      ret->offset = 0;
+      ret->dim[0].stride = 1;
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = ncopies - 1;
+    }
+  else
+    {
+      if (ncopies - 1 > (ret->dim[0].ubound - ret->dim[0].lbound)
+			   / ret->dim[0].stride)
+	runtime_error ("dim too large in spread()");
+    }
+
+  dest = ret->data;
+  stride = ret->dim[0].stride;
+
+  for (n = 0; n < ncopies; n++)
+    {
+      *dest = *source;
+      dest += stride;
+    }
+}
+
+#endif