[PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

[-- Attachment #1: Type: text/plain, Size: 1520 bytes --]

The attached patch provides a partial implementation for
the simplification for CSHIFT.  It is partial in that it
only applies to rank 1 arrays.  For arrays with rank > 1,
gfc_simplify_cshift will issue an error.  Here, the intent
is that hopefully someone that knows what they are doing
with supply a patch for rank > 1.

The meat of the patch for rank = 1 may not be the most
efficient.  It copies the array elements from 'a' to 'result'
in the circularly shifted order.  It inefficiently always
starts with the first element in 'a' to find the candidate
element for next 'result' element.

      cr = gfc_constructor_first (result->value.constructor);
      for (i = 0; i < sz; i++, cr = gfc_constructor_next (cr))
	{
	  j = (i + shft) % sz;
	  ca = gfc_constructor_first (a->value.constructor);
	  while (j-- > 0)
	    ca = gfc_constructor_next (ca);
	  cr->expr = gfc_copy_expr (ca->expr);
	}

As the values are storied in a splay tree, there may be
a more efficient way to split the splay and recombine
it into a new.

Anyway, I would like to commit the attached patch.
Built and tested on x86_64-*-freebsd?

2015-11-19  Steven G. Kargl  <kargl@gcc.gnu.org>

	* intrinsic.h: Prototype for gfc_simplify_cshift
	* intrinsic.c (add_functions): Use gfc_simplify_cshift.
	* simplify.c (gfc_simplify_cshift): Implement simplification of CSHIFT.
	(gfc_simplify_spread): Remove a FIXME and add error condition.
 
2015-11-19  Steven G. Kargl  <kargl@gcc.gnu.org>

	* gfortran.dg/simplify_cshift_1.f90: New test.

-- 
Steve

[-- Attachment #2: cshift.diff --]
[-- Type: text/x-diff, Size: 5109 bytes --]

Index: gcc/fortran/intrinsic.c
===================================================================
--- gcc/fortran/intrinsic.c	(revision 230585)
+++ gcc/fortran/intrinsic.c	(working copy)
@@ -1659,9 +1659,11 @@ add_functions (void)
 
   make_generic ("count", GFC_ISYM_COUNT, GFC_STD_F95);
 
-  add_sym_3 ("cshift", GFC_ISYM_CSHIFT, CLASS_TRANSFORMATIONAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_F95,
-	     gfc_check_cshift, NULL, gfc_resolve_cshift,
-	     ar, BT_REAL, dr, REQUIRED, sh, BT_INTEGER, di, REQUIRED,
+  add_sym_3 ("cshift", GFC_ISYM_CSHIFT, CLASS_TRANSFORMATIONAL, ACTUAL_NO,
+	     BT_REAL, dr, GFC_STD_F95,
+	     gfc_check_cshift, gfc_simplify_cshift, gfc_resolve_cshift,
+	     ar, BT_REAL, dr, REQUIRED,
+	     sh, BT_INTEGER, di, REQUIRED,
 	     dm, BT_INTEGER, ii, OPTIONAL);
 
   make_generic ("cshift", GFC_ISYM_CSHIFT, GFC_STD_F95);
Index: gcc/fortran/intrinsic.h
===================================================================
--- gcc/fortran/intrinsic.h	(revision 230585)
+++ gcc/fortran/intrinsic.h	(working copy)
@@ -271,6 +271,7 @@ gfc_expr *gfc_simplify_conjg (gfc_expr *
 gfc_expr *gfc_simplify_cos (gfc_expr *);
 gfc_expr *gfc_simplify_cosh (gfc_expr *);
 gfc_expr *gfc_simplify_count (gfc_expr *, gfc_expr *, gfc_expr *);
+gfc_expr *gfc_simplify_cshift (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_dcmplx (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_dble (gfc_expr *);
 gfc_expr *gfc_simplify_digits (gfc_expr *);
Index: gcc/fortran/simplify.c
===================================================================
--- gcc/fortran/simplify.c	(revision 230585)
+++ gcc/fortran/simplify.c	(working copy)
@@ -1789,6 +1789,88 @@ gfc_simplify_count (gfc_expr *mask, gfc_
 
 
 gfc_expr *
+gfc_simplify_cshift (gfc_expr *array, gfc_expr *shift, gfc_expr *dim)
+{
+  gfc_expr *a;
+
+  a = gfc_copy_expr (array);
+
+  switch (a->expr_type)
+    {
+      case EXPR_VARIABLE:
+      case EXPR_ARRAY:
+	gfc_simplify_expr (a, 0);
+	if (!is_constant_array_expr (a))
+	  {
+	    gfc_free_expr (a);
+	    return NULL;
+	  }
+	break;
+      default:
+	gcc_unreachable ();
+    }
+
+  if (a->rank == 1)
+    {
+      gfc_constructor *ca, *cr;
+      gfc_expr *result;
+      mpz_t size;
+      int i, j, shft, sz;
+
+      if (!gfc_is_constant_expr (shift))
+	return NULL;
+
+      shft = mpz_get_si (shift->value.integer);
+
+      /* Special case: rank 1 array with no shift!  */
+      if (shft == 0)
+	return a;
+
+      /*  Case (i):  If ARRAY has rank one, element i of the result is
+	  ARRAY (1 + MODULO (i + SHIFT ­ 1, SIZE (ARRAY))).  */
+
+      result = gfc_copy_expr (a);
+      mpz_init (size);
+      gfc_array_size (a, &size);
+      sz = mpz_get_si (size);
+      mpz_clear (size);
+      shft = shft < 0 ? 1 - shft : shft;
+      cr = gfc_constructor_first (result->value.constructor);
+      for (i = 0; i < sz; i++, cr = gfc_constructor_next (cr))
+	{
+	  j = (i + shft) % sz;
+	  ca = gfc_constructor_first (a->value.constructor);
+	  while (j-- > 0)
+	    ca = gfc_constructor_next (ca);
+	  cr->expr = gfc_copy_expr (ca->expr);
+	}
+
+      gfc_free_expr (a);
+      return result;
+    }
+  else
+    {
+      int dm;
+
+      if (dim)
+	{
+	  if (!gfc_is_constant_expr (dim))
+	    return NULL;
+
+	  dm = mpz_get_si (dim->value.integer);
+	}
+      else
+	dm = 1;
+
+      gfc_error ("Simplification of CSHIFT with an array with rank > 1 "
+	         "no yet support");
+    }
+
+  return NULL;
+}
+
+
+gfc_expr *
 gfc_simplify_dcmplx (gfc_expr *x, gfc_expr *y)
 {
   return simplify_cmplx ("DCMPLX", x, y, gfc_default_double_kind);
@@ -6089,10 +6171,11 @@ gfc_simplify_spread (gfc_expr *source, g
 	}
     }
   else
-    /* FIXME: Returning here avoids a regression in array_simplify_1.f90.
-       Replace NULL with gcc_unreachable() after implementing
-       gfc_simplify_cshift().  */
-    return NULL;
+    {
+      gfc_error ("Simplification of SPREAD at %L not yet implemented",
+		 &source->where);
+      return &gfc_bad_expr;
+    }
 
   if (source->ts.type == BT_CHARACTER)
     result->ts.u.cl = source->ts.u.cl;
Index: gcc/testsuite/gfortran.dg/simplify_cshift_1.f90
===================================================================
--- gcc/testsuite/gfortran.dg/simplify_cshift_1.f90	(nonexistent)
+++ gcc/testsuite/gfortran.dg/simplify_cshift_1.f90	(working copy)
@@ -0,0 +1,38 @@
+! { dg-do compile }
+program foo
+
+   implicit none
+   
+   type t
+      integer i
+   end type t
+
+   type(t), parameter :: d(5) = [t(1), t(2), t(3), t(4), t(5)]
+   type(t) e(5), q(5)
+
+   integer, parameter :: a(5) = [1, 2, 3, 4, 5]
+   integer i, b(5), c(5), v(5)
+
+   c = [1, 2, 3, 4, 5]
+
+   b = cshift(a, -2)
+   v = cshift(c, -2)
+   if (any(b /= v)) call abort
+
+   b = cshift(a, 2)
+   v = cshift(c,2)
+   if (any(b /= v)) call abort
+
+   b = cshift([1, 2, 3, 4, 5], 0)
+   if (any(b /= a)) call abort
+   b = cshift(2*a, 0)
+   if (any(b /= 2*a)) call abort
+ 
+   e = [t(1), t(2), t(3), t(4), t(5)]
+   e = cshift(e, 3)
+   q = cshift(d, 3)
+   do i = 1, 5
+      if (e(i)%i /= q(i)%i) call abort
+   end do
+
+end program foo

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

On Thu, Nov 19, 2015 at 04:58:36PM -0800, Steve Kargl wrote:
> +  else
> +    {
> +      int dm;
> +
> +      if (dim)
> +	{
> +	  if (!gfc_is_constant_expr (dim))
> +	    return NULL;
> +
> +	  dm = mpz_get_si (dim->value.integer);
> +	}
> +      else
> +	dm = 1;
> +
> +      gfc_error ("Simplification of CSHIFT with an array with rank > 1 "
> +	         "no yet support");
> +    }
> +

To save some time, the dim portion of the patch isn't
correct.  dim can be scalar or rank 1 array.  I'll
#if 0 ... #endif this section unless I persevere with
the rank > 1 case.

-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

On Thu, Nov 19, 2015 at 05:31:32PM -0800, Steve Kargl wrote:
> On Thu, Nov 19, 2015 at 04:58:36PM -0800, Steve Kargl wrote:
> > +  else
> > +    {
> > +      int dm;
> > +
> > +      if (dim)
> > +	{
> > +	  if (!gfc_is_constant_expr (dim))
> > +	    return NULL;
> > +
> > +	  dm = mpz_get_si (dim->value.integer);
> > +	}
> > +      else
> > +	dm = 1;
> > +
> > +      gfc_error ("Simplification of CSHIFT with an array with rank > 1 "
> > +	         "no yet support");
> > +    }
> > +
> 
> To save some time, the dim portion of the patch isn't
> correct.  dim can be scalar or rank 1 array.  I'll
> #if 0 ... #endif this section unless I persevere with
> the rank > 1 case.

Ugh.  Too much gdb today.  The above is correct.  I conflated
SHIFT and DIM's requirements.

-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

[-- Attachment #1: Type: text/plain, Size: 728 bytes --]

On Thu, Nov 19, 2015 at 04:58:36PM -0800, Steve Kargl wrote:
> 
> 2015-11-19  Steven G. Kargl  <kargl@gcc.gnu.org>
> 
> 	* intrinsic.h: Prototype for gfc_simplify_cshift
> 	* intrinsic.c (add_functions): Use gfc_simplify_cshift.
> 	* simplify.c (gfc_simplify_cshift): Implement simplification of CSHIFT.
> 	(gfc_simplify_spread): Remove a FIXME and add error condition.
>  
> 2015-11-19  Steven G. Kargl  <kargl@gcc.gnu.org>
> 
> 	* gfortran.dg/simplify_cshift_1.f90: New test.
> 

I've attached an updated patch.  The changes consists of
1) better/more comments
2) re-organize code to reduce copying of the array.
3) add optimization for a left/right shift that 
   returns the original array.
4) Don't leak memory.

-- 
Steve

[-- Attachment #2: cshift.diff --]
[-- Type: text/x-diff, Size: 5926 bytes --]

Index: gcc/fortran/intrinsic.c
===================================================================
--- gcc/fortran/intrinsic.c	(revision 230585)
+++ gcc/fortran/intrinsic.c	(working copy)
@@ -1659,9 +1659,11 @@ add_functions (void)
 
   make_generic ("count", GFC_ISYM_COUNT, GFC_STD_F95);
 
-  add_sym_3 ("cshift", GFC_ISYM_CSHIFT, CLASS_TRANSFORMATIONAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_F95,
-	     gfc_check_cshift, NULL, gfc_resolve_cshift,
-	     ar, BT_REAL, dr, REQUIRED, sh, BT_INTEGER, di, REQUIRED,
+  add_sym_3 ("cshift", GFC_ISYM_CSHIFT, CLASS_TRANSFORMATIONAL, ACTUAL_NO,
+	     BT_REAL, dr, GFC_STD_F95,
+	     gfc_check_cshift, gfc_simplify_cshift, gfc_resolve_cshift,
+	     ar, BT_REAL, dr, REQUIRED,
+	     sh, BT_INTEGER, di, REQUIRED,
 	     dm, BT_INTEGER, ii, OPTIONAL);
 
   make_generic ("cshift", GFC_ISYM_CSHIFT, GFC_STD_F95);
Index: gcc/fortran/intrinsic.h
===================================================================
--- gcc/fortran/intrinsic.h	(revision 230585)
+++ gcc/fortran/intrinsic.h	(working copy)
@@ -271,6 +271,7 @@ gfc_expr *gfc_simplify_conjg (gfc_expr *
 gfc_expr *gfc_simplify_cos (gfc_expr *);
 gfc_expr *gfc_simplify_cosh (gfc_expr *);
 gfc_expr *gfc_simplify_count (gfc_expr *, gfc_expr *, gfc_expr *);
+gfc_expr *gfc_simplify_cshift (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_dcmplx (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_dble (gfc_expr *);
 gfc_expr *gfc_simplify_digits (gfc_expr *);
Index: gcc/fortran/simplify.c
===================================================================
--- gcc/fortran/simplify.c	(revision 230585)
+++ gcc/fortran/simplify.c	(working copy)
@@ -1789,6 +1789,99 @@ gfc_simplify_count (gfc_expr *mask, gfc_
 
 
 gfc_expr *
+gfc_simplify_cshift (gfc_expr *array, gfc_expr *shift, gfc_expr *dim)
+{
+  int dm;
+  gfc_expr *a;
+
+  /* DIM is only useful for rank > 1, but deal with it here as one can
+     set DIM = 1 for rank = 1.  */
+  if (dim)
+    {
+      if (!gfc_is_constant_expr (dim))
+	  return NULL;
+      dm = mpz_get_si (dim->value.integer);
+    }
+  else
+    dm = 1;
+
+  /* Copy array into 'a', simplify it, and then test for a constant array.
+     Unexpected expr_type cause an ICE.   */
+  switch (array->expr_type)
+    {
+      case EXPR_VARIABLE:
+      case EXPR_ARRAY:
+	a = gfc_copy_expr (array);
+	gfc_simplify_expr (a, 0);
+	if (!is_constant_array_expr (a))
+	  {
+	    gfc_free_expr (a);
+	    return NULL;
+	  }
+	break;
+      default:
+	gcc_unreachable ();
+    }
+
+  if (a->rank == 1)
+    {
+      gfc_constructor *ca, *cr;
+      gfc_expr *result;
+      mpz_t size;
+      int i, j, shft, sz;
+
+      if (!gfc_is_constant_expr (shift))
+	{
+	  gfc_free_expr (a);
+	  return NULL;
+	}
+
+      shft = mpz_get_si (shift->value.integer);
+
+      /*  Case (i):  If ARRAY has rank one, element i of the result is
+	  ARRAY (1 + MODULO (i + SHIFT ­ 1, SIZE (ARRAY))).  */
+
+      mpz_init (size);
+      gfc_array_size (a, &size);
+      sz = mpz_get_si (size);
+      mpz_clear (size);
+
+      /* Special case: rank 1 array with no shift or a complete shift to
+	 the original order!  */
+      if (shft == 0 || shft == sz || shft == 1 - sz)
+	return a;
+
+      /* Adjust shft to deal with right or left shifts. */
+      shft = shft < 0 ? 1 - shft : shft;
+
+      result = gfc_copy_expr (a);
+      cr = gfc_constructor_first (result->value.constructor);
+      for (i = 0; i < sz; i++, cr = gfc_constructor_next (cr))
+	{
+	  j = (i + shft) % sz;
+	  ca = gfc_constructor_first (a->value.constructor);
+	  while (j-- > 0)
+	    ca = gfc_constructor_next (ca);
+	  cr->expr = gfc_copy_expr (ca->expr);
+	}
+
+      gfc_free_expr (a);
+      return result;
+    }
+  else
+    {
+      /* FIXME: Deal with rank > 1 arrays.  For now, don't leak memory
+	 and exit with an error message.  */
+      gfc_free_expr (a);
+      gfc_error ("Simplification of CSHIFT with an array with rank > 1 "
+	         "no yet support");
+    }
+
+  return NULL;
+}
+
+
+gfc_expr *
 gfc_simplify_dcmplx (gfc_expr *x, gfc_expr *y)
 {
   return simplify_cmplx ("DCMPLX", x, y, gfc_default_double_kind);
@@ -6089,10 +6182,11 @@ gfc_simplify_spread (gfc_expr *source, g
 	}
     }
   else
-    /* FIXME: Returning here avoids a regression in array_simplify_1.f90.
-       Replace NULL with gcc_unreachable() after implementing
-       gfc_simplify_cshift().  */
-    return NULL;
+    {
+      gfc_error ("Simplification of SPREAD at %L not yet implemented",
+		 &source->where);
+      return &gfc_bad_expr;
+    }
 
   if (source->ts.type == BT_CHARACTER)
     result->ts.u.cl = source->ts.u.cl;
Index: gcc/testsuite/gfortran.dg/simplify_cshift_1.f90
===================================================================
--- gcc/testsuite/gfortran.dg/simplify_cshift_1.f90	(nonexistent)
+++ gcc/testsuite/gfortran.dg/simplify_cshift_1.f90	(working copy)
@@ -0,0 +1,46 @@
+! { dg-do run }
+program foo
+
+   implicit none
+   
+   type t
+      integer i
+   end type t
+
+   type(t), parameter :: d(5) = [t(1), t(2), t(3), t(4), t(5)]
+   type(t) e(5), q(5)
+
+   integer, parameter :: a(5) = [1, 2, 3, 4, 5]
+   integer i, b(5), c(5), v(5)
+
+   c = [1, 2, 3, 4, 5]
+
+   b = cshift(a, -2)
+   v = cshift(c, -2)
+   if (any(b /= v)) call abort
+
+   b = cshift(a, 2)
+   v = cshift(c, 2)
+   if (any(b /= v)) call abort
+
+   ! Special cases shift = 0, size(a), 1-size(a)
+   b = cshift([1, 2, 3, 4, 5], 0)
+   if (any(b /= a)) call abort
+   b = cshift([1, 2, 3, 4, 5], size(a))
+   if (any(b /= a)) call abort
+   b = cshift([1, 2, 3, 4, 5], 1-size(a))
+   if (any(b /= a)) call abort
+
+   ! simplification of array arg.
+   b = cshift(2 * a, 0)
+   if (any(b /= 2 * a)) call abort
+
+   ! An array of derived types workd too.
+   e = [t(1), t(2), t(3), t(4), t(5)]
+   e = cshift(e, 3)
+   q = cshift(d, 3)
+   do i = 1, 5
+      if (e(i)%i /= q(i)%i) call abort
+   end do
+
+end program foo

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Paul Richard Thomas]

Hi Steve,

Just a couple of small typos:
"Unexpected expr_type cause an ICE" ;  causes?
"! An array of derived types workd too." ; works?

Apart from that it's OK for trunk.

Thanks for the patch

Cheers

Paul


On 20 November 2015 at 21:09, Steve Kargl
<sgk@troutmask.apl.washington.edu> wrote:
> On Thu, Nov 19, 2015 at 04:58:36PM -0800, Steve Kargl wrote:
>>
>> 2015-11-19  Steven G. Kargl  <kargl@gcc.gnu.org>
>>
>>       * intrinsic.h: Prototype for gfc_simplify_cshift
>>       * intrinsic.c (add_functions): Use gfc_simplify_cshift.
>>       * simplify.c (gfc_simplify_cshift): Implement simplification of CSHIFT.
>>       (gfc_simplify_spread): Remove a FIXME and add error condition.
>>
>> 2015-11-19  Steven G. Kargl  <kargl@gcc.gnu.org>
>>
>>       * gfortran.dg/simplify_cshift_1.f90: New test.
>>
>
> I've attached an updated patch.  The changes consists of
> 1) better/more comments
> 2) re-organize code to reduce copying of the array.
> 3) add optimization for a left/right shift that
>    returns the original array.
> 4) Don't leak memory.
>
> --
> Steve



-- 
Outside of a dog, a book is a man's best friend. Inside of a dog it's
too dark to read.

Groucho Marx

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

On Sat, Nov 21, 2015 at 11:41:51AM +0100, Paul Richard Thomas wrote:
> 
> Just a couple of small typos:
> "Unexpected expr_type cause an ICE" ;  causes?
> "! An array of derived types workd too." ; works?
> 
> Apart from that it's OK for trunk.
> 
> Thanks for the patch
> 

Thanks for the the review.  I don't have a clue as
to how to do simplification for rank > 2. :(

-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[H.J. Lu]

On Sat, Nov 21, 2015 at 8:26 AM, Steve Kargl
<sgk@troutmask.apl.washington.edu> wrote:
> On Sat, Nov 21, 2015 at 11:41:51AM +0100, Paul Richard Thomas wrote:
>>
>> Just a couple of small typos:
>> "Unexpected expr_type cause an ICE" ;  causes?
>> "! An array of derived types workd too." ; works?
>>
>> Apart from that it's OK for trunk.
>>
>> Thanks for the patch
>>
>
> Thanks for the the review.  I don't have a clue as
> to how to do simplification for rank > 2. :(
>

It breaks bootstrap:

  int dm;

  /* DIM is only useful for rank > 1, but deal with it here as one can
     set DIM = 1 for rank = 1.  */
  if (dim)
    {
      if (!gfc_is_constant_expr (dim))
return NULL;
      dm = mpz_get_si (dim->value.integer);
    }
  else
    dm = 1;

dm is set, but never used.

H.J.

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

On Sat, Nov 21, 2015 at 10:07:35AM -0800, H.J. Lu wrote:
> On Sat, Nov 21, 2015 at 8:26 AM, Steve Kargl
> <sgk@troutmask.apl.washington.edu> wrote:
> > On Sat, Nov 21, 2015 at 11:41:51AM +0100, Paul Richard Thomas wrote:
> >>
> >> Just a couple of small typos:
> >> "Unexpected expr_type cause an ICE" ;  causes?
> >> "! An array of derived types workd too." ; works?
> >>
> >> Apart from that it's OK for trunk.
> >>
> >> Thanks for the patch
> >>
> >
> > Thanks for the the review.  I don't have a clue as
> > to how to do simplification for rank > 2. :(
> >
> 
> It breaks bootstrap:
> 
>   int dm;
> 
>   /* DIM is only useful for rank > 1, but deal with it here as one can
>      set DIM = 1 for rank = 1.  */
>   if (dim)
>     {
>       if (!gfc_is_constant_expr (dim))
> return NULL;
>       dm = mpz_get_si (dim->value.integer);
>     }
>   else
>     dm = 1;
> 
> dm is set, but never used.
> 

Perhaps, bootstrap needs to set appropriate warning levels.

-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[H.J. Lu]

On Sat, Nov 21, 2015 at 10:20 AM, Steve Kargl
<sgk@troutmask.apl.washington.edu> wrote:
> On Sat, Nov 21, 2015 at 10:07:35AM -0800, H.J. Lu wrote:
>> On Sat, Nov 21, 2015 at 8:26 AM, Steve Kargl
>> <sgk@troutmask.apl.washington.edu> wrote:
>> > On Sat, Nov 21, 2015 at 11:41:51AM +0100, Paul Richard Thomas wrote:
>> >>
>> >> Just a couple of small typos:
>> >> "Unexpected expr_type cause an ICE" ;  causes?
>> >> "! An array of derived types workd too." ; works?
>> >>
>> >> Apart from that it's OK for trunk.
>> >>
>> >> Thanks for the patch
>> >>
>> >
>> > Thanks for the the review.  I don't have a clue as
>> > to how to do simplification for rank > 2. :(
>> >
>>
>> It breaks bootstrap:
>>
>>   int dm;
>>
>>   /* DIM is only useful for rank > 1, but deal with it here as one can
>>      set DIM = 1 for rank = 1.  */
>>   if (dim)
>>     {
>>       if (!gfc_is_constant_expr (dim))
>> return NULL;
>>       dm = mpz_get_si (dim->value.integer);
>>     }
>>   else
>>     dm = 1;
>>
>> dm is set, but never used.
>>
>
> Perhaps, bootstrap needs to set appropriate warning levels.

https://gcc.gnu.org/ml/gcc-regression/2015-11/msg00648.html

-- 
H.J.

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

On Sat, Nov 21, 2015 at 11:19:22AM -0800, H.J. Lu wrote:
> On Sat, Nov 21, 2015 at 10:20 AM, Steve Kargl
> <sgk@troutmask.apl.washington.edu> wrote:
> > On Sat, Nov 21, 2015 at 10:07:35AM -0800, H.J. Lu wrote:
> >> On Sat, Nov 21, 2015 at 8:26 AM, Steve Kargl
> >> <sgk@troutmask.apl.washington.edu> wrote:
> >> > On Sat, Nov 21, 2015 at 11:41:51AM +0100, Paul Richard Thomas wrote:
> >> >>
> >> >> Just a couple of small typos:
> >> >> "Unexpected expr_type cause an ICE" ;  causes?
> >> >> "! An array of derived types workd too." ; works?
> >> >>
> >> >> Apart from that it's OK for trunk.
> >> >>
> >> >> Thanks for the patch
> >> >>
> >> >
> >> > Thanks for the the review.  I don't have a clue as
> >> > to how to do simplification for rank > 2. :(
> >> >
> >>
> >> It breaks bootstrap:
> >>
> >>   int dm;
> >>
> >>   /* DIM is only useful for rank > 1, but deal with it here as one can
> >>      set DIM = 1 for rank = 1.  */
> >>   if (dim)
> >>     {
> >>       if (!gfc_is_constant_expr (dim))
> >> return NULL;
> >>       dm = mpz_get_si (dim->value.integer);
> >>     }
> >>   else
> >>     dm = 1;
> >>
> >> dm is set, but never used.
> >>
> >
> > Perhaps, bootstrap needs to set appropriate warning levels.
> 
> https://gcc.gnu.org/ml/gcc-regression/2015-11/msg00648.html
> 

See 5 lines up.

-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

> > >
> > > Perhaps, bootstrap needs to set appropriate warning levels.
> > 
> > https://gcc.gnu.org/ml/gcc-regression/2015-11/msg00648.html
> > 
> 
> See 5 lines up.
> 

Committed.

Index: ChangeLog
===================================================================
--- ChangeLog	(revision 230709)
+++ ChangeLog	(working copy)
@@ -1,5 +1,10 @@
 2015-11-21  Steven G. Kargl  <kargl@gcc.gnu.org>
 
+	* simplify.c (gfc_simplify_cshift): Work around bootstrap issues
+	due to inappropriate warning options. 
+
+2015-11-21  Steven G. Kargl  <kargl@gcc.gnu.org>
+
 	* simplify.c (gfc_simplify_cshift): Implement simplification of
 	CSHIFT for rank=1 arrays.
 	(gfc_simplify_spread): Remove a FIXME and add error condition.
Index: simplify.c
===================================================================
--- simplify.c	(revision 230709)
+++ simplify.c	(working copy)
@@ -1869,6 +1869,15 @@ gfc_simplify_cshift (gfc_expr *array, gf
   else
     {
       /* FIXME: Deal with rank > 1 arrays.  For now, don't leak memory.  */
+
+      /* GCC bootstrap is too stupid to realize that the above code for dm
+	 is correct.  First, dim can be specified for a rank 1 array.  It is
+	 not needed in this nor used here.  Second, the code is simply waiting
+	 for someone to implement rank > 1 simplification.   For now, add a
+	 pessimization to the code that has a zero valid reason to be here.  */
+      if (dm > array->rank)
+	gcc_unreachable ();
+
       gfc_free_expr (a);
     }
 
-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Eric Botcazou]

> +	* simplify.c (gfc_simplify_cshift): Work around bootstrap issues
> +	due to inappropriate warning options.

The warning options are appropriate, any dead code can potentially hide a bug 
and should be flagged; every static analyzer will also do it and we would soon 
have PRs opened with bugzilla if we tolerated it.

-- 
Eric Botcazou

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

On Sat, Nov 21, 2015 at 11:26:17PM +0100, Eric Botcazou wrote:
> > +	* simplify.c (gfc_simplify_cshift): Work around bootstrap issues
> > +	due to inappropriate warning options.
> 
> The warning options are appropriate, any dead code can potentially hide a bug 
> and should be flagged; every static analyzer will also do it and we would soon 
> have PRs opened with bugzilla if we tolerated it.
> 

I disgree.

If bootstrap is going to enforce -Werror -Wunused-*, then 
this should be the default for any possible invocation of
make.

-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Gerald Pfeifer]

On Sat, 21 Nov 2015, Steve Kargl wrote:
>  2015-11-21  Steven G. Kargl  <kargl@gcc.gnu.org>
>  
> 	* simplify.c (gfc_simplify_cshift): Work around bootstrap issues
> 	due to inappropriate warning options. 

> Index: simplify.c
> ===================================================================
> +      /* GCC bootstrap is too stupid to realize that the above code for dm
> +	 is correct.  First, dim can be specified for a rank 1 array.  It is
> +	 not needed in this nor used here.  Second, the code is simply waiting
> +	 for someone to implement rank > 1 simplification.   For now, add a
> +	 pessimization to the code that has a zero valid reason to be here.  */
> +      if (dm > array->rank)
> +	gcc_unreachable ();

I'm not sure this comment is appropriate as is.  At a minimum, it
should list the version of GCC this was introduced for/with.  So,
something like

  "As of GCC 6, when bootstrapping we do not realize that..."

Gerald

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

On Thu, Dec 31, 2015 at 09:57:10PM +0800, Gerald Pfeifer wrote:
> On Sat, 21 Nov 2015, Steve Kargl wrote:
> >  2015-11-21  Steven G. Kargl  <kargl@gcc.gnu.org>
> >  
> > 	* simplify.c (gfc_simplify_cshift): Work around bootstrap issues
> > 	due to inappropriate warning options. 
> 
> > Index: simplify.c
> > ===================================================================
> > +      /* GCC bootstrap is too stupid to realize that the above code for dm
> > +	 is correct.  First, dim can be specified for a rank 1 array.  It is
> > +	 not needed in this nor used here.  Second, the code is simply waiting
> > +	 for someone to implement rank > 1 simplification.   For now, add a
> > +	 pessimization to the code that has a zero valid reason to be here.  */
> > +      if (dm > array->rank)
> > +	gcc_unreachable ();
> 
> I'm not sure this comment is appropriate as is.  At a minimum, it
> should list the version of GCC this was introduced for/with.  So,
> something like
> 

I have no intention to change the comment.

-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Dominique d'Humières]

‘dm’ is actually not used, the building problem is fixed by the patch (I did not rearrange the nested ‘if’s)

--- ../_clean/gcc/fortran/simplify.c	2015-11-21 20:59:57.000000000 +0100
+++ gcc/fortran/simplify.c	2015-11-21 21:06:30.000000000 +0100
@@ -1792,7 +1792,6 @@ gfc_expr *
 gfc_simplify_cshift (gfc_expr *array, gfc_expr *shift, gfc_expr *dim)
 {
   gfc_expr *a, *result;
-  int dm;
 
   /* DIM is only useful for rank > 1, but deal with it here as one can
      set DIM = 1 for rank = 1.  */
@@ -1800,10 +1799,7 @@ gfc_simplify_cshift (gfc_expr *array, gf
     {
       if (!gfc_is_constant_expr (dim))
 	return NULL;
-      dm = mpz_get_si (dim->value.integer);
     }
-  else
-    dm = 1;
 
   /* Copy array into 'a', simplify it, and then test for a constant array.
      An unexpected expr_type causes an ICE.   */

Dominique

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

On Sat, Nov 21, 2015 at 09:22:40PM +0100, Dominique d'Humi??res wrote:
> ???dm??? is actually not used, the building problem is fixed by the patch (I did not rearrange the nested ???if???s)
> 
> --- ../_clean/gcc/fortran/simplify.c	2015-11-21 20:59:57.000000000 +0100
> +++ gcc/fortran/simplify.c	2015-11-21 21:06:30.000000000 +0100
> @@ -1792,7 +1792,6 @@ gfc_expr *
>  gfc_simplify_cshift (gfc_expr *array, gfc_expr *shift, gfc_expr *dim)
>  {
>    gfc_expr *a, *result;
> -  int dm;
>  
>    /* DIM is only useful for rank > 1, but deal with it here as one can
>       set DIM = 1 for rank = 1.  */
> @@ -1800,10 +1799,7 @@ gfc_simplify_cshift (gfc_expr *array, gf
>      {
>        if (!gfc_is_constant_expr (dim))
>  	return NULL;
> -      dm = mpz_get_si (dim->value.integer);
>      }
> -  else
> -    dm = 1;
>  
>    /* Copy array into 'a', simplify it, and then test for a constant array.
>       An unexpected expr_type causes an ICE.   */
> 

I already fix the problem.

Your patch unfixes the valid code that is needed for the rank>2 case.

-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Dominique d'Humières]

[-- Attachment #1: Type: text/plain, Size: 213 bytes --]

Hi Steve,

Compiling the attached code after revision r230710 gives the ICE

f951: internal compiler error: in gfc_simplify_cshift, at fortran/simplify.c:1823

while it compiled before.

TIA

Dominique

[-- Attachment #2: mhd.f90 --]
[-- Type: application/octet-stream, Size: 11017 bytes --]

! TVD split MHD code
! copyright (C) 2001,2003, Ue-Li Pen
! written November 2001 by Ue-Li Pen, pen@cita.utoronto.ca
!This program is free software; you can redistribute it and/or
!modify it under the terms of the GNU General Public License
!as published by the Free Software Foundation; either version 2
!of the License, or (at your option) any later version.
!
!This program is distributed in the hope that it will be useful,
!but WITHOUT ANY WARRANTY; without even the implied warranty of
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!GNU General Public License for more details.
!
!You should have received a copy of the GNU General Public License
!along with this program; if not, write to the Free Software
!Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!
! debugged to pass alven test March 29, 2002
!
! release 0.1 May 2003
! release 0.2 Sept 2009 : fixed b2 interpolation in cfl
!
! General Notes:
! restrictions: requires nx=nz or nz=1
! see http://arxiv.org/astroph/abs/astro-ph/0305088
! or http://www.cita.utoronto.ca/~pen/MHD
! for questions contact pen@cita.utoronto.ca
!
implicit none
integer, parameter :: n=100,nx=n,ny=1,nz=1
real, dimension(5,nx,ny,nz) :: u
real, dimension(3,nx,ny,nz) :: b
real t,dt,tf
integer iter,i
! B field is stored on the left side of each cell


call init(u,b,nx,ny,nz)
tf=nx*4
t=0
iter=0
do
   iter=iter+1
   if (t>=tf) exit
   dt=0.9*cfl(u,b)
   dt=min(dt,(tf-t)/2)
   t=t+2*dt
   call fluidx(u,b,nx,ny,nz,dt)
   call advectbyzx(u,b,nx,ny,nz,dt)
! the y sweep
   call transpose12(u,b,u,b,nx,ny,nz)
   call fluidx(u,b,ny,nx,nz,dt)
   call advectbyzx(u,b,ny,nx,nz,dt)
! z sweep
   call transpose13(u,b,u,b,ny,nx,nz)
   call fluidx(u,b,nz,nx,ny,dt)
   call advectbyzx(u,b,nz,nx,ny,dt)
   call advectbyzx(u,b,nz,nx,ny,dt)
   call fluidx(u,b,nz,nx,ny,dt)
   
! back
   call transpose13(u,b,u,b,nz,nx,ny)
   call advectbyzx(u,b,ny,nx,nz,dt)
   call fluidx(u,b,ny,nx,nz,dt)
! x again
   call transpose12(u,b,u,b,ny,nx,nz)
   call advectbyzx(u,b,nx,ny,nz,dt)
   call fluidx(u,b,nx,ny,nz,dt)
   if (mod(iter,10) .eq. 1) write(*,*) 't=',t,iter,u(2,nx/4,1,1)
end do
call output

contains
  subroutine output
    integer i,j
    open(10,file='u.dat')
    open(20,file='e.dat')
    do j=1,1
    do i=1,nx
       write(10,'(i5,7(1x,e9.3))') i,u(1,i,j,1)-1,u(2:4,i,j,1)/u(1,i,j,1),b(:,i,j,1)-(/1,0,0/)
    end do
    do i=1,nx
       write(20,'(i5,7(1x,e30.20))') i,u(5,i,j,1)-sum(u(2:4,i,j,1)**2)/u(1,i,j,1)/2-sum(b(:,i,j,1)**2)/2
    end do
    end do
    close(20)
    close(10)
  end subroutine output


  function cfl(u,b)
    real, dimension(:,:,:,:)::u,b
    real cfl
! locals
    real v,ps,p,c,gamma,bx,by,bz,b2
    integer i,j,k,ip,jp,kp

    gamma=5./3.
    c=0
    !$omp parallel do private(j,i,kp,jp,ip,bx,by,bz,v,ps,p,b2) reduction(max:c)
    do k=1,nz
       do j=1,ny
          do i=1,nx
             kp=mod(k,nz)+1
             jp=mod(j,ny)+1
             ip=mod(i,nx)+1
             bx=(b(1,i,j,k)+b(1,ip,j,k))/2
             by=(b(2,i,j,k)+b(2,i,jp,k))/2
             bz=(b(3,i,j,k)+b(3,i,j,kp))/2
             v=maxval(abs(u(2:4,i,j,k)/u(1,i,j,k)))
             b2=bx*bx+by*by+bz*bz 
             ps=(u(5,i,j,k)-sum(u(2:4,i,j,k)**2,1)/u(1,i,j,k)/2)*(gamma-1)+(2-gamma)*b2/2
             p=ps-b2/2
!   c=max(c,v+sqrt(abs(  (sum(b(:,i,j,k)**2)*(gamma-1)+gamma*p)/u(1,i,j,k))))
             c=max(c,v+sqrt(abs(  (b2*2+gamma*p)/u(1,i,j,k))))
          end do
       end do
    end do
    cfl=1/c
  end function cfl


end

subroutine init(u,b,nx,ny,nz)
  implicit none
  integer nx,ny,nz
  real u(5,nx,ny,nz),b(3,nx,ny,nz)
  real p0
  integer i
  real, parameter :: epsilon=.1
  
  p0=3./5.
  u=0
  b=0
  b(1,:,:,:)=1
  u(1,:,:,:)=1
  u(5,:,:,:)=1.5*p0
  u(2,:,1,1)=0.0001*sin( (/ (2*3.14159*i/nx,i=1,nx) /) )
  u(1,:,1,1)=u(1,:,1,1)+u(2,:,1,1)
  u(5,:,:,:)=u(5,:,:,:)+sum(b**2,1)/2+u(2,:,:,:)**2/u(1,:,:,:)/2
  u(5,:,1,1)=u(5,:,1,1)+p0*u(2,:,1,1)*5./3./(2./3.)
!    return
! circularly polarized alven wave:
! background : rho=B_x=1 all others zero
  u=0
  b=0
  b(1,:,:,:)=1
  u(1,:,:,:)=1
  u(5,:,:,:)=0.1  ! to keep things stable
  do i=1,nx
     u(3,i,:,:)=epsilon*sin( 2*3.14159*i/nx )
     u(4,i,:,:)=epsilon*cos( 2*3.14159*i/nx )
  enddo
  b(2,:,:,:)=-u(3,:,:,:)
  b(3,:,:,:)=-u(4,:,:,:)
  b(2,:,:,:)=(b(2,:,:,:)+cshift(b(2,:,:,:),-1,2))/2
  b(3,:,:,:)=(b(3,:,:,:)+cshift(b(3,:,:,:),-1,3))/2
  
  u(5,:,:,:)=u(5,:,:,:)+sum(b**2,1)/2+sum(u(2:4,:,:,:)**2,1)/u(1,:,:,:)/2
!  return
  
! alven wave:
! background : rho=B_x=1 all others zero
! \dot{v_y}+\div_x(-B_y)=0
! \dot{B_y}+\div_x (       -   v_y)     = 0
! let v_y=\epsilon sin(2 \pi (x-t)/L)
! then B_y=-v_y
  u=0
  b=0
  b(1,:,:,:)=1
  u(1,:,:,:)=1
  u(5,:,:,:)=.001  ! to keep things stable
  do i=1,nx
     u(3,i,:,:)=0.1*sin( 2*3.14159*i/nx )
  enddo
  b(2,:,:,:)=-u(3,:,:,:)
  b(2,:,:,:)=(b(2,:,:,:)+cshift(b(2,:,:,:),-1,1))/2
  u(5,:,:,:)=u(5,:,:,:)+sum(b**2,1)/2+u(3,:,:,:)**2/u(1,:,:,:)/2
 ! return
  ! magnetic
  u=0
  b=0
  b(2,:,:,:)=1
  u(1,:,:,:)=1
  u(5,:,:,:)=.000000001  ! to keep things stable
  do i=1,nx
     u(2,i,:,:)=0.001*sin( 2*3.14159*i/nx )
  enddo
  b(2,:,:,:)=b(2,:,:,:)+u(2,:,:,:)
  u(1,:,:,:)=b(2,:,:,:)
  u(5,:,:,:)=u(5,:,:,:)+sum(b**2,1)/2+sum(u(2:4,:,:,:)**2,1)/u(1,:,:,:)/2
  write(*,*) 'magnetic'
!  return
  ! magnetosonic
  p0=3./5./2.
  u=0
  b=0
  b(2,:,:,:)=1./sqrt(2.)
  u(1,:,:,:)=1
  u(5,:,:,:)=p0*1.5
  do i=1,nx
     u(2,i,:,:)=0.0001*sin( 2*3.14159*i/nx )
  enddo
  b(2,:,:,:)=b(2,:,:,:)+u(2,:,:,:)/sqrt(2.)
  u(1,:,:,:)=u(1,:,:,:)+u(2,:,:,:)
  u(5,:,:,:)=u(5,:,:,:)+sum(b**2,1)/2+sum(u(2:4,:,:,:)**2,1)/u(1,:,:,:)/2
  u(5,:,:,:)=u(5,:,:,:)+p0*u(2,:,:,:)*5./3./(2./3.)
  write(*,*) 'magnetosonic'
end subroutine init

subroutine fluidx(u,b,nx,ny,nz,dt)
  implicit none
  integer nx,ny,nz
  real u(5,nx,ny,nz),b(3,nx,ny,nz),dt
  real, dimension(3,nx) :: b3x
  real, dimension(5,nx) :: u1x
  integer j,k,jp,kp

  !$omp parallel do private(j,u1x,b3x,jp,kp)
  do k=1,nz
     do j=1,ny
        b3x=b(:,:,j,k)/2
        jp=mod(j,ny)+1
        kp=mod(k,nz)+1
        b3x(1,:)=b3x(1,:)+cshift(b3x(1,:),1)
        b3x(2,:)=b3x(2,:)+b(2,:,jp,k)/2
        b3x(3,:)=b3x(3,:)+b(3,:,j,kp)/2
        call tvd1(u(1,1,j,k),b3x,nx,dt)
     end do
  end do
end subroutine fluidx

subroutine advectbyzx(u,b,nx,ny,nz,dt)
  implicit none
  integer nx,ny,nz
  real u(5,nx,ny,nz),b(3,nx,ny,nz),dt
  real, dimension(nx) :: fluxbx,b1x,vx
  integer j,k,jm,km

  !$omp parallel do private(j,jm,vx,b1x,fluxbx)
  do k=1,nz
     do j=1,ny
        jm=mod(j+ny-2,ny)+1
        vx=(u(2,:,jm,k)+u(2,:,j,k))/(u(1,:,jm,k)+u(1,:,j,k))
        vx=(cshift(vx,-1)+cshift(vx,1)+2*vx)/4
        b1x=b(2,:,j,k)
        call tvdb(fluxbx,b1x,vx,nx,dt)
        b(2,:,j,k)=b1x
        b(1,:,j,k)=b(1,:,j,k)-cshift(fluxbx,-1)
        b(1,:,jm,k)=b(1,:,jm,k)+cshift(fluxbx,-1)
     end do
  end do
  do k=1,nz
     !$omp parallel do private(km,vx,b1x,fluxbx)
     do j=1,ny
        km=mod(k+nz-2,nz)+1
        vx=(u(2,:,j,km)+u(2,:,j,k))/(u(1,:,j,km)+u(1,:,j,k))
        vx=(cshift(vx,-1)+cshift(vx,1)+2*vx)/4
        b1x=b(3,:,j,k)
        call tvdb(fluxbx,b1x,vx,nx,dt)
        b(3,:,j,k)=b1x
        b(1,:,j,k)=b(1,:,j,k)-cshift(fluxbx,-1)
        b(1,:,j,km)=b(1,:,j,km)+cshift(fluxbx,-1)
     end do
  end do
end subroutine advectbyzx

recursive subroutine tvdb(flux,b,vg,n,dt)
  integer i,n,ip,ipp,im
  real dt
  real, dimension(n) :: flux,b,vg
! locals
  real, dimension(n) :: b1,flux1,vh
  real w,wp,wm,dw,v
  
  ! unlike the B field, the flux lives on the right cell boundary
  vh=(vg+cshift(vg,1))/2
  where(vh>0)
     flux1=b*vg
  elsewhere
     flux1=cshift(b*vg,1)
  end where
  b1=b-(flux1-cshift(flux1,-1))*dt/2
  do i=1,n
     ip=mod(i,n)+1
     ipp=mod(ip,n)+1
     im=mod(i+n-2,n)+1
     v=vh(i)
     if (v>0) then
        w=vg(i)*b1(i)
        wp=(vg(ip)*b1(ip)-w)/2
        wm=(w-vg(im)*b1(im))/2
     else
        w=vg(ip)*b1(ip)
        wp=(w-vg(ipp)*b1(ipp))/2
        wm=(vg(i)*b1(i)-w)/2
     end if
     dw=0
     if(wm*wp>0) dw=2*wm*wp/(wm+wp)
     flux(i)=(w+dw)*dt
  end do
  b=b-(flux-cshift(flux,-1))
end subroutine tvdb

recursive subroutine tvd1(u,b,n,dt)
  implicit none
  integer n
  real dt
  real u(5,n),b(3,n)
!locals    
  real, dimension(5,n) :: v, u1,wr,wl,fr,fl,flux,dfrp,dfrm&
       ,dflm,dflp,dfl,dfr
  real c

  call mhdflux(v,c,u,b,n)
  wr=u+v
  wl=u-v
  fr=c*wr
  fl=cshift(c*wl,1,2)
  flux=(fr-fl)/2
  u1=u-(flux-cshift(flux,-1,2))*dt/2
  call mhdflux(v,c,u1,b,n)
  wr=u1+v
  wl=u1-v
  fr=c*wr
  dfrp=(cshift(fr,1,2)-fr)/2
  dfrm=(fr-cshift(fr,-1,2))/2
  dfr=0
  where(dfrp*dfrm>0)
     dfr=2*dfrp*dfrm/(dfrp+dfrm)
  end where
  fl=cshift(c*wl,1,2)
  dflp=(fl-cshift(fl,1,2))/2
  dflm=(cshift(fl,-1,2)-fl)/2
  dfl=0
  where(dflp*dflm>0)
     dfl=2*dflp*dflm/(dflp+dflm)
  end where
  flux=(fr-fl+(dfr-dfl))/2
  u=u-(flux-cshift(flux,-1,2))*dt
  return
end subroutine tvd1



recursive subroutine mhdflux(v,c,u,b,n)
  implicit none
  integer n
  real, dimension(5,n)::v,u
  real b(3,n)
  real, intent(out) :: c
  ! locals
  real, dimension(n) :: vx,ps,p
  real gamma

  gamma=5./3.

  vx=u(2,:)/u(1,:)
  ps=(u(5,:)-sum(u(2:4,:)**2,1)/u(1,:)/2)*(gamma-1)+(2-gamma)*sum(b**2,1)/2
  v(1,:)=u(2,:)
  v(2,:)=u(2,:)*vx+ps-b(1,:)**2
  v(3,:)=u(3,:)*vx-b(2,:)*b(1,:)
  v(4,:)=u(4,:)*vx-b(3,:)*b(1,:)
  v(5,:)=(u(5,:)+ps)*vx-b(1,:)*sum(b*u(2:4,:),1)/u(1,:)
  p=ps-sum(b**2,1)/2
  c=maxval(abs(vx)+sqrt(abs( (sum(b**2,1)+gamma*p)/u(1,:))))
  if (c>0) v=v/c
end subroutine mhdflux



subroutine transpose12(ut,bt,u,b,nx,ny,nz)
  implicit none
  integer nx,ny,nz
  real u(5,nx,ny,nz),b(3,nx,ny,nz)
  real ut(5,ny,nx,nz),bt(3,ny,nx,nz)

  integer i,j,k
  real u2(5,ny,nx),b2(3,ny,nx)

!$omp parallel do default(none) shared(u,b,ut,bt,ny,nx,nz) private(i,j,u2,b2)
  do k=1,nz
     do j=1,ny
        do i=1,nx
           u2(:,j,i)=u((/1,3,2,4,5/),i,j,k)
           b2(:,j,i)=b((/2,1,3/),i,j,k)
        end do
     end do
     ut(:,:,:,k)=u2
     bt(:,:,:,k)=b2
  end do
end subroutine transpose12


subroutine transpose13(ut,bt,u,b,nx,ny,nz)
  implicit none
  integer nx,ny,nz
  real u(5,nx,ny,nz),b(3,nx,ny,nz)
  real ut(5,nz,ny,nx),bt(3,nz,ny,nx)

  integer i,j,k
  real u2(5,nz,nx),b2(3,nz,nx)

  if (nx .eq. nz) then
!$omp parallel do default(none) shared(u,b,ut,bt,ny,nx,nz) private(i,k,u2,b2)
  do j=1,ny
     do k=1,nz
        do i=1,nx
           u2(:,k,i)=u((/1,4,3,2,5/),i,j,k)
           b2(:,k,i)=b((/3,2,1/),i,j,k)
        end do
     end do
     ut(:,:,j,:)=u2
     bt(:,:,j,:)=b2
  end do
  else if (nz .eq. 1) then
     call transpose12(ut,bt,u,b,nx,ny,nz)
     do i=1,nx
        do j=1,ny
           ut(:,1,j,i)=ut((/1,4,2,3,5/),1,j,i)
           bt(:,1,j,i)=bt((/3,1,2/),1,j,i)
        end do
     end do
  else if (nx .eq. 1) then
     call transpose12(ut,bt,u,b,ny,nz,nx)
     do i=1,ny
        do j=1,nz
           ut(:,j,i,1)=ut((/1,4,2,3,5/),j,i,1)
           bt(:,j,i,1)=bt((/3,1,2/),j,i,1)
        end do
     end do
  else
     write(*,*) 'nz<>nx not supported'
     stop
  endif
end subroutine transpose13

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Steve Kargl]

On Sun, Nov 22, 2015 at 09:51:09AM +0100, Dominique d'Humi??res wrote:
> 
> Compiling the attached code after revision r230710 gives the ICE
> 
> f951: internal compiler error: in gfc_simplify_cshift, at fortran/simplify.c:1823
> 
> while it compiled before.
> 

File a bug report, then fix rank > 2.

-- 
Steve

Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT

[Dominique d'Humières]

> Le 22 nov. 2015 à 17:21, Steve Kargl <sgk@troutmask.apl.washington.edu> a écrit :
> 
> On Sun, Nov 22, 2015 at 09:51:09AM +0100, Dominique d'Humi??res wrote:
>> 
>> Compiling the attached code after revision r230710 gives the ICE
>> 
>> f951: internal compiler error: in gfc_simplify_cshift, at fortran/simplify.c:1823
>> 
>> while it compiled before.
>> 
> File a bug report, then fix rank > 2.

Already done by H.J. Lu (pr68486), see also my comment 2.

Dominique

> 
> -- 
> Steve