public inbox for gcc-patches@gcc.gnu.org
 help / color / mirror / Atom feed
* [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
@ 2015-11-20  0:58 Steve Kargl
  2015-11-20  1:31 ` Steve Kargl
  2015-11-20 20:09 ` Steve Kargl
  0 siblings, 2 replies; 20+ messages in thread
From: Steve Kargl @ 2015-11-20  0:58 UTC (permalink / raw)
  To: fortran, gcc-patches

[-- 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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-20  0:58 [PATCH] (Partial) Implementation of simplificaiton of CSHIFT Steve Kargl
@ 2015-11-20  1:31 ` Steve Kargl
  2015-11-20  3:16   ` Steve Kargl
  2015-11-20 20:09 ` Steve Kargl
  1 sibling, 1 reply; 20+ messages in thread
From: Steve Kargl @ 2015-11-20  1:31 UTC (permalink / raw)
  To: fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-20  1:31 ` Steve Kargl
@ 2015-11-20  3:16   ` Steve Kargl
  0 siblings, 0 replies; 20+ messages in thread
From: Steve Kargl @ 2015-11-20  3:16 UTC (permalink / raw)
  To: fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-20  0:58 [PATCH] (Partial) Implementation of simplificaiton of CSHIFT Steve Kargl
  2015-11-20  1:31 ` Steve Kargl
@ 2015-11-20 20:09 ` Steve Kargl
  2015-11-21 11:45   ` Paul Richard Thomas
  1 sibling, 1 reply; 20+ messages in thread
From: Steve Kargl @ 2015-11-20 20:09 UTC (permalink / raw)
  To: fortran, gcc-patches

[-- 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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-20 20:09 ` Steve Kargl
@ 2015-11-21 11:45   ` Paul Richard Thomas
  2015-11-21 18:07     ` Steve Kargl
  0 siblings, 1 reply; 20+ messages in thread
From: Paul Richard Thomas @ 2015-11-21 11:45 UTC (permalink / raw)
  To: Steve Kargl; +Cc: fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 11:45   ` Paul Richard Thomas
@ 2015-11-21 18:07     ` Steve Kargl
  2015-11-21 18:12       ` H.J. Lu
  0 siblings, 1 reply; 20+ messages in thread
From: Steve Kargl @ 2015-11-21 18:07 UTC (permalink / raw)
  To: Paul Richard Thomas; +Cc: fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 18:07     ` Steve Kargl
@ 2015-11-21 18:12       ` H.J. Lu
  2015-11-21 18:21         ` Steve Kargl
  0 siblings, 1 reply; 20+ messages in thread
From: H.J. Lu @ 2015-11-21 18:12 UTC (permalink / raw)
  To: Steve Kargl; +Cc: Paul Richard Thomas, fortran, gcc-patches

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.

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 18:12       ` H.J. Lu
@ 2015-11-21 18:21         ` Steve Kargl
  2015-11-21 19:26           ` H.J. Lu
  0 siblings, 1 reply; 20+ messages in thread
From: Steve Kargl @ 2015-11-21 18:21 UTC (permalink / raw)
  To: H.J. Lu; +Cc: Paul Richard Thomas, fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 18:21         ` Steve Kargl
@ 2015-11-21 19:26           ` H.J. Lu
  2015-11-21 20:00             ` Steve Kargl
  0 siblings, 1 reply; 20+ messages in thread
From: H.J. Lu @ 2015-11-21 19:26 UTC (permalink / raw)
  To: Steve Kargl; +Cc: Paul Richard Thomas, fortran, gcc-patches

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.

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 19:26           ` H.J. Lu
@ 2015-11-21 20:00             ` Steve Kargl
  2015-11-21 20:22               ` Steve Kargl
  0 siblings, 1 reply; 20+ messages in thread
From: Steve Kargl @ 2015-11-21 20:00 UTC (permalink / raw)
  To: H.J. Lu; +Cc: Paul Richard Thomas, fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 20:00             ` Steve Kargl
@ 2015-11-21 20:22               ` Steve Kargl
  2015-11-21 22:39                 ` Eric Botcazou
  2015-12-31 14:00                 ` Gerald Pfeifer
  0 siblings, 2 replies; 20+ messages in thread
From: Steve Kargl @ 2015-11-21 20:22 UTC (permalink / raw)
  To: H.J. Lu; +Cc: Paul Richard Thomas, fortran, gcc-patches

> > >
> > > 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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 20:22               ` Steve Kargl
@ 2015-11-21 22:39                 ` Eric Botcazou
  2015-11-22  1:53                   ` Steve Kargl
  2015-12-31 14:00                 ` Gerald Pfeifer
  1 sibling, 1 reply; 20+ messages in thread
From: Eric Botcazou @ 2015-11-21 22:39 UTC (permalink / raw)
  To: Steve Kargl; +Cc: gcc-patches, H.J. Lu, Paul Richard Thomas, fortran

> +	* 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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 22:39                 ` Eric Botcazou
@ 2015-11-22  1:53                   ` Steve Kargl
  0 siblings, 0 replies; 20+ messages in thread
From: Steve Kargl @ 2015-11-22  1:53 UTC (permalink / raw)
  To: Eric Botcazou; +Cc: gcc-patches, H.J. Lu, Paul Richard Thomas, fortran

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 20:22               ` Steve Kargl
  2015-11-21 22:39                 ` Eric Botcazou
@ 2015-12-31 14:00                 ` Gerald Pfeifer
  2015-12-31 16:13                   ` Steve Kargl
  1 sibling, 1 reply; 20+ messages in thread
From: Gerald Pfeifer @ 2015-12-31 14:00 UTC (permalink / raw)
  To: Steve Kargl; +Cc: H.J. Lu, Paul Richard Thomas, fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-12-31 14:00                 ` Gerald Pfeifer
@ 2015-12-31 16:13                   ` Steve Kargl
  0 siblings, 0 replies; 20+ messages in thread
From: Steve Kargl @ 2015-12-31 16:13 UTC (permalink / raw)
  To: Gerald Pfeifer; +Cc: H.J. Lu, Paul Richard Thomas, fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-22 16:57     ` Steve Kargl
@ 2015-11-22 18:41       ` Dominique d'Humières
  0 siblings, 0 replies; 20+ messages in thread
From: Dominique d'Humières @ 2015-11-22 18:41 UTC (permalink / raw)
  To: Steve Kargl; +Cc: fortran, gcc-patches


> 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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-22  9:14   ` Dominique d'Humières
@ 2015-11-22 16:57     ` Steve Kargl
  2015-11-22 18:41       ` Dominique d'Humières
  0 siblings, 1 reply; 20+ messages in thread
From: Steve Kargl @ 2015-11-22 16:57 UTC (permalink / raw)
  To: Dominique d'Humi??res; +Cc: fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 21:50 ` Steve Kargl
@ 2015-11-22  9:14   ` Dominique d'Humières
  2015-11-22 16:57     ` Steve Kargl
  0 siblings, 1 reply; 20+ messages in thread
From: Dominique d'Humières @ 2015-11-22  9:14 UTC (permalink / raw)
  To: Steve Kargl; +Cc: fortran, gcc-patches

[-- 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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
  2015-11-21 20:30 Dominique d'Humières
@ 2015-11-21 21:50 ` Steve Kargl
  2015-11-22  9:14   ` Dominique d'Humières
  0 siblings, 1 reply; 20+ messages in thread
From: Steve Kargl @ 2015-11-21 21:50 UTC (permalink / raw)
  To: Dominique d'Humi??res; +Cc: fortran, gcc-patches

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

^ permalink raw reply	[flat|nested] 20+ messages in thread

* Re: [PATCH] (Partial) Implementation of simplificaiton of CSHIFT
@ 2015-11-21 20:30 Dominique d'Humières
  2015-11-21 21:50 ` Steve Kargl
  0 siblings, 1 reply; 20+ messages in thread
From: Dominique d'Humières @ 2015-11-21 20:30 UTC (permalink / raw)
  To: Steve Kargl; +Cc: fortran, gcc-patches

‘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

^ permalink raw reply	[flat|nested] 20+ messages in thread

end of thread, other threads:[~2015-12-31 16:13 UTC | newest]

Thread overview: 20+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2015-11-20  0:58 [PATCH] (Partial) Implementation of simplificaiton of CSHIFT Steve Kargl
2015-11-20  1:31 ` Steve Kargl
2015-11-20  3:16   ` Steve Kargl
2015-11-20 20:09 ` Steve Kargl
2015-11-21 11:45   ` Paul Richard Thomas
2015-11-21 18:07     ` Steve Kargl
2015-11-21 18:12       ` H.J. Lu
2015-11-21 18:21         ` Steve Kargl
2015-11-21 19:26           ` H.J. Lu
2015-11-21 20:00             ` Steve Kargl
2015-11-21 20:22               ` Steve Kargl
2015-11-21 22:39                 ` Eric Botcazou
2015-11-22  1:53                   ` Steve Kargl
2015-12-31 14:00                 ` Gerald Pfeifer
2015-12-31 16:13                   ` Steve Kargl
2015-11-21 20:30 Dominique d'Humières
2015-11-21 21:50 ` Steve Kargl
2015-11-22  9:14   ` Dominique d'Humières
2015-11-22 16:57     ` Steve Kargl
2015-11-22 18:41       ` Dominique d'Humières

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).