Merge tree-ssa-20020619-branch into mainline.

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

39 files changed, 3091 insertions, 0 deletions

diff --git a/libgfortran/m4/all.m4 b/libgfortran/m4/all.m4
new file mode 100644
index 00000000000..61a45ac345a
--- /dev/null
+++ b/libgfortran/m4/all.m4
@@ -0,0 +1,37 @@
+`/* Implementation of the ALL intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+ARRAY_FUNCTION(1,
+`  /* Return true only if all the elements are set.  */
+  result = 1;',
+`  if (! *src)
+    {
+      result = 0;
+      break;
+    }')
+
diff --git a/libgfortran/m4/any.m4 b/libgfortran/m4/any.m4
new file mode 100644
index 00000000000..cb79022a9e5
--- /dev/null
+++ b/libgfortran/m4/any.m4
@@ -0,0 +1,37 @@
+`/* Implementation of the ANY intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+ARRAY_FUNCTION(0,
+`  result = 0;',
+`  /* Return true if any of the elements are set.  */
+  if (*src)
+    {
+      result = 1;
+      break;
+    }')
+
diff --git a/libgfortran/m4/cexp.m4 b/libgfortran/m4/cexp.m4
new file mode 100644
index 00000000000..1d22b089a92
--- /dev/null
+++ b/libgfortran/m4/cexp.m4
@@ -0,0 +1,140 @@
+`/* Complex exponential functions
+   Copyright 2002, 2004 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+#include <math.h>
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+/* z = a + ib  */
+/* Absolute value.  */
+real_type
+cabs`'q (complex_type z)
+{
+  return hypot`'q (REALPART (z), IMAGPART (z));
+}
+
+/* Complex argument.  The angle made with the +ve real axis.  Range 0-2pi.  */
+real_type
+carg`'q (complex_type z)
+{
+  real_type arg;
+
+  arg = atan2`'q (IMAGPART (z), REALPART (z));
+  if (arg < 0)
+    return arg + 2 * M_PI;
+  else
+    return arg;
+}
+
+/* exp(z) = exp(a)*(cos(b) + isin(b))  */
+complex_type
+cexp`'q (complex_type z)
+{
+  real_type a;
+  real_type b;
+  complex_type v;
+
+  a = REALPART (z);
+  b = IMAGPART (z);
+  COMPLEX_ASSIGN (v, cos`'q (b), sin`'q (b));
+  return exp`'q (a) * v;
+}
+
+/* log(z) = log (cabs(z)) + i*carg(z)  */
+complex_type
+clog`'q (complex_type z)
+{
+  complex_type v;
+
+  COMPLEX_ASSIGN (v, log`'q (cabs`'q (z)), carg`'q (z));
+  return v;
+}
+
+/* log10(z) = log10 (cabs(z)) + i*carg(z)  */
+complex_type
+clog10`'q (complex_type z)
+{
+  complex_type v;
+
+  COMPLEX_ASSIGN (v, log10`'q (cabs`'q (z)), carg`'q (z));
+  return v;
+}
+
+/* pow(base, power) = cexp (power * clog (base))  */
+complex_type
+cpow`'q (complex_type base, complex_type power)
+{
+  return cexp`'q (power * clog`'q (base));
+}
+
+/* sqrt(z).  Algorithm pulled from glibc.  */
+complex_type
+csqrt`'q (complex_type z)
+{
+  real_type re;
+  real_type im;
+  complex_type v;
+
+  re = REALPART (z);
+  im = IMAGPART (z);
+  if (im == 0.0)
+    {
+      if (re < 0.0)
+        {
+          COMPLEX_ASSIGN (v, 0.0, copysign`'q (sqrt`'q (-re), im));
+        }
+      else
+        {
+          COMPLEX_ASSIGN (v, fabs`'q (sqrt (re)),
+                          copysign`'q (0.0, im));
+        }
+    }
+  else if (re == 0.0)
+    {
+      real_type r;
+
+      r = sqrt`'q (0.5 * fabs (im));
+
+      COMPLEX_ASSIGN (v, copysign`'q (r, im), r);
+    }
+  else
+    {
+      real_type d, r, s;
+
+      d = hypot`'q (re, im);
+      /* Use the identity   2  Re res  Im res = Im x
+         to avoid cancellation error in  d +/- Re x.  */
+      if (re > 0)
+        {
+          r = sqrt`'q (0.5 * d + 0.5 * re);
+          s = (0.5 * im) / r;
+        }
+      else
+        {
+          s = sqrt`'q (0.5 * d - 0.5 * re);
+          r = fabs`'q ((0.5 * im) / s);
+        }
+
+      COMPLEX_ASSIGN (v, r, copysign`'q (s, im));
+    }
+  return v;
+}
+
diff --git a/libgfortran/m4/chyp.m4 b/libgfortran/m4/chyp.m4
new file mode 100644
index 00000000000..f6ee972ea52
--- /dev/null
+++ b/libgfortran/m4/chyp.m4
@@ -0,0 +1,72 @@
+`/* Complex hyperbolic functions
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+#include <math.h>
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+/* Complex number z = a + ib.  */
+
+/* sinh(z) = sinh(a)cos(b) + icosh(a)sin(b)  */
+complex_type
+csinh`'q (complex_type a)
+{
+  real_type r;
+  real_type i;
+  complex_type v;
+
+  r = REALPART (a);
+  i = IMAGPART (a);
+  COMPLEX_ASSIGN (v, sinh`'q (r) * cos`'q (i), cosh`'q (r) * sin`'q (i));
+  return v;
+}
+
+/* cosh(z) = cosh(a)cos(b) - isinh(a)sin(b)  */
+complex_type
+ccosh`'q (complex_type a)
+{
+  real_type r;
+  real_type i;
+  complex_type v;
+
+  r = REALPART (a);
+  i = IMAGPART (a);
+  COMPLEX_ASSIGN (v, cosh`'q (r) * cos`'q (i), - (sinh`'q (r) * sin`'q (i)));
+  return v;
+}
+
+/* tanh(z) = (tanh(a) + itan(b)) / (1 - itanh(a)tan(b))  */
+complex_type
+ctanh`'q (complex_type a)
+{
+  real_type rt;
+  real_type it;
+  complex_type n;
+  complex_type d;
+
+  rt = tanh`'q (REALPART (a));
+  it = tan`'q (IMAGPART (a));
+  COMPLEX_ASSIGN (n, rt, it);
+  COMPLEX_ASSIGN (d, 1, - (rt * it));
+
+  return n / d;
+}
+
diff --git a/libgfortran/m4/count.m4 b/libgfortran/m4/count.m4
new file mode 100644
index 00000000000..a5ea4735a9d
--- /dev/null
+++ b/libgfortran/m4/count.m4
@@ -0,0 +1,33 @@
+`/* Implementation of the COUNT intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+ARRAY_FUNCTION(0,
+`  result = 0;',
+`  if (*src)
+    result++;')
+
diff --git a/libgfortran/m4/cshift1.m4 b/libgfortran/m4/cshift1.m4
new file mode 100644
index 00000000000..382537b7029
--- /dev/null
+++ b/libgfortran/m4/cshift1.m4
@@ -0,0 +1,175 @@
+`/* Implementation of the CSHIFT intrinsic
+   Copyright 2003 Free Software Foundation, Inc.
+   Contributed by Feng Wang <wf_cs@yahoo.com>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(htype_kind, regexp(file, `_\([0-9]+\)\.', `\1'))dnl
+define(htype_code,`i'rtype_name)dnl
+define(htype,get_arraytype(i,htype_kind))dnl
+define(htype_name, get_typename(i,htype_kind))dnl
+
+void
+`__cshift1_'htype_kind (const gfc_array_char * ret, const gfc_array_char * array,
+    const htype * h, const htype_name * pwhich)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS - 1];
+  index_type rstride0;
+  index_type roffset;
+  char *rptr;
+  char *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS - 1];
+  index_type sstride0;
+  index_type soffset;
+  const char *sptr;
+  const char *src;
+`  /* h.* indicates the shift array.  */'
+  index_type hstride[GFC_MAX_DIMENSIONS - 1];
+  index_type hstride0;
+  const htype_name *hptr;
+
+  index_type count[GFC_MAX_DIMENSIONS - 1];
+  index_type extent[GFC_MAX_DIMENSIONS - 1];
+  index_type dim;
+  index_type size;
+  index_type len;
+  index_type n;
+  int which;
+  htype_name sh;
+
+  if (pwhich)
+    which = *pwhich - 1;
+  else
+    which = 0;
+
+  if (which < 0 || (which + 1) > GFC_DESCRIPTOR_RANK (array))
+    runtime_error ("Argument 'DIM' is out of range in call to 'CSHIFT'");
+
+  size = GFC_DESCRIPTOR_SIZE (ret);
+
+  extent[0] = 1;
+  count[0] = 0;
+  size = GFC_DESCRIPTOR_SIZE (array);
+  n = 0;
+
+`/* Initialized for avoiding compiler warnings.  */'
+  roffset = size;
+  soffset = size;
+  len = 0;
+
+  for (dim = 0; dim < GFC_DESCRIPTOR_RANK (array); dim++)
+    {
+      if (dim == which)
+        {
+          roffset = ret->dim[dim].stride * size;
+          if (roffset == 0)
+            roffset = size;
+          soffset = array->dim[dim].stride * size;
+          if (soffset == 0)
+            soffset = size;
+          len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
+        }
+      else
+        {
+          count[n] = 0;
+          extent[n] = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
+          rstride[n] = ret->dim[dim].stride * size;
+          sstride[n] = array->dim[dim].stride * size;
+
+          hstride[n] = h->dim[n].stride;
+          n++;
+        }
+    }
+  if (sstride[0] == 0)
+    sstride[0] = size;
+  if (rstride[0] == 0)
+    rstride[0] = size;
+  if (hstride[0] == 0)
+    hstride[0] = 1;
+
+  dim = GFC_DESCRIPTOR_RANK (array);
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+  hstride0 = hstride[0];
+  rptr = ret->data;
+  sptr = array->data;
+  hptr = h->data;
+
+  while (rptr)
+    {
+`      /* Do the shift for this dimension.  */'
+      sh = *hptr;
+      sh = (div (sh, len)).rem;
+      if (sh < 0)
+        sh += len;
+
+      src = &sptr[sh * soffset];
+      dest = rptr;
+
+      for (n = 0; n < len; n++)
+        {
+          memcpy (dest, src, size);
+          dest += roffset;
+          if (n == len - sh - 1)
+            src = sptr;
+          else
+            src += soffset;
+        }
+
+      /* Advance to the next section.  */
+      rptr += rstride0;
+      sptr += sstride0;
+      hptr += hstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          rptr -= rstride[n] * extent[n];
+          sptr -= sstride[n] * extent[n];
+	  hptr -= hstride[n] * extent[n];
+          n++;
+          if (n >= dim - 1)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              rptr += rstride[n];
+              sptr += sstride[n];
+	      hptr += hstride[n];
+            }
+        }
+    }
+}
+
diff --git a/libgfortran/m4/ctrig.m4 b/libgfortran/m4/ctrig.m4
new file mode 100644
index 00000000000..d35e4500934
--- /dev/null
+++ b/libgfortran/m4/ctrig.m4
@@ -0,0 +1,72 @@
+`/* Complex trig functions
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+#include <math.h>
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+/* Complex number z = a + ib.  */
+
+/* sin(z) = sin(a)cosh(b) + icos(a)sinh(b)  */
+complex_type
+csin`'q (complex_type a)
+{
+  real_type r;
+  real_type i;
+  complex_type v;
+
+  r = REALPART (a);
+  i = IMAGPART (a);
+  COMPLEX_ASSIGN (v, sin`'q (r) * cosh`'q (i), cos`'q (r) * sinh`'q (i));
+  return v;
+}
+
+/* cos(z) = cos(a)cosh(b) - isin(a)sinh(b)  */
+complex_type
+ccos`'q (complex_type a)
+{
+  real_type r;
+  real_type i;
+  complex_type v;
+
+  r = REALPART (a);
+  i = IMAGPART (a);
+  COMPLEX_ASSIGN (v, cos`'q (r) * cosh`'q (i), - (sin`'q (r) * sinh`'q (i)));
+  return v;
+}
+
+/* tan(z) = (tan(a) + itanh(b)) / (1 - itan(a)tanh(b))  */
+complex_type
+ctan`'q (complex_type a)
+{
+  real_type rt;
+  real_type it;
+  complex_type n;
+  complex_type d;
+
+  rt = tan`'q (REALPART (a));
+  it = tanh`'q (IMAGPART (a));
+  COMPLEX_ASSIGN (n, rt, it);
+  COMPLEX_ASSIGN (d , 1, - (rt * it));
+
+  return n / d;
+}
+
diff --git a/libgfortran/m4/dotprod.m4 b/libgfortran/m4/dotprod.m4
new file mode 100644
index 00000000000..051475f6c00
--- /dev/null
+++ b/libgfortran/m4/dotprod.m4
@@ -0,0 +1,71 @@
+`/* Implementation of the DOT_PRODUCT intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_code, regexp(file, `_\([ir][0-9]+\)\.', `\1'))dnl
+define(rtype_letter,substr(rtype_code, 0, 1))dnl
+define(rtype_kind, substr(rtype_code, 1))dnl
+define(rtype,get_arraytype(rtype_letter,rtype_kind))dnl
+define(rtype_name, get_typename(rtype_letter, rtype_kind))dnl
+
+typedef GFC_ARRAY_DESCRIPTOR(GFC_MAX_DIMENSIONS, char) char_array;
+
+/* Both parameters will already have been converted to the result type.  */
+rtype_name
+`__dot_product_'rtype_code (rtype * a, rtype * b)
+{
+  rtype_name *pa;
+  rtype_name *pb;
+  rtype_name res;
+  index_type count;
+  index_type astride;
+  index_type bstride;
+
+  assert (GFC_DESCRIPTOR_RANK (a) == 1
+          && GFC_DESCRIPTOR_RANK (b) == 1);
+
+  if (a->dim[0].stride == 0)
+    a->dim[0].stride = 1;
+  if (b->dim[0].stride == 0)
+    b->dim[0].stride = 1;
+
+  astride = a->dim[0].stride;
+  bstride = b->dim[0].stride;
+  count = a->dim[0].ubound + 1 - a->dim[0].lbound;
+  res = 0;
+  pa = a->data;
+  pb = b->data;
+sinclude(`dotprod_asm_'rtype_code`.m4')dnl
+
+  while (count--)
+    {
+      res += *pa * *pb;
+      pa += astride;
+      pb += bstride;
+    }
+
+  return res;
+}
+
diff --git a/libgfortran/m4/dotprodc.m4 b/libgfortran/m4/dotprodc.m4
new file mode 100644
index 00000000000..0e77c0abfed
--- /dev/null
+++ b/libgfortran/m4/dotprodc.m4
@@ -0,0 +1,74 @@
+`/* Implementation of the DOT_PRODUCT intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+   and Feng Wang <fengwang@nudt.edu.cn>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_code, regexp(file, `_\(c[0-9]+\)\.', `\1'))dnl
+define(rtype_letter,substr(rtype_code, 0, 1))dnl
+define(rtype_kind, substr(rtype_code, 1))dnl
+define(rtype,get_arraytype(rtype_letter,rtype_kind))dnl
+define(rtype_name, get_typename(rtype_letter, rtype_kind))dnl
+
+typedef GFC_ARRAY_DESCRIPTOR(GFC_MAX_DIMENSIONS, char) char_array;
+
+/* Both parameters will already have been converted to the result type.  */
+rtype_name
+`__dot_product_'rtype_code (rtype * a, rtype * b)
+{
+  rtype_name *pa;
+  rtype_name *pb;
+  rtype_name res;
+  rtype_name conjga;
+  index_type count;
+  index_type astride;
+  index_type bstride;
+
+  assert (GFC_DESCRIPTOR_RANK (a) == 1
+          && GFC_DESCRIPTOR_RANK (b) == 1);
+
+  if (a->dim[0].stride == 0)
+    a->dim[0].stride = 1;
+  if (b->dim[0].stride == 0)
+    b->dim[0].stride = 1;
+
+  astride = a->dim[0].stride;
+  bstride = b->dim[0].stride;
+  count = a->dim[0].ubound + 1 - a->dim[0].lbound;
+  res = 0;
+  pa = a->data;
+  pb = b->data;
+sinclude(`dotprod_asm_'rtype_code`.m4')dnl
+
+  while (count--)
+    {
+      COMPLEX_ASSIGN(conjga, REALPART (*pa), -IMAGPART (*pa));
+      res += conjga * *pb;
+      pa += astride;
+      pb += bstride;
+    }
+
+  return res;
+}
+
diff --git a/libgfortran/m4/dotprodl.m4 b/libgfortran/m4/dotprodl.m4
new file mode 100644
index 00000000000..7cbe6000471
--- /dev/null
+++ b/libgfortran/m4/dotprodl.m4
@@ -0,0 +1,79 @@
+`/* Implementation of the DOT_PRODUCT intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_kind, regexp(file, `_l\([0-9]+\)\.', `\1'))dnl
+define(rtype_code,`l'rtype_kind)dnl
+define(rtype,get_arraytype(l,rtype_kind))dnl
+define(rtype_name, get_typename(l, rtype_kind))dnl
+
+rtype_name
+`__dot_product_'rtype_code (gfc_array_l4 * a, gfc_array_l4 * b)
+{
+  GFC_LOGICAL_4 *pa;
+  GFC_LOGICAL_4 *pb;
+  index_type count;
+  index_type astride;
+  index_type bstride;
+
+  assert (GFC_DESCRIPTOR_RANK (a) == 1
+          && GFC_DESCRIPTOR_RANK (b) == 1);
+
+  if (a->dim[0].stride == 0)
+    a->dim[0].stride = 1;
+  if (b->dim[0].stride == 0)
+    b->dim[0].stride = 1;
+
+  astride = a->dim[0].stride;
+  bstride = b->dim[0].stride;
+  count = a->dim[0].ubound + 1 - a->dim[0].lbound;
+
+  pa = a->data;
+  if (GFC_DESCRIPTOR_SIZE (a) != 4)
+    {
+      assert (GFC_DESCRIPTOR_SIZE (a) == 8);
+      pa = GFOR_POINTER_L8_TO_L4 (pa);
+      astride <<= 1;
+    }
+  pb = b->data;
+  if (GFC_DESCRIPTOR_SIZE (b) != 4)
+    {
+      assert (GFC_DESCRIPTOR_SIZE (b) == 8);
+      pb = GFOR_POINTER_L8_TO_L4 (pb);
+      bstride <<= 1;
+    }
+
+  while (count--)
+    {
+      if (*pa && *pb)
+        return 1;
+
+      pa += astride;
+      pb += bstride;
+    }
+
+  return 0;
+}
+
diff --git a/libgfortran/m4/eoshift1.m4 b/libgfortran/m4/eoshift1.m4
new file mode 100644
index 00000000000..304d0038578
--- /dev/null
+++ b/libgfortran/m4/eoshift1.m4
@@ -0,0 +1,183 @@
+`/* Implementation of the EOSHIFT intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(htype_kind, regexp(file, `_\([0-9]+\)\.', `\1'))dnl
+define(htype_code,`i'rtype_name)dnl
+define(htype,get_arraytype(i,htype_kind))dnl
+define(htype_name, get_typename(i,htype_kind))dnl
+
+static const char zeros[16] =
+  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+
+void
+`__eoshift1_'htype_kind (const gfc_array_char * ret, const gfc_array_char * array,
+    const htype * h, const char * pbound, const htype_name * pwhich)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS - 1];
+  index_type rstride0;
+  index_type roffset;
+  char *rptr;
+  char *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS - 1];
+  index_type sstride0;
+  index_type soffset;
+  const char *sptr;
+  const char *src;
+`  /* h.* indicates the shift array.  */'
+  index_type hstride[GFC_MAX_DIMENSIONS - 1];
+  index_type hstride0;
+  const htype_name *hptr;
+
+  index_type count[GFC_MAX_DIMENSIONS - 1];
+  index_type extent[GFC_MAX_DIMENSIONS - 1];
+  index_type dim;
+  index_type size;
+  index_type len;
+  index_type n;
+  int which;
+  htype_name sh;
+  htype_name delta;
+
+  if (pwhich)
+    which = *pwhich - 1;
+  else
+    which = 0;
+
+  if (!pbound)
+    pbound = zeros;
+
+  size = GFC_DESCRIPTOR_SIZE (ret);
+
+  extent[0] = 1;
+  count[0] = 0;
+  size = GFC_DESCRIPTOR_SIZE (array);
+  n = 0;
+  for (dim = 0; dim < GFC_DESCRIPTOR_RANK (array); dim++)
+    {
+      if (dim == which)
+        {
+          roffset = ret->dim[dim].stride * size;
+          if (roffset == 0)
+            roffset = size;
+          soffset = array->dim[dim].stride * size;
+          if (soffset == 0)
+            soffset = size;
+          len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
+        }
+      else
+        {
+          count[n] = 0;
+          extent[n] = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
+          rstride[n] = ret->dim[dim].stride * size;
+          sstride[n] = array->dim[dim].stride * size;
+
+          hstride[n] = h->dim[n].stride;
+          n++;
+        }
+    }
+  if (sstride[0] == 0)
+    sstride[0] = size;
+  if (rstride[0] == 0)
+    rstride[0] = size;
+  if (hstride[0] == 0)
+    hstride[0] = 1;
+
+  dim = GFC_DESCRIPTOR_RANK (array);
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+  hstride0 = hstride[0];
+  rptr = ret->data;
+  sptr = array->data;
+  hptr = h->data;
+
+  while (rptr)
+    {
+`      /* Do the shift for this dimension.  */'
+      sh = *hptr;
+      delta = (sh >= 0) ? sh: -sh;
+      if (sh > 0)
+        {
+          src = &sptr[delta * soffset];
+          dest = rptr;
+        }
+      else
+        {
+          src = sptr;
+          dest = &rptr[delta * roffset];
+        }
+      for (n = 0; n < len - delta; n++)
+        {
+          memcpy (dest, src, size);
+          dest += roffset;
+          src += soffset;
+        }
+      if (sh < 0)
+        dest = rptr;
+      n = delta;
+
+      while (n--)
+        {
+          memcpy (dest, pbound, size);
+          dest += roffset;
+        }
+
+      /* Advance to the next section.  */
+      rptr += rstride0;
+      sptr += sstride0;
+      hptr += hstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          rptr -= rstride[n] * extent[n];
+          sptr -= sstride[n] * extent[n];
+	  hptr -= hstride[n] * extent[n];
+          n++;
+          if (n >= dim - 1)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              rptr += rstride[n];
+              sptr += sstride[n];
+	      hptr += hstride[n];
+            }
+        }
+    }
+}
+
diff --git a/libgfortran/m4/eoshift3.m4 b/libgfortran/m4/eoshift3.m4
new file mode 100644
index 00000000000..b86a80c41ca
--- /dev/null
+++ b/libgfortran/m4/eoshift3.m4
@@ -0,0 +1,198 @@
+`/* Implementation of the EOSHIFT intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(htype_kind, regexp(file, `_\([0-9]+\)\.', `\1'))dnl
+define(htype_code,`i'rtype_name)dnl
+define(htype,get_arraytype(i,htype_kind))dnl
+define(htype_name, get_typename(i,htype_kind))dnl
+
+static const char zeros[16] =
+  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+
+void
+`__eoshift3_'htype_kind (gfc_array_char * ret, gfc_array_char * array,
+    htype * h, const gfc_array_char * bound, htype_name * pwhich)
+{
+  /* r.* indicates the return array.  */
+  index_type rstride[GFC_MAX_DIMENSIONS - 1];
+  index_type rstride0;
+  index_type roffset;
+  char *rptr;
+  char *dest;
+  /* s.* indicates the source array.  */
+  index_type sstride[GFC_MAX_DIMENSIONS - 1];
+  index_type sstride0;
+  index_type soffset;
+  const char *sptr;
+  const char *src;
+`  /* h.* indicates the shift array.  */'
+  index_type hstride[GFC_MAX_DIMENSIONS - 1];
+  index_type hstride0;
+  const htype_name *hptr;
+  /* b.* indicates the bound array.  */
+  index_type bstride[GFC_MAX_DIMENSIONS - 1];
+  index_type bstride0;
+  const char *bptr;
+
+  index_type count[GFC_MAX_DIMENSIONS - 1];
+  index_type extent[GFC_MAX_DIMENSIONS - 1];
+  index_type dim;
+  index_type size;
+  index_type len;
+  index_type n;
+  int which;
+  htype_name sh;
+  htype_name delta;
+
+  if (pwhich)
+    which = *pwhich - 1;
+  else
+    which = 0;
+
+  size = GFC_DESCRIPTOR_SIZE (ret);
+
+  extent[0] = 1;
+  count[0] = 0;
+  size = GFC_DESCRIPTOR_SIZE (array);
+  n = 0;
+  for (dim = 0; dim < GFC_DESCRIPTOR_RANK (array); dim++)
+    {
+      if (dim == which)
+        {
+          roffset = ret->dim[dim].stride * size;
+          if (roffset == 0)
+            roffset = size;
+          soffset = array->dim[dim].stride * size;
+          if (soffset == 0)
+            soffset = size;
+          len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
+        }
+      else
+        {
+          count[n] = 0;
+          extent[n] = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
+          rstride[n] = ret->dim[dim].stride * size;
+          sstride[n] = array->dim[dim].stride * size;
+
+          hstride[n] = h->dim[n].stride;
+          if (bound)
+            bstride[n] = bound->dim[n].stride;
+          else
+            bstride[n] = 0;
+          n++;
+        }
+    }
+  if (sstride[0] == 0)
+    sstride[0] = size;
+  if (rstride[0] == 0)
+    rstride[0] = size;
+  if (hstride[0] == 0)
+    hstride[0] = 1;
+  if (bound && bstride[0] == 0)
+    bstride[0] = size;
+
+  dim = GFC_DESCRIPTOR_RANK (array);
+  rstride0 = rstride[0];
+  sstride0 = sstride[0];
+  hstride0 = hstride[0];
+  bstride0 = bstride[0];
+  rptr = ret->data;
+  sptr = array->data;
+  hptr = h->data;
+  if (bound)
+    bptr = bound->data;
+  else
+    bptr = zeros;
+
+  while (rptr)
+    {
+`      /* Do the shift for this dimension.  */'
+      sh = *hptr;
+      delta = (sh >= 0) ? sh: -sh;
+      if (sh > 0)
+        {
+          src = &sptr[delta * soffset];
+          dest = rptr;
+        }
+      else
+        {
+          src = sptr;
+          dest = &rptr[delta * roffset];
+        }
+      for (n = 0; n < len - delta; n++)
+        {
+          memcpy (dest, src, size);
+          dest += roffset;
+          src += soffset;
+        }
+      if (sh < 0)
+        dest = rptr;
+      n = delta;
+
+      while (n--)
+        {
+          memcpy (dest, bptr, size);
+          dest += roffset;
+        }
+
+      /* Advance to the next section.  */
+      rptr += rstride0;
+      sptr += sstride0;
+      hptr += hstride0;
+      bptr += bstride0;
+      count[0]++;
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          rptr -= rstride[n] * extent[n];
+          sptr -= sstride[n] * extent[n];
+	  hptr -= hstride[n] * extent[n];
+          bptr -= bstride[n] * extent[n];
+          n++;
+          if (n >= dim - 1)
+            {
+              /* Break out of the loop.  */
+              rptr = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              rptr += rstride[n];
+              sptr += sstride[n];
+	      hptr += hstride[n];
+              bptr += bstride[n];
+            }
+        }
+    }
+}
+
diff --git a/libgfortran/m4/exponent.m4 b/libgfortran/m4/exponent.m4
new file mode 100644
index 00000000000..510f763292d
--- /dev/null
+++ b/libgfortran/m4/exponent.m4
@@ -0,0 +1,32 @@
+`/* Implementation of the EXPONENT intrinsic
+   Copyright 2003 Free Software Foundation, Inc.
+   Contributed by Richard Henderson <rth@redhat.com>.
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+#include <math.h>
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+GFC_INTEGER_4
+prefix(exponent_r`'kind) (real_type s)
+{
+  int ret;
+  frexp`'q (s, &ret);
+  return ret;
+}
diff --git a/libgfortran/m4/fraction.m4 b/libgfortran/m4/fraction.m4
new file mode 100644
index 00000000000..c453e78ac41
--- /dev/null
+++ b/libgfortran/m4/fraction.m4
@@ -0,0 +1,31 @@
+`/* Implementation of the FRACTION intrinsic
+   Copyright 2003 Free Software Foundation, Inc.
+   Contributed by Richard Henderson <rth@redhat.com>.
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+#include <math.h>
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+real_type
+prefix(fraction_r`'kind) (real_type s)
+{
+  int dummy_exp;
+  return frexp`'q (s, &dummy_exp);
+}
diff --git a/libgfortran/m4/head.m4 b/libgfortran/m4/head.m4
new file mode 100644
index 00000000000..0b2f9eeaa63
--- /dev/null
+++ b/libgfortran/m4/head.m4
@@ -0,0 +1,21 @@
+`!   Copyright 2002 Free Software Foundation, Inc.
+!   Contributed by Paul Brook <paul@nowt.org>
+!
+!This file is part of the GNU Fortran 95 runtime library (libgfor).
+!
+!GNU libgfor is free software; you can redistribute it and/or
+!modify it under the terms of the GNU Lesser General Public
+!License as published by the Free Software Foundation; either
+!version 2.1 of the License, or (at your option) any later version.
+!
+!GNU libgfor 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 Lesser General Public License for more details.
+!
+!You should have received a copy of the GNU Lesser General Public
+!License along with libgfor; see the file COPYING.  If not,
+!write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+!Boston, MA 02111-1307, USA.
+!
+!This file is machine generated.'
diff --git a/libgfortran/m4/iforeach.m4 b/libgfortran/m4/iforeach.m4
new file mode 100644
index 00000000000..2397036ba72
--- /dev/null
+++ b/libgfortran/m4/iforeach.m4
@@ -0,0 +1,196 @@
+dnl Support macro file for intrinsic functions.
+dnl Contains the generic sections of the array functions.
+dnl This file is part of the GNU Fortran 95 Runtime Library (libgfortran)
+dnl Distributed under the GNU LGPL.  See COPYING for details.
+define(START_FOREACH_FUNCTION,
+`void
+`__'name`'rtype_qual`_'type_code (rtype * retarray, atype *array)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type dstride;
+  type_name *base;
+  rtype_name *dest;
+  index_type rank;
+  index_type n;
+
+  rank = GFC_DESCRIPTOR_RANK (array);
+  assert (rank > 0);
+  assert (GFC_DESCRIPTOR_RANK (retarray) == 1);
+  assert (retarray->dim[0].ubound + 1 - retarray->dim[0].lbound == rank);
+  if (array->dim[0].stride == 0)
+    array->dim[0].stride = 1;
+  if (retarray->dim[0].stride == 0)
+    retarray->dim[0].stride = 1;
+
+  dstride = retarray->dim[0].stride;
+  dest = retarray->data;
+  for (n = 0; n < rank; n++)
+    {
+      sstride[n] = array->dim[n].stride;
+      extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
+      count[n] = 0;
+      if (extent[n] <= 0)
+	{
+	  /* Set the return value.  */
+	  for (n = 0; n < rank; n++)
+	    dest[n * dstride] = 0;
+	  return;
+	}
+    }
+
+  base = array->data;
+
+  /* Initialize the return value.  */
+  for (n = 0; n < rank; n++)
+    dest[n * dstride] = 1;
+  {
+')dnl
+define(START_FOREACH_BLOCK,
+`  while (base)
+    {
+      {
+        /* Implementation start.  */
+')dnl
+define(FINISH_FOREACH_FUNCTION,
+`        /* Implementation end.  */
+      }
+      /* Advance to the next element.  */
+      count[0]++;
+      base += sstride[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 proabably not worth it.  */
+          base -= sstride[n] * extent[n];
+          n++;
+          if (n == rank)
+            {
+              /* Break out of the loop.  */
+              base = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              base += sstride[n];
+            }
+        }
+    }
+  }
+}')dnl
+define(START_MASKED_FOREACH_FUNCTION,
+`void
+`__m'name`'rtype_qual`_'type_code (rtype * retarray, atype *array, gfc_array_l4 * mask)
+{
+  index_type count[GFC_MAX_DIMENSIONS];
+  index_type extent[GFC_MAX_DIMENSIONS];
+  index_type sstride[GFC_MAX_DIMENSIONS];
+  index_type mstride[GFC_MAX_DIMENSIONS];
+  index_type dstride;
+  rtype_name *dest;
+  type_name *base;
+  GFC_LOGICAL_4 *mbase;
+  int rank;
+  index_type n;
+
+  rank = GFC_DESCRIPTOR_RANK (array);
+  assert (rank > 0);
+  assert (GFC_DESCRIPTOR_RANK (retarray) == 1);
+  assert (retarray->dim[0].ubound + 1 - retarray->dim[0].lbound == rank);
+  assert (GFC_DESCRIPTOR_RANK (mask) == rank);
+
+  if (array->dim[0].stride == 0)
+    array->dim[0].stride = 1;
+  if (retarray->dim[0].stride == 0)
+    retarray->dim[0].stride = 1;
+  if (retarray->dim[0].stride == 0)
+    retarray->dim[0].stride = 1;
+
+  dstride = retarray->dim[0].stride;
+  dest = retarray->data;
+  for (n = 0; n < rank; n++)
+    {
+      sstride[n] = array->dim[n].stride;
+      mstride[n] = mask->dim[n].stride;
+      extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
+      count[n] = 0;
+      if (extent[n] <= 0)
+	{
+	  /* Set the return value.  */
+	  for (n = 0; n < rank; n++)
+	    dest[n * dstride] = 0;
+	  return;
+	}
+    }
+
+  base = array->data;
+  mbase = mask->data;
+
+  if (GFC_DESCRIPTOR_SIZE (mask) != 4)
+    {
+      /* This allows the same loop to be used for all logical types.  */
+      assert (GFC_DESCRIPTOR_SIZE (mask) == 8);
+      for (n = 0; n < rank; n++)
+        mstride[n] <<= 1;
+      mbase = (GFOR_POINTER_L8_TO_L4 (mbase));
+    }
+
+
+  /* Initialize the return value.  */
+  for (n = 0; n < rank; n++)
+    dest[n * dstride] = 1;
+  {
+')dnl
+define(START_MASKED_FOREACH_BLOCK, `START_FOREACH_BLOCK')dnl
+define(FINISH_MASKED_FOREACH_FUNCTION,
+`        /* Implementation end.  */
+      }
+      /* Advance to the next element.  */
+      count[0]++;
+      base += sstride[0];
+      mbase += mstride[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 proabably not worth it.  */
+          base -= sstride[n] * extent[n];
+          mbase -= mstride[n] * extent[n];
+          n++;
+          if (n == rank)
+            {
+              /* Break out of the loop.  */
+              base = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              base += sstride[n];
+              mbase += mstride[n];
+            }
+        }
+    }
+  }
+}')dnl
+define(FOREACH_FUNCTION,
+`START_FOREACH_FUNCTION
+$1
+START_FOREACH_BLOCK
+$2
+FINISH_FOREACH_FUNCTION')dnl
+define(MASKED_FOREACH_FUNCTION,
+`START_MASKED_FOREACH_FUNCTION
+$1
+START_MASKED_FOREACH_BLOCK
+$2
+FINISH_MASKED_FOREACH_FUNCTION')dnl
diff --git a/libgfortran/m4/ifunction.m4 b/libgfortran/m4/ifunction.m4
new file mode 100644
index 00000000000..95445846185
--- /dev/null
+++ b/libgfortran/m4/ifunction.m4
@@ -0,0 +1,256 @@
+dnl Support macro file for intrinsic functions.
+dnl Contains the generic sections of the array functions.
+dnl This file is part of the GNU Fortran 95 Runtime Library (libgfortran)
+dnl Distributed under the GNU LGPL.  See COPYING for details.
+dnl
+dnl Pass the implementation for a single section as the parameter to
+dnl {MASK_}ARRAY_FUNCTION.
+dnl The variables base, delta, and len describe the input section.
+dnl For masked section the mask is described by mbase and mdelta.
+dnl These should not be modified. The result should be stored in *dest.
+dnl The names count, extent, sstride, dstride, base, dest, rank, dim
+dnl retarray, array, pdim and mstride should not be used.
+dnl The variable n is declared as index_type and may be used.
+dnl Other variable declarations may be placed at the start of the code,
+dnl The types of the array parameter and the return value are
+dnl type_name and rtype_name respectively.
+dnl Execution should be allowed to continue to the end of the block.
+dnl You should not return or break from the inner loop of the implementation.
+dnl Care should also be taken to avoid using the names defined in iparm.m4
+define(START_ARRAY_FUNCTION,
+`void
+`__'name`'rtype_qual`_'type_code (rtype * retarray, atype *array, index_type *pdim)
+{
+  index_type count[GFC_MAX_DIMENSIONS - 1];
+  index_type extent[GFC_MAX_DIMENSIONS - 1];
+  index_type sstride[GFC_MAX_DIMENSIONS - 1];
+  index_type dstride[GFC_MAX_DIMENSIONS - 1];
+  type_name *base;
+  rtype_name *dest;
+  index_type rank;
+  index_type n;
+  index_type len;
+  index_type delta;
+  index_type dim;
+
+  /* Make dim zero based to avoid confusion.  */
+  dim = (*pdim) - 1;
+  rank = GFC_DESCRIPTOR_RANK (array) - 1;
+  assert (rank == GFC_DESCRIPTOR_RANK (retarray));
+  if (array->dim[0].stride == 0)
+    array->dim[0].stride = 1;
+  if (retarray->dim[0].stride == 0)
+    retarray->dim[0].stride = 1;
+
+  len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
+  delta = array->dim[dim].stride;
+
+  for (n = 0; n < dim; n++)
+    {
+      sstride[n] = array->dim[n].stride;
+      extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
+    }
+  for (n = dim; n < rank; n++)
+    {
+      sstride[n] = array->dim[n + 1].stride;
+      extent[n] =
+        array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
+    }
+
+  for (n = 0; n < rank; n++)
+    {
+      count[n] = 0;
+      dstride[n] = retarray->dim[n].stride;
+      if (extent[n] <= 0)
+        len = 0;
+    }
+
+  base = array->data;
+  dest = retarray->data;
+
+  while (base)
+    {
+      type_name *src;
+      rtype_name result;
+      src = base;
+      {
+')dnl
+define(START_ARRAY_BLOCK,
+`        if (len <= 0)
+	  *dest = '$1`;
+	else
+	  {
+	    for (n = 0; n < len; n++, src += delta)
+	      {
+')dnl
+define(FINISH_ARRAY_FUNCTION,
+    `          }
+	    *dest = result;
+	  }
+      }
+      /* Advance to the next element.  */
+      count[0]++;
+      base += sstride[0];
+      dest += dstride[0];
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          base -= sstride[n] * extent[n];
+          dest -= dstride[n] * extent[n];
+          n++;
+          if (n == rank)
+            {
+              /* Break out of the look.  */
+              base = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              base += sstride[n];
+              dest += dstride[n];
+            }
+        }
+    }
+}')dnl
+define(START_MASKED_ARRAY_FUNCTION,
+`void
+`__m'name`'rtype_qual`_'type_code (rtype * retarray, atype * array, index_type *pdim, gfc_array_l4 * mask)
+{
+  index_type count[GFC_MAX_DIMENSIONS - 1];
+  index_type extent[GFC_MAX_DIMENSIONS - 1];
+  index_type sstride[GFC_MAX_DIMENSIONS - 1];
+  index_type dstride[GFC_MAX_DIMENSIONS - 1];
+  index_type mstride[GFC_MAX_DIMENSIONS - 1];
+  rtype_name *dest;
+  type_name *base;
+  GFC_LOGICAL_4 *mbase;
+  int rank;
+  int dim;
+  index_type n;
+  index_type len;
+  index_type delta;
+  index_type mdelta;
+
+  dim = (*pdim) - 1;
+  rank = GFC_DESCRIPTOR_RANK (array) - 1;
+  assert (rank == GFC_DESCRIPTOR_RANK (retarray));
+  if (array->dim[0].stride == 0)
+    array->dim[0].stride = 1;
+  if (retarray->dim[0].stride == 0)
+    retarray->dim[0].stride = 1;
+
+  len = array->dim[dim].ubound + 1 - array->dim[dim].lbound;
+  if (len <= 0)
+    return;
+  delta = array->dim[dim].stride;
+  mdelta = mask->dim[dim].stride;
+
+  for (n = 0; n < dim; n++)
+    {
+      sstride[n] = array->dim[n].stride;
+      mstride[n] = mask->dim[n].stride;
+      extent[n] = array->dim[n].ubound + 1 - array->dim[n].lbound;
+    }
+  for (n = dim; n < rank; n++)
+    {
+      sstride[n] = array->dim[n + 1].stride;
+      mstride[n] = mask->dim[n + 1].stride;
+      extent[n] =
+        array->dim[n + 1].ubound + 1 - array->dim[n + 1].lbound;
+    }
+
+  for (n = 0; n < rank; n++)
+    {
+      count[n] = 0;
+      dstride[n] = retarray->dim[n].stride;
+      if (extent[n] <= 0)
+        return;
+    }
+
+  dest = retarray->data;
+  base = array->data;
+  mbase = mask->data;
+
+  if (GFC_DESCRIPTOR_SIZE (mask) != 4)
+    {
+      /* This allows the same loop to be used for all logical types.  */
+      assert (GFC_DESCRIPTOR_SIZE (mask) == 8);
+      for (n = 0; n < rank; n++)
+        mstride[n] <<= 1;
+      mdelta <<= 1;
+      mbase = (GFOR_POINTER_L8_TO_L4 (mbase));
+    }
+
+  while (base)
+    {
+      type_name *src;
+      GFC_LOGICAL_4 *msrc;
+      rtype_name result;
+      src = base;
+      msrc = mbase;
+      {
+')dnl
+define(START_MASKED_ARRAY_BLOCK,
+`        if (len <= 0)
+	  *dest = '$1`;
+	else
+	  {
+	    for (n = 0; n < len; n++, src += delta, msrc += mdelta)
+	      {
+')dnl
+define(FINISH_MASKED_ARRAY_FUNCTION,
+`              }
+	    *dest = result;
+	  }
+      }
+      /* Advance to the next element.  */
+      count[0]++;
+      base += sstride[0];
+      mbase += mstride[0];
+      dest += dstride[0];
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          base -= sstride[n] * extent[n];
+          mbase -= mstride[n] * extent[n];
+          dest -= dstride[n] * extent[n];
+          n++;
+          if (n == rank)
+            {
+              /* Break out of the look.  */
+              base = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              base += sstride[n];
+              mbase += mstride[n];
+              dest += dstride[n];
+            }
+        }
+    }
+}')dnl
+define(ARRAY_FUNCTION,
+`START_ARRAY_FUNCTION
+$2
+START_ARRAY_BLOCK($1)
+$3
+FINISH_ARRAY_FUNCTION')dnl
+define(MASKED_ARRAY_FUNCTION,
+`START_MASKED_ARRAY_FUNCTION
+$2
+START_MASKED_ARRAY_BLOCK($1)
+$3
+FINISH_MASKED_ARRAY_FUNCTION')dnl
diff --git a/libgfortran/m4/in_pack.m4 b/libgfortran/m4/in_pack.m4
new file mode 100644
index 00000000000..4998ed5f2db
--- /dev/null
+++ b/libgfortran/m4/in_pack.m4
@@ -0,0 +1,122 @@
+`/* Helper function for repacking arrays.
+   Copyright 2003 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfortran).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfortran; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_kind, regexp(file, `_.\([0-9]+\)\.', `\1'))dnl
+define(rtype_letter, regexp(file, `_\(.\)[0-9]+\.', `\1'))dnl
+define(rtype_code,rtype_letter`'rtype_name)dnl
+define(rtype,get_arraytype(rtype_letter,rtype_kind))dnl
+define(rtype_name, get_typename(rtype_letter, rtype_kind))dnl
+
+
+/* Allocates a block of memory with internal_malloc if the array needs
+   repacking.  */
+
+dnl Only the kind (ie size) is used to name the function.
+rtype_name *
+`internal_pack_'rtype_kind (rtype * source)
+{
+  index_type count[GFC_MAX_DIMENSIONS - 1];
+  index_type extent[GFC_MAX_DIMENSIONS - 1];
+  index_type stride[GFC_MAX_DIMENSIONS - 1];
+  index_type stride0;
+  index_type dim;
+  index_type ssize;
+  const rtype_name *src;
+  rtype_name *dest;
+  rtype_name *destptr;
+  int n;
+  int packed;
+
+  if (source->dim[0].stride == 0)
+    {
+      source->dim[0].stride = 1;
+      return source->data;
+    }
+
+  dim = GFC_DESCRIPTOR_RANK (source);
+  ssize = 1;
+  packed = 1;
+  for (n = 0; n < dim; n++)
+    {
+      count[n] = 0;
+      stride[n] = source->dim[n].stride;
+      extent[n] = source->dim[n].ubound + 1 - source->dim[n].lbound;
+      if (extent[n] <= 0)
+        {
+          /* Do nothing.  */
+          packed = 1;
+          break;
+        }
+
+      if (ssize != stride[n])
+        packed = 0;
+
+      ssize *= extent[n];
+    }
+
+  if (packed)
+    return source->data;
+
+  /* Allocate storage for the destination.  */
+  destptr = (rtype_name *)internal_malloc_size (ssize * rtype_kind);
+  dest = destptr;
+  src = source->data;
+  stride0 = stride[0];
+
+
+  while (src)
+    {
+      /* Copy the data.  */
+      *(dest++) = *src;
+      /* Advance to the next element.  */
+      src += stride0;
+      count[0]++;
+      /* Advance to the next source element.  */
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          src -= stride[n] * extent[n];
+          n++;
+          if (n == dim)
+            {
+              src = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              src += stride[n];
+            }
+        }
+    }
+  return destptr;
+}
+
diff --git a/libgfortran/m4/in_unpack.m4 b/libgfortran/m4/in_unpack.m4
new file mode 100644
index 00000000000..fe344ca7d31
--- /dev/null
+++ b/libgfortran/m4/in_unpack.m4
@@ -0,0 +1,109 @@
+`/* Helper function for repacking arrays.
+   Copyright 2003 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_kind, regexp(file, `_.\([0-9]+\)\.', `\1'))dnl
+define(rtype_letter, regexp(file, `_\(.\)[0-9]+\.', `\1'))dnl
+define(rtype_code,rtype_letter`'rtype_name)dnl
+define(rtype,get_arraytype(rtype_letter,rtype_kind))dnl
+define(rtype_name, get_typename(rtype_letter, rtype_kind))dnl
+
+dnl Only the kind (ie size) is used to name the function.
+void
+`internal_unpack_'rtype_kind (rtype * d, const rtype_name * src)
+{
+  index_type count[GFC_MAX_DIMENSIONS - 1];
+  index_type extent[GFC_MAX_DIMENSIONS - 1];
+  index_type stride[GFC_MAX_DIMENSIONS - 1];
+  index_type stride0;
+  index_type dim;
+  index_type dsize;
+  rtype_name *dest;
+  int n;
+
+  dest = d->data;
+  if (src == dest || !src)
+    return;
+
+  if (d->dim[0].stride == 0)
+    d->dim[0].stride = 1;
+
+  dim = GFC_DESCRIPTOR_RANK (d);
+  dsize = 1;
+  for (n = 0; n < dim; n++)
+    {
+      count[n] = 0;
+      stride[n] = d->dim[n].stride;
+      extent[n] = d->dim[n].ubound + 1 - d->dim[n].lbound;
+      if (extent[n] <= 0)
+        abort ();
+
+      if (dsize == stride[n])
+        dsize *= extent[n];
+      else
+        dsize = 0;
+    }
+
+  if (dsize != 0)
+    {
+      memcpy (dest, src, dsize * rtype_kind);
+      return;
+    }
+
+  stride0 = stride[0];
+
+  while (dest)
+    {
+      /* Copy the data.  */
+      *dest = *(src++);
+      /* Advance to the next element.  */
+      dest += stride0;
+      count[0]++;
+      /* Advance to the next source element.  */
+      n = 0;
+      while (count[n] == extent[n])
+        {
+          /* When we get to the end of a dimension, reset it and increment
+             the next dimension.  */
+          count[n] = 0;
+          /* We could precalculate these products, but this is a less
+             frequently used path so proabably not worth it.  */
+          dest -= stride[n] * extent[n];
+          n++;
+          if (n == dim)
+            {
+              dest = NULL;
+              break;
+            }
+          else
+            {
+              count[n]++;
+              dest += stride[n];
+            }
+        }
+    }
+}
+
diff --git a/libgfortran/m4/iparm.m4 b/libgfortran/m4/iparm.m4
new file mode 100644
index 00000000000..6cbd7b247a8
--- /dev/null
+++ b/libgfortran/m4/iparm.m4
@@ -0,0 +1,26 @@
+dnl Support macro file for intrinsic functions.
+dnl Works out all the function types from the filename.
+dnl This file is part of the GNU Fortran 95 Runtime Library (libgfortran)
+dnl Distributed under the GNU LGPL.  See COPYING for details.
+dnl M4 macro file to get type names from filenames
+include(`types.m4')
+define(type_letter, regexp(file, `_\([irlc]\)[^_]*$', \1))dnl
+define(type_kind, regexp(file, `_[irlc]\([0-9]*\)[^_]*$', \1))dnl
+define(rtype_kind, regexp(file, `_\([0-9]*\)_[irlc][0-9]*[^_]*$', `\1'))dnl
+define(type_code, type_letter`'type_kind)dnl
+define(type_name, get_typename(type_letter,type_kind))dnl
+define(atype, get_arraytype(type_letter,type_kind))dnl
+ifelse(rtype_kind,,
+`define(rtype_letter,type_letter)dnl
+define(rtype_name, type_name)dnl
+define(rtype_code, type_code)dnl
+define(rtype, atype)dnl
+define(rtype_qual,`')dnl
+',
+`define(rtype_letter,i)dnl
+define(rtype_name, get_typename(rtype_letter,rtype_kind))dnl
+define(rtype, get_arraytype(rtype_letter,rtype_kind))dnl
+define(rtype_qual,`_'rtype_kind)dnl
+')dnl
+define(type_max, type_name`_HUGE')dnl
+define(type_min, `-'type_max)dnl
diff --git a/libgfortran/m4/matmul.m4 b/libgfortran/m4/matmul.m4
new file mode 100644
index 00000000000..fcf63a6504e
--- /dev/null
+++ b/libgfortran/m4/matmul.m4
@@ -0,0 +1,145 @@
+`/* Implementation of the MATMUL intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_code, regexp(file, `_\([irc][0-9]+\)\.', `\1'))dnl
+define(rtype_letter,substr(rtype_code, 0, 1))dnl
+define(rtype_kind, substr(rtype_code, 1))dnl
+define(rtype,get_arraytype(rtype_letter,rtype_kind))dnl
+define(rtype_name, get_typename(rtype_letter, rtype_kind))dnl
+
+/* Dimensions: retarray(x,y) a(x, count) b(count,y).
+   Either a or b can be rank 1.  In this case x or y is 1.  */
+void
+`__matmul_'rtype_code (rtype * retarray, rtype * a, rtype * b)
+{
+  rtype_name *abase;
+  rtype_name *bbase;
+  rtype_name *dest;
+  rtype_name res;
+  index_type rxstride;
+  index_type rystride;
+  index_type xcount;
+  index_type ycount;
+  index_type xstride;
+  index_type ystride;
+  index_type x;
+  index_type y;
+
+  rtype_name *pa;
+  rtype_name *pb;
+  index_type astride;
+  index_type bstride;
+  index_type count;
+  index_type n;
+
+  assert (GFC_DESCRIPTOR_RANK (a) == 2
+          || GFC_DESCRIPTOR_RANK (b) == 2);
+  abase = a->data;
+  bbase = b->data;
+  dest = retarray->data;
+
+  if (retarray->dim[0].stride == 0)
+    retarray->dim[0].stride = 1;
+  if (a->dim[0].stride == 0)
+    a->dim[0].stride = 1;
+  if (b->dim[0].stride == 0)
+    b->dim[0].stride = 1;
+
+sinclude(`matmul_asm_'rtype_code`.m4')dnl
+
+  if (GFC_DESCRIPTOR_RANK (retarray) == 1)
+    {
+      rxstride = retarray->dim[0].stride;
+      rystride = rxstride;
+    }
+  else
+    {
+      rxstride = retarray->dim[0].stride;
+      rystride = retarray->dim[1].stride;
+    }
+
+  /* If we have rank 1 parameters, zero the absent stride, and set the size to
+     one.  */
+  if (GFC_DESCRIPTOR_RANK (a) == 1)
+    {
+      astride = a->dim[0].stride;
+      count = a->dim[0].ubound + 1 - a->dim[0].lbound;
+      xstride = 0;
+      rxstride = 0;
+      xcount = 1;
+    }
+  else
+    {
+      astride = a->dim[1].stride;
+      count = a->dim[1].ubound + 1 - a->dim[1].lbound;
+      xstride = a->dim[0].stride;
+      xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
+    }
+  if (GFC_DESCRIPTOR_RANK (b) == 1)
+    {
+      bstride = b->dim[0].stride;
+      assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
+      ystride = 0;
+      rystride = 0;
+      ycount = 1;
+    }
+  else
+    {
+      bstride = b->dim[0].stride;
+      assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
+      ystride = b->dim[1].stride;
+      ycount = b->dim[1].ubound + 1 - b->dim[1].lbound;
+    }
+
+  for (y = 0; y < ycount; y++)
+    {
+      for (x = 0; x < xcount; x++)
+        {
+          /* Do the summation for this element.  For real and integer types
+             this is the same as DOT_PRODUCT.  For complex types we use do
+             a*b, not conjg(a)*b.  */
+          pa = abase;
+          pb = bbase;
+          res = 0;
+
+          for (n = 0; n < count; n++)
+            {
+              res += *pa * *pb;
+              pa += astride;
+              pb += bstride;
+            }
+
+          *dest = res;
+
+          dest += rxstride;
+          abase += xstride;
+        }
+      abase -= xstride * xcount;
+      bbase += ystride;
+      dest += rystride - (rxstride * xcount);
+    }
+}
+
diff --git a/libgfortran/m4/matmull.m4 b/libgfortran/m4/matmull.m4
new file mode 100644
index 00000000000..e522a93faef
--- /dev/null
+++ b/libgfortran/m4/matmull.m4
@@ -0,0 +1,157 @@
+`/* Implementation of the MATMUL intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_kind, regexp(file, `_l\([0-9]+\)\.', `\1'))dnl
+define(rtype_code,`l'rtype_kind)dnl
+define(rtype,get_arraytype(l,rtype_kind))dnl
+define(rtype_name, get_typename(l, rtype_kind))dnl
+
+/* Dimensions: retarray(x,y) a(x, count) b(count,y).
+   Either a or b can be rank 1.  In this case x or y is 1.  */
+void
+`__matmul_'rtype_code (rtype * retarray, gfc_array_l4 * a, gfc_array_l4 * b)
+{
+  GFC_INTEGER_4 *abase;
+  GFC_INTEGER_4 *bbase;
+  rtype_name *dest;
+  index_type rxstride;
+  index_type rystride;
+  index_type xcount;
+  index_type ycount;
+  index_type xstride;
+  index_type ystride;
+  index_type x;
+  index_type y;
+
+  GFC_INTEGER_4 *pa;
+  GFC_INTEGER_4 *pb;
+  index_type astride;
+  index_type bstride;
+  index_type count;
+  index_type n;
+
+  assert (GFC_DESCRIPTOR_RANK (a) == 2
+          || GFC_DESCRIPTOR_RANK (b) == 2);
+  abase = a->data;
+  if (GFC_DESCRIPTOR_SIZE (a) != 4)
+    {
+      assert (GFC_DESCRIPTOR_SIZE (a) == 8);
+      abase = GFOR_POINTER_L8_TO_L4 (abase);
+      astride <<= 1;
+    }
+  bbase = b->data;
+  if (GFC_DESCRIPTOR_SIZE (b) != 4)
+    {
+      assert (GFC_DESCRIPTOR_SIZE (b) == 8);
+      bbase = GFOR_POINTER_L8_TO_L4 (bbase);
+      bstride <<= 1;
+    }
+  dest = retarray->data;
+
+  if (retarray->dim[0].stride == 0)
+    retarray->dim[0].stride = 1;
+  if (a->dim[0].stride == 0)
+    a->dim[0].stride = 1;
+  if (b->dim[0].stride == 0)
+    b->dim[0].stride = 1;
+
+sinclude(`matmul_asm_'rtype_code`.m4')dnl
+
+  if (GFC_DESCRIPTOR_RANK (retarray) == 1)
+    {
+      rxstride = retarray->dim[0].stride;
+      rystride = rxstride;
+    }
+  else
+    {
+      rxstride = retarray->dim[0].stride;
+      rystride = retarray->dim[1].stride;
+    }
+
+  /* If we have rank 1 parameters, zero the absent stride, and set the size to
+     one.  */
+  if (GFC_DESCRIPTOR_RANK (a) == 1)
+    {
+      astride = a->dim[0].stride;
+      count = a->dim[0].ubound + 1 - a->dim[0].lbound;
+      xstride = 0;
+      rxstride = 0;
+      xcount = 1;
+    }
+  else
+    {
+      astride = a->dim[1].stride;
+      count = a->dim[1].ubound + 1 - a->dim[1].lbound;
+      xstride = a->dim[0].stride;
+      xcount = a->dim[0].ubound + 1 - a->dim[0].lbound;
+    }
+  if (GFC_DESCRIPTOR_RANK (b) == 1)
+    {
+      bstride = b->dim[0].stride;
+      assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
+      ystride = 0;
+      rystride = 0;
+      ycount = 1;
+    }
+  else
+    {
+      bstride = b->dim[0].stride;
+      assert(count == b->dim[0].ubound + 1 - b->dim[0].lbound);
+      ystride = b->dim[1].stride;
+      ycount = b->dim[1].ubound + 1 - b->dim[1].lbound;
+    }
+
+  for (y = 0; y < ycount; y++)
+    {
+      for (x = 0; x < xcount; x++)
+        {
+          /* Do the summation for this element.  For real and integer types
+             this is the same as DOT_PRODUCT.  For complex types we use do
+             a*b, not conjg(a)*b.  */
+          pa = abase;
+          pb = bbase;
+          *dest = 0;
+
+          for (n = 0; n < count; n++)
+            {
+              if (*pa && *pb)
+                {
+                  *dest = 1;
+                  break;
+                }
+              pa += astride;
+              pb += bstride;
+            }
+
+          dest += rxstride;
+          abase += xstride;
+        }
+      abase -= xstride * xcount;
+      bbase += ystride;
+      dest += rystride - (rxstride * xcount);
+    }
+}
+
diff --git a/libgfortran/m4/maxloc0.m4 b/libgfortran/m4/maxloc0.m4
new file mode 100644
index 00000000000..89ecacb9656
--- /dev/null
+++ b/libgfortran/m4/maxloc0.m4
@@ -0,0 +1,54 @@
+`/* Implementation of the MAXLOC intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <float.h>
+#include <limits.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(iforeach.m4)dnl
+
+FOREACH_FUNCTION(
+`  type_name maxval;
+
+  maxval = type_min;'
+,
+`  if (*base > maxval)
+    {
+      maxval = *base;
+      for (n = 0; n < rank; n++)
+        dest[n * dstride] = count[n] + 1;
+    }')
+
+MASKED_FOREACH_FUNCTION(
+`  type_name maxval;
+
+  maxval = type_min;'
+,
+`  if (*mbase && *base > maxval)
+    {
+      maxval = *base;
+      for (n = 0; n < rank; n++)
+        dest[n * dstride] = count[n] + 1;
+    }')
diff --git a/libgfortran/m4/maxloc1.m4 b/libgfortran/m4/maxloc1.m4
new file mode 100644
index 00000000000..0eb259f9c6d
--- /dev/null
+++ b/libgfortran/m4/maxloc1.m4
@@ -0,0 +1,50 @@
+`/* Implementation of the MAXLOC intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <float.h>
+#include <limits.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+ARRAY_FUNCTION(0,
+`  type_name maxval;
+  maxval = type_min;
+  result = 1;',
+`  if (*src > maxval)
+    {
+      maxval = *src;
+      result = (rtype_name)n + 1;
+    }')
+
+MASKED_ARRAY_FUNCTION(0,
+`  type_name maxval;
+  maxval = type_min;
+  result = 1;',
+`  if (*msrc && *src > maxval)
+    {
+      maxval = *src;
+      result = (rtype_name)n + 1;
+    }')
+
diff --git a/libgfortran/m4/maxval.m4 b/libgfortran/m4/maxval.m4
new file mode 100644
index 00000000000..b6a5666aae1
--- /dev/null
+++ b/libgfortran/m4/maxval.m4
@@ -0,0 +1,39 @@
+`/* Implementation of the MAXVAL intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <float.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+ARRAY_FUNCTION(type_min,
+`  result = type_min;',
+`  if (*src > result)
+    result = *src;')
+
+MASKED_ARRAY_FUNCTION(type_min,
+`  result = type_min;',
+`  if (*msrc && *src > result)
+    result = *src;')
+
diff --git a/libgfortran/m4/minloc0.m4 b/libgfortran/m4/minloc0.m4
new file mode 100644
index 00000000000..5411418b68a
--- /dev/null
+++ b/libgfortran/m4/minloc0.m4
@@ -0,0 +1,54 @@
+`/* Implementation of the MINLOC intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <float.h>
+#include <limits.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(iforeach.m4)dnl
+
+FOREACH_FUNCTION(
+`  type_name minval;
+
+  minval = type_max;'
+,
+`  if (*base < minval)
+    {
+      minval = *base;
+      for (n = 0; n < rank; n++)
+        dest[n * dstride] = count[n] + 1;
+    }')
+
+MASKED_FOREACH_FUNCTION(
+`  type_name minval;
+
+  minval = type_max;'
+,
+`  if (*mbase && *base < minval)
+    {
+      minval = *base;
+      for (n = 0; n < rank; n++)
+        dest[n * dstride] = count[n] + 1;
+    }')
diff --git a/libgfortran/m4/minloc1.m4 b/libgfortran/m4/minloc1.m4
new file mode 100644
index 00000000000..e3101c63861
--- /dev/null
+++ b/libgfortran/m4/minloc1.m4
@@ -0,0 +1,50 @@
+`/* Implementation of the MINLOC intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <float.h>
+#include <limits.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+ARRAY_FUNCTION(0,
+`  type_name minval;
+  minval = type_max;
+  result = 1;',
+`  if (*src < minval)
+    {
+      minval = *src;
+      result = (rtype_name)n + 1;
+    }')
+
+MASKED_ARRAY_FUNCTION(0,
+`  type_name minval;
+  minval = type_max;
+  result = 1;',
+`  if (*msrc && *src < minval)
+    {
+      minval = *src;
+      result = (rtype_name)n + 1;
+    }')
+
diff --git a/libgfortran/m4/minval.m4 b/libgfortran/m4/minval.m4
new file mode 100644
index 00000000000..2c1be2b88f9
--- /dev/null
+++ b/libgfortran/m4/minval.m4
@@ -0,0 +1,39 @@
+`/* Implementation of the MINVAL intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <float.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+ARRAY_FUNCTION(type_max,
+`  result = type_max;',
+`  if (*src < result)
+    result = *src;')
+
+MASKED_ARRAY_FUNCTION(type_max,
+`  result = type_max;',
+`  if (*msrc && *src < result)
+    result = *src;')
+
diff --git a/libgfortran/m4/mtype.m4 b/libgfortran/m4/mtype.m4
new file mode 100644
index 00000000000..84bf39f3561
--- /dev/null
+++ b/libgfortran/m4/mtype.m4
@@ -0,0 +1,5 @@
+dnl Get type kind from filename.
+define(kind,regexp(file, `_.\([0-9]+\).c$', `\1'))dnl
+define(complex_type, `GFC_COMPLEX_'kind)dnl
+define(real_type, `GFC_REAL_'kind)dnl
+define(q,ifelse(kind,4,f,ifelse(kind,8,`',`_'kind)))dnl
diff --git a/libgfortran/m4/nearest.m4 b/libgfortran/m4/nearest.m4
new file mode 100644
index 00000000000..5168d99c8e2
--- /dev/null
+++ b/libgfortran/m4/nearest.m4
@@ -0,0 +1,39 @@
+`/* Implementation of the NEAREST intrinsic
+   Copyright 2003 Free Software Foundation, Inc.
+   Contributed by Richard Henderson <rth@redhat.com>.
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+#include <math.h>
+#include <float.h>
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+real_type
+prefix(nearest_r`'kind) (real_type s, real_type dir)
+{
+  dir = copysign`'q (__builtin_inf`'q (), dir);
+  if (FLT_EVAL_METHOD != 0)
+    {
+      /* ??? Work around glibc bug on x86.  */
+      volatile real_type r = nextafter`'q (s, dir);
+      return r;
+    }
+  else
+    return nextafter`'q (s, dir);
+}
diff --git a/libgfortran/m4/product.m4 b/libgfortran/m4/product.m4
new file mode 100644
index 00000000000..bef9b8184eb
--- /dev/null
+++ b/libgfortran/m4/product.m4
@@ -0,0 +1,37 @@
+`/* Implementation of the PRODUCT intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+ARRAY_FUNCTION(1,
+`  result = 1;',
+`  result *= *src;')
+
+MASKED_ARRAY_FUNCTION(1,
+`  result = 1;',
+`  if (*msrc)
+    result *= *src;')
+
diff --git a/libgfortran/m4/reshape.m4 b/libgfortran/m4/reshape.m4
new file mode 100644
index 00000000000..b8143fe24af
--- /dev/null
+++ b/libgfortran/m4/reshape.m4
@@ -0,0 +1,232 @@
+`/* Implementation of the RESHAPE
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_kind, regexp(file, `_.\([0-9]+\)\.', `\1'))dnl
+define(rtype_letter, regexp(file, `_\(.\)[0-9]+\.', `\1'))dnl
+define(rtype_code,rtype_letter`'rtype_name)dnl
+define(rtype,get_arraytype(rtype_letter,rtype_kind))dnl
+define(rtype_name, get_typename(rtype_letter, rtype_kind))dnl
+
+typedef GFC_ARRAY_DESCRIPTOR(1, index_type) shape_type;
+
+/* The shape parameter is ignored. We can currently deduce the shape from the
+   return array.  */
+dnl Only the kind (ie size) is used to name the function.
+void
+`__reshape_'rtype_kind (rtype * ret, rtype * source, shape_type * shape,
+                      rtype * pad, shape_type * order)
+{
+  /* r.* indicates the return array.  */
+  index_type rcount[GFC_MAX_DIMENSIONS - 1];
+  index_type rextent[GFC_MAX_DIMENSIONS - 1];
+  index_type rstride[GFC_MAX_DIMENSIONS - 1];
+  index_type rstride0;
+  index_type rdim;
+  index_type rsize;
+  rtype_name *rptr;
+  /* s.* indicates the source array.  */
+  index_type scount[GFC_MAX_DIMENSIONS - 1];
+  index_type sextent[GFC_MAX_DIMENSIONS - 1];
+  index_type sstride[GFC_MAX_DIMENSIONS - 1];
+  index_type sstride0;
+  index_type sdim;
+  index_type ssize;
+  const rtype_name *sptr;
+  /* p.* indicates the pad array.  */
+  index_type pcount[GFC_MAX_DIMENSIONS - 1];
+  index_type pextent[GFC_MAX_DIMENSIONS - 1];
+  index_type pstride[GFC_MAX_DIMENSIONS - 1];
+  index_type pdim;
+  index_type psize;
+  const rtype_name *pptr;
+
+  const rtype_name *src;
+  int n;
+  int dim;
+
+  if (ret->dim[0].stride == 0)
+    ret->dim[0].stride = 1;
+  if (source->dim[0].stride == 0)
+    source->dim[0].stride = 1;
+  if (shape->dim[0].stride == 0)
+    shape->dim[0].stride = 1;
+  if (pad && pad->dim[0].stride == 0)
+    pad->dim[0].stride = 1;
+  if (order && order->dim[0].stride == 0)
+    order->dim[0].stride = 1;
+
+  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[dim] <= 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)
+    {
+      if (pad->dim[0].stride == 0)
+        pad->dim[0].stride = 1;
+      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 *= rtype_kind;
+      ssize *= rtype_kind;
+      psize *= rtype_kind;
+      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];
+            }
+        }
+    }
+}
+
diff --git a/libgfortran/m4/set_exponent.m4 b/libgfortran/m4/set_exponent.m4
new file mode 100644
index 00000000000..352a129865e
--- /dev/null
+++ b/libgfortran/m4/set_exponent.m4
@@ -0,0 +1,31 @@
+`/* Implementation of the SET_EXPONENT intrinsic
+   Copyright 2003 Free Software Foundation, Inc.
+   Contributed by Richard Henderson <rth@redhat.com>.
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+#include <math.h>
+#include "libgfortran.h"'
+
+include(`mtype.m4')dnl
+
+real_type
+prefix(set_exponent_r`'kind) (real_type s, GFC_INTEGER_4 i)
+{
+  int dummy_exp;
+  return scalbn`'q (frexp`'q (s, &dummy_exp), i);
+}
diff --git a/libgfortran/m4/shape.m4 b/libgfortran/m4/shape.m4
new file mode 100644
index 00000000000..826077eabe7
--- /dev/null
+++ b/libgfortran/m4/shape.m4
@@ -0,0 +1,48 @@
+`/* Implementation of the SHAPE intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_kind, regexp(file, `_i\([0-9]+\)\.', `\1'))dnl
+define(rtype_code,`i'rtype_name)dnl
+define(rtype,get_arraytype(i,rtype_kind))dnl
+define(rtype_name, get_typename(i, rtype_kind))dnl
+
+void
+`__shape_'rtype_kind (rtype * ret, const rtype * array)
+{
+  int n;
+  index_type stride;
+
+  stride = ret->dim[0].stride;
+  if (stride == 0)
+    stride = 1;
+
+  for (n = 0; n < GFC_DESCRIPTOR_RANK (array); n++)
+    {
+      ret->data[n * stride] =
+        array->dim[n].ubound + 1 - array->dim[n].lbound;
+    }
+}
+
diff --git a/libgfortran/m4/specific.m4 b/libgfortran/m4/specific.m4
new file mode 100644
index 00000000000..bf88857ccf1
--- /dev/null
+++ b/libgfortran/m4/specific.m4
@@ -0,0 +1,16 @@
+include(head.m4)
+define(type_code,regexp(file,`_\([ircl][0-9]+\).f90',`\1'))dnl
+define(type_letter,substr(type_code, 0, 1))dnl
+define(type_kind,substr(type_code, 1))dnl
+define(get_typename2, `$1 (kind=$2)')dnl
+define(get_typename, `get_typename2(ifelse($1,i,integer,ifelse($1,r,real,ifelse($1,l,logical,ifelse($1,c,complex,unknown)))),`$2')')dnl
+define(type_name, get_typename(type_letter,type_kind))dnl
+define(name, regexp(regexp(file, `[^/]*$', `\&'), `^_\([^_]*\)_', `\1'))dnl
+define(function_name,`specific__'name`_'type_code)dnl
+
+elemental function function_name (parm)
+   type_name, intent (in) :: parm
+   type_name :: function_name
+
+   function_name = name (parm)
+end function
diff --git a/libgfortran/m4/specific2.m4 b/libgfortran/m4/specific2.m4
new file mode 100644
index 00000000000..ca0d831f5e0
--- /dev/null
+++ b/libgfortran/m4/specific2.m4
@@ -0,0 +1,16 @@
+include(head.m4)
+define(type_code,regexp(file,`_\([ircl][0-9]+\).f90',`\1'))dnl
+define(type_letter,substr(type_code, 0, 1))dnl
+define(type_kind,substr(type_code, 1))dnl
+define(get_typename2, `$1 (kind=$2)')dnl
+define(get_typename, `get_typename2(ifelse($1,i,integer,ifelse($1,r,real,ifelse($1,l,logical,ifelse($1,c,complex,unknown)))),`$2')')dnl
+define(type_name, get_typename(type_letter,type_kind))dnl
+define(name, regexp(regexp(file, `[^/]*$', `\&'), `^_\([^_]*\)_', `\1'))dnl
+define(function_name,`specific__'name`_'type_code)dnl
+
+elemental function function_name (p1, p2)
+   type_name, intent (in) :: p1, p2
+   type_name :: function_name
+
+   function_name = name (p1, p2)
+end function
diff --git a/libgfortran/m4/sum.m4 b/libgfortran/m4/sum.m4
new file mode 100644
index 00000000000..0ea3477bf87
--- /dev/null
+++ b/libgfortran/m4/sum.m4
@@ -0,0 +1,36 @@
+`/* Implementation of the SUM intrinsic
+   Copyright 2002 Free Software Foundation, Inc.
+   Contributed by Paul Brook <paul@nowt.org>
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfortran is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Libgfortran is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <stdlib.h>
+#include <assert.h>
+#include "libgfortran.h"'
+
+include(iparm.m4)dnl
+include(ifunction.m4)dnl
+ARRAY_FUNCTION(0,
+`  result = 0;',
+`  result += *src;')
+
+MASKED_ARRAY_FUNCTION(0,
+`  result = 0;',
+`  if (*msrc)
+    result += *src;')
diff --git a/libgfortran/m4/transpose.m4 b/libgfortran/m4/transpose.m4
new file mode 100644
index 00000000000..35df64bf347
--- /dev/null
+++ b/libgfortran/m4/transpose.m4
@@ -0,0 +1,75 @@
+`/* Implementation of the TRANSPOSE intrinsic
+   Copyright 2003 Free Software Foundation, Inc.
+   Contributed by Tobias Schlüter
+
+This file is part of the GNU Fortran 95 runtime library (libgfor).
+
+Libgfor is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+Ligbfor 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with libgfor; see the file COPYING.LIB.  If not,
+write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+#include "config.h"
+#include <assert.h>
+#include "libgfortran.h"'
+include(types.m4)dnl
+define(rtype_kind, regexp(file, `_.\([0-9]+\)\.', `\1'))dnl
+define(rtype_letter, regexp(file, `_\(.\)[0-9]+\.', `\1'))dnl
+define(rtype_code,rtype_letter`'rtype_name)dnl
+define(rtype,get_arraytype(rtype_letter,rtype_kind))dnl
+define(rtype_name, get_typename(rtype_letter, rtype_kind))dnl
+
+void
+`__transpose_'rtype_kind (rtype * ret, rtype * source)
+{
+  /* r.* indicates the return array.  */
+  index_type rxstride, rystride;
+  rtype_name *rptr;
+  /* s.* indicates the source array.  */
+  index_type sxstride, systride;
+  const rtype_name *sptr;
+
+  index_type xcount, ycount;
+  index_type x, y;
+
+  assert (GFC_DESCRIPTOR_RANK (source) == 2);
+
+  if (ret->dim[0].stride == 0)
+    ret->dim[0].stride = 1;
+  if (source->dim[0].stride == 0)
+    source->dim[0].stride = 1;
+
+  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);
+    }
+}
diff --git a/libgfortran/m4/types.m4 b/libgfortran/m4/types.m4
new file mode 100644
index 00000000000..cb808290c08
--- /dev/null
+++ b/libgfortran/m4/types.m4
@@ -0,0 +1,4 @@
+define(get_typename2, `GFC_$1_$2')dnl
+define(get_typename, `get_typename2(ifelse($1,i,INTEGER,ifelse($1,r,REAL,ifelse($1,l,LOGICAL,ifelse($1,c,COMPLEX,unknown)))),`$2')')dnl
+define(get_arraytype, `gfc_array_$1$2')dnl
+define(name, regexp(regexp(file, `[^/]*$', `\&'), `^\([^_]*\)_', `\1'))dnl