re PR fortran/26769 (Implement transpose() and reshape() for real instead of using integer)

	PR libfortran/26769

	* iresolve.c (gfc_resolve_reshape): Call reshape_r4 and
	reshape_r8 instead of reshape_4 and reshape_8.
	(gfc_resolve_transpose): Likewise for transpose.

	* Makefile.am: Add r4 and r8 versions of reshape and transpose.
	* Makefile.in: Regenerate.
	* generated/reshape_r4.c: New file.
	* generated/reshape_r8.c: New file.
	* generated/transpose_r4.c: New file.
	* generated/transpose_r8.c: New file.

From-SVN: r114880

8 files changed, 752 insertions, 30 deletions

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index a614f1d..55e66cd 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,3 +1,10 @@
+2006-06-22  Francois-Xavier Coudert  <coudert@clipper.ens.fr>
+
+	PR libfortran/26769
+	* iresolve.c (gfc_resolve_reshape): Call reshape_r4 and
+	reshape_r8 instead of reshape_4 and reshape_8.
+	(gfc_resolve_transpose): Likewise for transpose.
+
 2006-06-21  Francois-Xavier Coudert  <coudert@clipper.ens.fr>
 
 	* trans-expr.c (gfc_conv_missing_dummy, gfc_conv_unary_op,
diff --git a/gcc/fortran/iresolve.c b/gcc/fortran/iresolve.c
index 3cf84db..f81488a 100644
--- a/gcc/fortran/iresolve.c
+++ b/gcc/fortran/iresolve.c
@@ -1574,14 +1574,10 @@ gfc_resolve_reshape (gfc_expr * f, gfc_expr * source, gfc_expr * shape,
     case 8:
     case 10:
     case 16:
-      if (source->ts.type == BT_COMPLEX)
+      if (source->ts.type == BT_COMPLEX || source->ts.type == BT_REAL)
 	f->value.function.name =
 	  gfc_get_string (PREFIX("reshape_%c%d"),
-			  gfc_type_letter (BT_COMPLEX), source->ts.kind);
-      else if (source->ts.type == BT_REAL && (kind == 10 || kind == 16))
-	f->value.function.name =
-	  gfc_get_string (PREFIX("reshape_%c%d"),
-			  gfc_type_letter (BT_REAL), source->ts.kind);
+			  gfc_type_letter (source->ts.type), source->ts.kind);
       else
 	f->value.function.name =
 	  gfc_get_string (PREFIX("reshape_%d"), source->ts.kind);
@@ -2025,8 +2021,6 @@ gfc_resolve_transfer (gfc_expr * f, gfc_expr * source ATTRIBUTE_UNUSED,
 void
 gfc_resolve_transpose (gfc_expr * f, gfc_expr * matrix)
 {
-  int kind;
-
   f->ts = matrix->ts;
   f->rank = 2;
   if (matrix->shape)
@@ -2036,9 +2030,7 @@ gfc_resolve_transpose (gfc_expr * f, gfc_expr * matrix)
       mpz_init_set (f->shape[1], matrix->shape[0]);
     }
 
-  kind = matrix->ts.kind;
-
-  switch (kind)
+  switch (matrix->ts.kind)
     {
     case 4:
     case 8:
@@ -2046,30 +2038,20 @@ gfc_resolve_transpose (gfc_expr * f, gfc_expr * matrix)
     case 16:
       switch (matrix->ts.type)
         {
+        case BT_REAL:
         case BT_COMPLEX:
           f->value.function.name =
-            gfc_get_string (PREFIX("transpose_c%d"), kind);
+            gfc_get_string (PREFIX("transpose_%c%d"),
+			    gfc_type_letter (matrix->ts.type),
+			    matrix->ts.kind);
           break;
 
-        case BT_REAL:
-	  /* There is no kind=10 integer type and on 32-bit targets
-	     there is usually no kind=16 integer type.  We need to
-	     call the real version.  */
-	  if (kind == 10 || kind == 16)
-	    {
-	      f->value.function.name =
-		gfc_get_string (PREFIX("transpose_r%d"), kind);
-	      break;
-	    }
-
-	  /* Fall through */
-
         case BT_INTEGER:
         case BT_LOGICAL:
 	  /* Use the integer routines for real and logical cases.  This
 	     assumes they all have the same alignment requirements.  */
           f->value.function.name =
-            gfc_get_string (PREFIX("transpose_i%d"), kind);
+            gfc_get_string (PREFIX("transpose_i%d"), matrix->ts.kind);
           break;
 
         default:
diff --git a/libgfortran/Makefile.am b/libgfortran/Makefile.am
index 9c71f4f..22989c8 100644
--- a/libgfortran/Makefile.am
+++ b/libgfortran/Makefile.am
@@ -312,6 +312,8 @@ i_transpose_c= \
 generated/transpose_i4.c \
 generated/transpose_i8.c \
 generated/transpose_i16.c \
+generated/transpose_r4.c \
+generated/transpose_r8.c \
 generated/transpose_r10.c \
 generated/transpose_r16.c \
 generated/transpose_c4.c \
@@ -328,6 +330,8 @@ i_reshape_c= \
 generated/reshape_i4.c \
 generated/reshape_i8.c \
 generated/reshape_i16.c \
+generated/reshape_r4.c \
+generated/reshape_r8.c \
 generated/reshape_r10.c \
 generated/reshape_r16.c \
 generated/reshape_c4.c \
diff --git a/libgfortran/Makefile.in b/libgfortran/Makefile.in
index d7984aa..f8e0bd3 100644
--- a/libgfortran/Makefile.in
+++ b/libgfortran/Makefile.in
@@ -122,15 +122,16 @@ am__objects_16 = matmul_i4.lo matmul_i8.lo matmul_i16.lo matmul_r4.lo \
 	matmul_c8.lo matmul_c10.lo matmul_c16.lo
 am__objects_17 = matmul_l4.lo matmul_l8.lo matmul_l16.lo
 am__objects_18 = transpose_i4.lo transpose_i8.lo transpose_i16.lo \
-	transpose_r10.lo transpose_r16.lo transpose_c4.lo \
-	transpose_c8.lo transpose_c10.lo transpose_c16.lo
+	transpose_r4.lo transpose_r8.lo transpose_r10.lo \
+	transpose_r16.lo transpose_c4.lo transpose_c8.lo \
+	transpose_c10.lo transpose_c16.lo
 am__objects_19 = shape_i4.lo shape_i8.lo shape_i16.lo
 am__objects_20 = eoshift1_4.lo eoshift1_8.lo eoshift1_16.lo
 am__objects_21 = eoshift3_4.lo eoshift3_8.lo eoshift3_16.lo
 am__objects_22 = cshift1_4.lo cshift1_8.lo cshift1_16.lo
 am__objects_23 = reshape_i4.lo reshape_i8.lo reshape_i16.lo \
-	reshape_r10.lo reshape_r16.lo reshape_c4.lo reshape_c8.lo \
-	reshape_c10.lo reshape_c16.lo
+	reshape_r4.lo reshape_r8.lo reshape_r10.lo reshape_r16.lo \
+	reshape_c4.lo reshape_c8.lo reshape_c10.lo reshape_c16.lo
 am__objects_24 = in_pack_i4.lo in_pack_i8.lo in_pack_i16.lo \
 	in_pack_c4.lo in_pack_c8.lo in_pack_c10.lo in_pack_c16.lo
 am__objects_25 = in_unpack_i4.lo in_unpack_i8.lo in_unpack_i16.lo \
@@ -659,6 +660,8 @@ i_transpose_c = \
 generated/transpose_i4.c \
 generated/transpose_i8.c \
 generated/transpose_i16.c \
+generated/transpose_r4.c \
+generated/transpose_r8.c \
 generated/transpose_r10.c \
 generated/transpose_r16.c \
 generated/transpose_c4.c \
@@ -675,6 +678,8 @@ i_reshape_c = \
 generated/reshape_i4.c \
 generated/reshape_i8.c \
 generated/reshape_i16.c \
+generated/reshape_r4.c \
+generated/reshape_r8.c \
 generated/reshape_r10.c \
 generated/reshape_r16.c \
 generated/reshape_c4.c \
@@ -1933,6 +1938,12 @@ transpose_i8.lo: generated/transpose_i8.c
 transpose_i16.lo: generated/transpose_i16.c
 	$(LIBTOOL) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o transpose_i16.lo `test -f 'generated/transpose_i16.c' || echo '$(srcdir)/'`generated/transpose_i16.c
 
+transpose_r4.lo: generated/transpose_r4.c
+	$(LIBTOOL) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o transpose_r4.lo `test -f 'generated/transpose_r4.c' || echo '$(srcdir)/'`generated/transpose_r4.c
+
+transpose_r8.lo: generated/transpose_r8.c
+	$(LIBTOOL) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o transpose_r8.lo `test -f 'generated/transpose_r8.c' || echo '$(srcdir)/'`generated/transpose_r8.c
+
 transpose_r10.lo: generated/transpose_r10.c
 	$(LIBTOOL) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o transpose_r10.lo `test -f 'generated/transpose_r10.c' || echo '$(srcdir)/'`generated/transpose_r10.c
 
@@ -1996,6 +2007,12 @@ reshape_i8.lo: generated/reshape_i8.c
 reshape_i16.lo: generated/reshape_i16.c
 	$(LIBTOOL) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o reshape_i16.lo `test -f 'generated/reshape_i16.c' || echo '$(srcdir)/'`generated/reshape_i16.c
 
+reshape_r4.lo: generated/reshape_r4.c
+	$(LIBTOOL) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o reshape_r4.lo `test -f 'generated/reshape_r4.c' || echo '$(srcdir)/'`generated/reshape_r4.c
+
+reshape_r8.lo: generated/reshape_r8.c
+	$(LIBTOOL) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o reshape_r8.lo `test -f 'generated/reshape_r8.c' || echo '$(srcdir)/'`generated/reshape_r8.c
+
 reshape_r10.lo: generated/reshape_r10.c
 	$(LIBTOOL) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o reshape_r10.lo `test -f 'generated/reshape_r10.c' || echo '$(srcdir)/'`generated/reshape_r10.c
 
diff --git a/libgfortran/generated/reshape_r4.c b/libgfortran/generated/reshape_r4.c
new file mode 100644
index 0000000..acd73a4
--- /dev/null
+++ b/libgfortran/generated/reshape_r4.c
@@ -0,0 +1,257 @@
+/* Implementation of the RESHAPE
+   Copyright 2002, 2006 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public
+License as published by the Free Software Foundation; either
+version 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.)
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+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 "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"
+
+#if defined (HAVE_GFC_REAL_4)
+
+typedef GFC_ARRAY_DESCRIPTOR(1, index_type) shape_type;
+
+/* The shape parameter is ignored. We can currently deduce the shape from the
+   return array.  */
+
+extern void reshape_r4 (gfc_array_r4 * const restrict, 
+	gfc_array_r4 * const restrict, 
+	shape_type * const restrict,
+	gfc_array_r4 * const restrict, 
+	shape_type * const restrict);
+export_proto(reshape_r4);
+
+void
+reshape_r4 (gfc_array_r4 * const restrict ret, 
+	gfc_array_r4 * const restrict source, 
+	shape_type * const restrict shape,
+	gfc_array_r4 * const restrict pad, 
+	shape_type * const restrict order)
+{
+  /* r.* indicates the return array.  */
+  index_type rcount[GFC_MAX_DIMENSIONS];
+  index_type rextent[GFC_MAX_DIMENSIONS];
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type rdim;
+  index_type rsize;
+  index_type rs;
+  index_type rex;
+  GFC_REAL_4 *rptr;
+  /* s.* indicates the source array.  */
+  index_type scount[GFC_MAX_DIMENSIONS];
+  index_type sextent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type sdim;
+  index_type ssize;
+  const GFC_REAL_4 *sptr;
+  /* p.* indicates the pad array.  */
+  index_type pcount[GFC_MAX_DIMENSIONS];
+  index_type pextent[GFC_MAX_DIMENSIONS];
+  index_type pstride[GFC_MAX_DIMENSIONS];
+  index_type pdim;
+  index_type psize;
+  const GFC_REAL_4 *pptr;
+
+  const GFC_REAL_4 *src;
+  int n;
+  int dim;
+
+  if (ret->data == NULL)
+    {
+      rdim = shape->dim[0].ubound - shape->dim[0].lbound + 1;
+      rs = 1;
+      for (n=0; n < rdim; n++)
+	{
+	  ret->dim[n].lbound = 0;
+	  rex = shape->data[n * shape->dim[0].stride];
+	  ret->dim[n].ubound =  rex - 1;
+	  ret->dim[n].stride = rs;
+	  rs *= rex;
+	}
+      ret->offset = 0;
+      ret->data = internal_malloc_size ( rs * sizeof (GFC_REAL_4));
+      ret->dtype = (source->dtype & ~GFC_DTYPE_RANK_MASK) | rdim;
+    }
+  else
+    {
+      rdim = GFC_DESCRIPTOR_RANK (ret);
+    }
+
+  rsize = 1;
+  for (n = 0; n < rdim; n++)
+    {
+      if (order)
+        dim = order->data[n * order->dim[0].stride] - 1;
+      else
+        dim = n;
+
+      rcount[n] = 0;
+      rstride[n] = ret->dim[dim].stride;
+      rextent[n] = ret->dim[dim].ubound + 1 - ret->dim[dim].lbound;
+
+      if (rextent[n] != shape->data[dim * shape->dim[0].stride])
+        runtime_error ("shape and target do not conform");
+
+      if (rsize == rstride[n])
+        rsize *= rextent[n];
+      else
+        rsize = 0;
+      if (rextent[n] <= 0)
+        return;
+    }
+
+  sdim = GFC_DESCRIPTOR_RANK (source);
+  ssize = 1;
+  for (n = 0; n < sdim; n++)
+    {
+      scount[n] = 0;
+      sstride[n] = source->dim[n].stride;
+      sextent[n] = source->dim[n].ubound + 1 - source->dim[n].lbound;
+      if (sextent[n] <= 0)
+        abort ();
+
+      if (ssize == sstride[n])
+        ssize *= sextent[n];
+      else
+        ssize = 0;
+    }
+
+  if (pad)
+    {
+      pdim = GFC_DESCRIPTOR_RANK (pad);
+      psize = 1;
+      for (n = 0; n < pdim; n++)
+        {
+          pcount[n] = 0;
+          pstride[n] = pad->dim[n].stride;
+          pextent[n] = pad->dim[n].ubound + 1 - pad->dim[n].lbound;
+          if (pextent[n] <= 0)
+            abort ();
+          if (psize == pstride[n])
+            psize *= pextent[n];
+          else
+            psize = 0;
+        }
+      pptr = pad->data;
+    }
+  else
+    {
+      pdim = 0;
+      psize = 1;
+      pptr = NULL;
+    }
+
+  if (rsize != 0 && ssize != 0 && psize != 0)
+    {
+      rsize *= sizeof (GFC_REAL_4);
+      ssize *= sizeof (GFC_REAL_4);
+      psize *= sizeof (GFC_REAL_4);
+      reshape_packed ((char *)ret->data, rsize, (char *)source->data,
+		      ssize, pad ? (char *)pad->data : NULL, psize);
+      return;
+    }
+  rptr = ret->data;
+  src = sptr = source->data;
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+
+  while (rptr)
+    {
+      /* Select between the source and pad arrays.  */
+      *rptr = *src;
+      /* Advance to the next element.  */
+      rptr += rstride0;
+      src += sstride0;
+      rcount[0]++;
+      scount[0]++;
+      /* Advance to the next destination element.  */
+      n = 0;
+      while (rcount[n] == rextent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          rcount[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          rptr -= rstride[n] * rextent[n];
+          n++;
+          if (n == rdim)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              rcount[n]++;
+              rptr += rstride[n];
+            }
+        }
+      /* Advance to the next source element.  */
+      n = 0;
+      while (scount[n] == sextent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          scount[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          src -= sstride[n] * sextent[n];
+          n++;
+          if (n == sdim)
+            {
+              if (sptr && pad)
+                {
+                  /* Switch to the pad array.  */
+                  sptr = NULL;
+                  sdim = pdim;
+                  for (dim = 0; dim < pdim; dim++)
+                    {
+                      scount[dim] = pcount[dim];
+                      sextent[dim] = pextent[dim];
+                      sstride[dim] = pstride[dim];
+                      sstride0 = sstride[0];
+                    }
+                }
+              /* We now start again from the beginning of the pad array.  */
+              src = pptr;
+              break;
+            }
+          else
+            {
+              scount[n]++;
+              src += sstride[n];
+            }
+        }
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/reshape_r8.c b/libgfortran/generated/reshape_r8.c
new file mode 100644
index 0000000..5f80f4c
--- /dev/null
+++ b/libgfortran/generated/reshape_r8.c
@@ -0,0 +1,257 @@
+/* Implementation of the RESHAPE
+   Copyright 2002, 2006 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public
+License as published by the Free Software Foundation; either
+version 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.)
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+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 "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"
+
+#if defined (HAVE_GFC_REAL_8)
+
+typedef GFC_ARRAY_DESCRIPTOR(1, index_type) shape_type;
+
+/* The shape parameter is ignored. We can currently deduce the shape from the
+   return array.  */
+
+extern void reshape_r8 (gfc_array_r8 * const restrict, 
+	gfc_array_r8 * const restrict, 
+	shape_type * const restrict,
+	gfc_array_r8 * const restrict, 
+	shape_type * const restrict);
+export_proto(reshape_r8);
+
+void
+reshape_r8 (gfc_array_r8 * const restrict ret, 
+	gfc_array_r8 * const restrict source, 
+	shape_type * const restrict shape,
+	gfc_array_r8 * const restrict pad, 
+	shape_type * const restrict order)
+{
+  /* r.* indicates the return array.  */
+  index_type rcount[GFC_MAX_DIMENSIONS];
+  index_type rextent[GFC_MAX_DIMENSIONS];
+  index_type rstride[GFC_MAX_DIMENSIONS];
+  index_type rstride0;
+  index_type rdim;
+  index_type rsize;
+  index_type rs;
+  index_type rex;
+  GFC_REAL_8 *rptr;
+  /* s.* indicates the source array.  */
+  index_type scount[GFC_MAX_DIMENSIONS];
+  index_type sextent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type sstride0;
+  index_type sdim;
+  index_type ssize;
+  const GFC_REAL_8 *sptr;
+  /* p.* indicates the pad array.  */
+  index_type pcount[GFC_MAX_DIMENSIONS];
+  index_type pextent[GFC_MAX_DIMENSIONS];
+  index_type pstride[GFC_MAX_DIMENSIONS];
+  index_type pdim;
+  index_type psize;
+  const GFC_REAL_8 *pptr;
+
+  const GFC_REAL_8 *src;
+  int n;
+  int dim;
+
+  if (ret->data == NULL)
+    {
+      rdim = shape->dim[0].ubound - shape->dim[0].lbound + 1;
+      rs = 1;
+      for (n=0; n < rdim; n++)
+	{
+	  ret->dim[n].lbound = 0;
+	  rex = shape->data[n * shape->dim[0].stride];
+	  ret->dim[n].ubound =  rex - 1;
+	  ret->dim[n].stride = rs;
+	  rs *= rex;
+	}
+      ret->offset = 0;
+      ret->data = internal_malloc_size ( rs * sizeof (GFC_REAL_8));
+      ret->dtype = (source->dtype & ~GFC_DTYPE_RANK_MASK) | rdim;
+    }
+  else
+    {
+      rdim = GFC_DESCRIPTOR_RANK (ret);
+    }
+
+  rsize = 1;
+  for (n = 0; n < rdim; n++)
+    {
+      if (order)
+        dim = order->data[n * order->dim[0].stride] - 1;
+      else
+        dim = n;
+
+      rcount[n] = 0;
+      rstride[n] = ret->dim[dim].stride;
+      rextent[n] = ret->dim[dim].ubound + 1 - ret->dim[dim].lbound;
+
+      if (rextent[n] != shape->data[dim * shape->dim[0].stride])
+        runtime_error ("shape and target do not conform");
+
+      if (rsize == rstride[n])
+        rsize *= rextent[n];
+      else
+        rsize = 0;
+      if (rextent[n] <= 0)
+        return;
+    }
+
+  sdim = GFC_DESCRIPTOR_RANK (source);
+  ssize = 1;
+  for (n = 0; n < sdim; n++)
+    {
+      scount[n] = 0;
+      sstride[n] = source->dim[n].stride;
+      sextent[n] = source->dim[n].ubound + 1 - source->dim[n].lbound;
+      if (sextent[n] <= 0)
+        abort ();
+
+      if (ssize == sstride[n])
+        ssize *= sextent[n];
+      else
+        ssize = 0;
+    }
+
+  if (pad)
+    {
+      pdim = GFC_DESCRIPTOR_RANK (pad);
+      psize = 1;
+      for (n = 0; n < pdim; n++)
+        {
+          pcount[n] = 0;
+          pstride[n] = pad->dim[n].stride;
+          pextent[n] = pad->dim[n].ubound + 1 - pad->dim[n].lbound;
+          if (pextent[n] <= 0)
+            abort ();
+          if (psize == pstride[n])
+            psize *= pextent[n];
+          else
+            psize = 0;
+        }
+      pptr = pad->data;
+    }
+  else
+    {
+      pdim = 0;
+      psize = 1;
+      pptr = NULL;
+    }
+
+  if (rsize != 0 && ssize != 0 && psize != 0)
+    {
+      rsize *= sizeof (GFC_REAL_8);
+      ssize *= sizeof (GFC_REAL_8);
+      psize *= sizeof (GFC_REAL_8);
+      reshape_packed ((char *)ret->data, rsize, (char *)source->data,
+		      ssize, pad ? (char *)pad->data : NULL, psize);
+      return;
+    }
+  rptr = ret->data;
+  src = sptr = source->data;
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+
+  while (rptr)
+    {
+      /* Select between the source and pad arrays.  */
+      *rptr = *src;
+      /* Advance to the next element.  */
+      rptr += rstride0;
+      src += sstride0;
+      rcount[0]++;
+      scount[0]++;
+      /* Advance to the next destination element.  */
+      n = 0;
+      while (rcount[n] == rextent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          rcount[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          rptr -= rstride[n] * rextent[n];
+          n++;
+          if (n == rdim)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              rcount[n]++;
+              rptr += rstride[n];
+            }
+        }
+      /* Advance to the next source element.  */
+      n = 0;
+      while (scount[n] == sextent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          scount[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          src -= sstride[n] * sextent[n];
+          n++;
+          if (n == sdim)
+            {
+              if (sptr && pad)
+                {
+                  /* Switch to the pad array.  */
+                  sptr = NULL;
+                  sdim = pdim;
+                  for (dim = 0; dim < pdim; dim++)
+                    {
+                      scount[dim] = pcount[dim];
+                      sextent[dim] = pextent[dim];
+                      sstride[dim] = pstride[dim];
+                      sstride0 = sstride[0];
+                    }
+                }
+              /* We now start again from the beginning of the pad array.  */
+              src = pptr;
+              break;
+            }
+          else
+            {
+              scount[n]++;
+              src += sstride[n];
+            }
+        }
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/transpose_r4.c b/libgfortran/generated/transpose_r4.c
new file mode 100644
index 0000000..f441f37
--- /dev/null
+++ b/libgfortran/generated/transpose_r4.c
@@ -0,0 +1,99 @@
+/* Implementation of the TRANSPOSE intrinsic
+   Copyright 2003, 2005, 2006 Free Software Foundation, Inc.
+   Contributed by Tobias Schlüter
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public
+License as published by the Free Software Foundation; either
+version 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.)
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+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 "config.h"
+#include <assert.h>
+#include "libgfortran.h"
+
+#if defined (HAVE_GFC_REAL_4)
+
+extern void transpose_r4 (gfc_array_r4 * const restrict ret, 
+	gfc_array_r4 * const restrict source);
+export_proto(transpose_r4);
+
+void
+transpose_r4 (gfc_array_r4 * const restrict ret, 
+	gfc_array_r4 * const restrict source)
+{
+  /* r.* indicates the return array.  */
+  index_type rxstride, rystride;
+  GFC_REAL_4 *rptr;
+  /* s.* indicates the source array.  */
+  index_type sxstride, systride;
+  const GFC_REAL_4 *sptr;
+
+  index_type xcount, ycount;
+  index_type x, y;
+
+  assert (GFC_DESCRIPTOR_RANK (source) == 2);
+
+  if (ret->data == NULL)
+    {
+      assert (GFC_DESCRIPTOR_RANK (ret) == 2);
+      assert (ret->dtype == source->dtype);
+
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = source->dim[1].ubound - source->dim[1].lbound;
+      ret->dim[0].stride = 1;
+
+      ret->dim[1].lbound = 0;
+      ret->dim[1].ubound = source->dim[0].ubound - source->dim[0].lbound;
+      ret->dim[1].stride = ret->dim[0].ubound+1;
+
+      ret->data = internal_malloc_size (sizeof (GFC_REAL_4) * size0 ((array_t *) ret));
+      ret->offset = 0;
+    }
+
+  sxstride = source->dim[0].stride;
+  systride = source->dim[1].stride;
+  xcount = source->dim[0].ubound + 1 - source->dim[0].lbound;
+  ycount = source->dim[1].ubound + 1 - source->dim[1].lbound;
+
+  rxstride = ret->dim[0].stride;
+  rystride = ret->dim[1].stride;
+
+  rptr = ret->data;
+  sptr = source->data;
+
+  for (y=0; y < ycount; y++)
+    {
+      for (x=0; x < xcount; x++)
+        {
+          *rptr = *sptr;
+
+          sptr += sxstride;
+          rptr += rystride;
+        }
+        sptr += systride - (sxstride * xcount);
+        rptr += rxstride - (rystride * xcount);
+    }
+}
+
+#endif
diff --git a/libgfortran/generated/transpose_r8.c b/libgfortran/generated/transpose_r8.c
new file mode 100644
index 0000000..a97206e
--- /dev/null
+++ b/libgfortran/generated/transpose_r8.c
@@ -0,0 +1,99 @@
+/* Implementation of the TRANSPOSE intrinsic
+   Copyright 2003, 2005, 2006 Free Software Foundation, Inc.
+   Contributed by Tobias Schlüter
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public
+License as published by the Free Software Foundation; either
+version 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.)
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+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 "config.h"
+#include <assert.h>
+#include "libgfortran.h"
+
+#if defined (HAVE_GFC_REAL_8)
+
+extern void transpose_r8 (gfc_array_r8 * const restrict ret, 
+	gfc_array_r8 * const restrict source);
+export_proto(transpose_r8);
+
+void
+transpose_r8 (gfc_array_r8 * const restrict ret, 
+	gfc_array_r8 * const restrict source)
+{
+  /* r.* indicates the return array.  */
+  index_type rxstride, rystride;
+  GFC_REAL_8 *rptr;
+  /* s.* indicates the source array.  */
+  index_type sxstride, systride;
+  const GFC_REAL_8 *sptr;
+
+  index_type xcount, ycount;
+  index_type x, y;
+
+  assert (GFC_DESCRIPTOR_RANK (source) == 2);
+
+  if (ret->data == NULL)
+    {
+      assert (GFC_DESCRIPTOR_RANK (ret) == 2);
+      assert (ret->dtype == source->dtype);
+
+      ret->dim[0].lbound = 0;
+      ret->dim[0].ubound = source->dim[1].ubound - source->dim[1].lbound;
+      ret->dim[0].stride = 1;
+
+      ret->dim[1].lbound = 0;
+      ret->dim[1].ubound = source->dim[0].ubound - source->dim[0].lbound;
+      ret->dim[1].stride = ret->dim[0].ubound+1;
+
+      ret->data = internal_malloc_size (sizeof (GFC_REAL_8) * size0 ((array_t *) ret));
+      ret->offset = 0;
+    }
+
+  sxstride = source->dim[0].stride;
+  systride = source->dim[1].stride;
+  xcount = source->dim[0].ubound + 1 - source->dim[0].lbound;
+  ycount = source->dim[1].ubound + 1 - source->dim[1].lbound;
+
+  rxstride = ret->dim[0].stride;
+  rystride = ret->dim[1].stride;
+
+  rptr = ret->data;
+  sptr = source->data;
+
+  for (y=0; y < ycount; y++)
+    {
+      for (x=0; x < xcount; x++)
+        {
+          *rptr = *sptr;
+
+          sptr += sxstride;
+          rptr += rystride;
+        }
+        sptr += systride - (sxstride * xcount);
+        rptr += rxstride - (rystride * xcount);
+    }
+}
+
+#endif