[matgen] Cleanup power with integer exponent

Add header.
Include header instead of forward declarations.
Rename function to no longer follow Fortran
naming scheme.
Call by value instead of reference.
Move float and double functions into a single file.
Add Doxygen documentation.

Author: gruenich

DOC/modules.txt

@@ -56,3 +56,9 @@
 
  These tests start with the letter z.
  */
+
+/*!
+ \defgroup TestingMatgen Matrix generation
+ \ingroup Testing
+ \brief Matrix generation for tests.
+ */

TESTING/MATGEN/CMakeLists.txt

@@ -30,7 +30,7 @@ if(enable_single)
     slarot.c
     slatm2.c
     slatm3.c
-    pow_ri.c
+    powi.c
   )
 endif()
 
@@ -48,7 +48,7 @@ if(enable_double)
     dlarot.c
     dlatm2.c
     dlatm3.c
-    pow_di.c
+    powi.c
   )
 endif()
 
@@ -72,7 +72,7 @@ if(enable_complex)
     claghe.c
     clarnd.c
 #    cdotc.c
-    pow_ri.c
+    powi.c
 )
 endif()
 
@@ -96,7 +96,7 @@ if(enable_complex16)
     zlaghe.c
     zlarnd.c
 #    zdotc.c
-    pow_di.c
+    powi.c
   )
 endif()
 

TESTING/MATGEN/Makefile

@@ -33,12 +33,12 @@ include ../../make.inc
 SCATGEN = slatm1.o slaran.o slarnd.o slaruv.o slabad.o slarnv.o
 SLASRC  = slatb4.o slaset.o slartg.o
 SMATGEN = slatms.o slagge.o slagsy.o slarge.o slaror.o slarot.o slatm2.o slatm3.o
-SINTRINSIC  = pow_ri.o
+SINTRINSIC  = powi.o
 
 DZATGEN = dlatm1.o dlaran.o dlarnd.o dlaruv.o dlabad.o dlarnv.o
 DLASRC  = dlatb4.o dlaset.o dlartg.o
 DMATGEN = dlatms.o dlagge.o dlagsy.o dlarge.o dlaror.o dlarot.o dlatm2.o dlatm3.o
-DINTRINSIC  = pow_di.o
+DINTRINSIC  = powi.o
 
 CLASRC  = clatb4.o claset.o clartg.o clarnv.o clacgv.o csymv.o
 CMATGEN = clatms.o clagge.o clagsy.o clarge.o claror.o clarot.o clatm2.o clatm3.o \

TESTING/MATGEN/dlartg.c

@@ -1,6 +1,8 @@
+#include "../../SRC/slu_ddefs.h"
+#include "powi.h"
+
 #include <stdbool.h>
 #include <math.h>
-#include "../../SRC/slu_ddefs.h"
 
 /* Subroutine */ int dlartg_slu(double *f, double *g, double *cs, double *sn, double *r)
 {
@@ -53,8 +55,7 @@
     /* System generated locals */
     int i__1;
     double d__1, d__2;
-    /* Builtin functions */
-    double pow_di(double *, int *);
+
     /* Local variables */
     static int i;
     static double scale;
@@ -71,7 +72,7 @@
 	eps = dmach("E");
 	d__1 = dmach("B");
 	i__1 = (int) (log(safmin / eps) / log(dmach("B")) / 2.);
-	safmn2 = pow_di(&d__1, &i__1);
+	safmn2 = powi(d__1, i__1);
 	safmx2 = 1. / safmn2;
     }
     if (*g == 0.) {

TESTING/MATGEN/dlatm1.c

@@ -1,5 +1,6 @@
 /*  -- translated by f2c (version 19940927).
 */
+#include "powi.h"
 
 #include <math.h>
 #include <stdlib.h>
@@ -12,9 +13,6 @@
     int i__1, i__2;
     double d__1;
 
-    /* Builtin functions */
-    double pow_di(double *, int *);
-
     /* Local variables */
     static double temp;
     static int i;
@@ -189,7 +187,7 @@
 	    i__1 = *n;
 	    for (i = 2; i <= i__1; ++i) {
 		i__2 = i - 1;
-		d[i] = pow_di(&alpha, &i__2);
+		d[i] = powi(alpha, i__2);
 /* L60: */
 	    }
 	}

TESTING/MATGEN/pow_di.c

@@ -0,0 +1,68 @@
+/*! \file
+ * \ingroup TestingMatgen
+ * \brief Power with integer exponent
+ */
+
+/*!
+ * Power with integer exponent
+ *
+ * Implemented as Exponentiation by Squaring for speed.
+ *
+ * \param[in] base base as double value
+ * \param[in] exp  exponent as integer value
+ * \return Computes the value of \p base raised to the power \p exp
+ */
+double powi(double base, int exp)
+{
+    double x = base;
+    int n = exp;
+    double pow = 1.0;
+
+    if (n != 0) {
+        if(n < 0) {
+            n = -n;
+            x = 1 / x;
+        }
+        for( ; ; ) {
+            if(n & 01)
+                pow *= x;
+            if(n >>= 1)
+                x *= x;
+            else
+                break;
+        }
+    }
+    return pow;
+}
+
+/*!
+ * Power with integer exponent
+ *
+ * Implemented as Exponentiation by Squaring for speed.
+ *
+ * \param[in] base base as float value
+ * \param[in] exp  exponent as integer value
+ * \return Computes the value of \p base raised to the power \p exp
+ */
+float powif(float base, int exp)
+{
+    float x = base;
+    int n = exp;
+    float pow = 1.0;
+
+    if (n != 0)  {
+        if (n < 0) {
+            n = -n;
+            x = 1 / x;
+        }
+        for ( ; ; ) {
+            if (n & 01)
+                pow *= x;
+            if (n >>= 1)
+                x *= x;
+            else
+                break;
+        }
+    }
+    return pow;
+}

TESTING/MATGEN/pow_ri.c

@@ -0,0 +1,3 @@
+double powi(double base, int exp);
+
+float powif(float base, int exp);

TESTING/MATGEN/powi.c

@@ -1,3 +1,5 @@
+#include "powi.h"
+
 #include <math.h>
 #include <stdbool.h>
 
@@ -52,8 +54,7 @@
     /* System generated locals */
     int i__1;
     float r__1, r__2;
-    /* Builtin functions */
-    double pow_ri(float *, int *);
+
     /* Local variables */
     static int i;
     static float scale;
@@ -69,7 +70,7 @@
 	eps = smach("E");
 	r__1 = smach("B");
 	i__1 = (int) (log(safmin / eps) / log(smach("B")) / 2.f);
-	safmn2 = pow_ri(&r__1, &i__1);
+	safmn2 = powif(r__1, i__1);
 	safmx2 = 1.f / safmn2;
     }
     if (*g == 0.f) {

TESTING/MATGEN/powi.h

@@ -1,5 +1,6 @@
 /*  -- translated by f2c (version 19940927).
 */
+#include "powi.h"
 
 #include <math.h>
 #include <stdlib.h>
@@ -11,9 +12,6 @@
     int i__1, i__2;
     double d__1, d__2;
 
-    /* Builtin functions */
-    double pow_ri(float *, int *);
-
     /* Local variables */
     static float temp;
     static int i;
@@ -190,7 +188,7 @@
 	    i__1 = *n;
 	    for (i = 2; i <= i__1; ++i) {
 		i__2 = i - 1;
-		d[i] = pow_ri(&alpha, &i__2);
+		d[i] = powif(alpha, i__2);
 /* L60: */
 	    }
 	}