row ops in ST

3 files changed, 92 insertions, 3 deletions

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index dcce1c5..e42889d 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -4,6 +4,12 @@
 #include <math.h>
 #include <time.h>
 #include <inttypes.h>
+#include <complex.h>
+
+typedef double complex TCD;
+typedef float  complex TCF;
+
+#undef complex
 
 #include "lapack-aux.h"
 
@@ -46,6 +52,10 @@
 #define NODEFPOS 2006
 #define NOSPRTD  2007
 
+inline int mod (int a, int b);
+
+inline int64_t mod_l (int64_t a, int64_t b);
+
 //---------------------------------------
 void asm_finit() {
 #ifdef i386
@@ -1310,6 +1320,83 @@ int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) {
 int multiplyI(int     m, KOIMAT(a), KOIMAT(b), OIMAT(r)) MULT_IMP
 int multiplyL(int64_t m, KOLMAT(a), KOLMAT(b), OLMAT(r)) MULT_IMP
 
+/////////////////////////////// inplace row ops ////////////////////////////////
+
+#define AXPY_IMP {                    \
+    int j;                            \
+    for(j=j1; j<=j2; j++) {           \
+        AT(r,i2,j) += a*AT(r,i1,j);   \
+    } OK }
+
+#define AXPY_MOD_IMP(M) {                                      \
+    int j;                                                     \
+    for(j=j1; j<=j2; j++) {                                    \
+        AT(r,i2,j) = M(AT(r,i2,j) + M(a*AT(r,i1,j), m) , m);   \
+    } OK }
+
+
+#define SCAL_IMP {                    \
+    int i,j;                          \
+    for(i=i1; i<=i2; i++) {           \
+        for(j=j1; j<=j2; j++) {       \
+            AT(r,i,j) = a*AT(r,i,j);  \
+            }                         \
+    } OK }
+
+#define SCAL_MOD_IMP(M) {                   \
+    int i,j;                                \
+    for(i=i1; i<=i2; i++) {                 \
+        for(j=j1; j<=j2; j++) {             \
+            AT(r,i,j) = M(a*AT(r,i,j) , m); \
+            }                               \
+    } OK }
+
+
+#define SWAP_IMP(T)   {               \
+    T aux;                            \
+    int k;                            \
+    if (i1 != i2) {                   \
+        for (k=j1; k<=j2; k++) {      \
+            aux = AT(r,i1,k);         \
+            AT(r,i1,k) = AT(r,i2,k);  \
+            AT(r,i2,k) = aux;         \
+        }                             \
+    } OK }
+
+
+#define ROWOP_IMP(T) {                \
+    T a = *pa;                        \
+    switch(code) {                    \
+        case 0:  AXPY_IMP             \
+        case 1:  SCAL_IMP             \
+        case 2:  SWAP_IMP(T)          \
+        default: ERROR(BAD_CODE);     \
+    }                                 \
+}
+
+#define ROWOP_MOD_IMP(T,M) {          \
+    T a = *pa;                        \
+    switch(code) {                    \
+        case 0:  AXPY_MOD_IMP(M)      \
+        case 1:  SCAL_MOD_IMP(M)      \
+        case 2:  SWAP_IMP(T)          \
+        default: ERROR(BAD_CODE);     \
+    }                                 \
+}
+
+
+#define ROWOP(T) int rowop_##T(int code, T* pa, int i1, int i2, int j1, int j2, MATG(T,r)) ROWOP_IMP(T)
+
+#define ROWOP_MOD(T,M) int rowop_mod_##T(T m, int code, T* pa, int i1, int i2, int j1, int j2, MATG(T,r)) ROWOP_MOD_IMP(T,M)
+
+ROWOP(double)
+ROWOP(float)
+ROWOP(TCD)
+ROWOP(TCF)
+ROWOP(int32_t)
+ROWOP(int64_t)
+ROWOP_MOD(int32_t,mod)
+ROWOP_MOD(int64_t,mod_l)
 
 ////////////////// sparse matrix-product ///////////////////////////////////////
 
diff --git a/packages/base/src/Internal/C/lapack-aux.h b/packages/base/src/Internal/C/lapack-aux.h
index 1549bb5..e4d95bc 100644
--- a/packages/base/src/Internal/C/lapack-aux.h
+++ b/packages/base/src/Internal/C/lapack-aux.h
@@ -59,6 +59,7 @@ typedef short ftnlen;
 #define OQMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, complex* A##p
 #define OCMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, doublecomplex* A##p
 
+#define MATG(T,A) int A##r, int A##c, int A##Xr, int A##Xc, T* A##p
 
 #define KIVEC(A) int A##n, const int*A##p
 #define KLVEC(A) int A##n, const int64_t*A##p
diff --git a/packages/base/src/Internal/C/vector-aux.c b/packages/base/src/Internal/C/vector-aux.c
index c161556..5528a9d 100644
--- a/packages/base/src/Internal/C/vector-aux.c
+++ b/packages/base/src/Internal/C/vector-aux.c
@@ -716,6 +716,7 @@ int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) {
     }
 }
 
+inline
 int mod (int a, int b) {
     int m = a % b;
     if (b>0) {
@@ -741,7 +742,7 @@ int mapValI(int code, int* pval, KIVEC(x), IVEC(r)) {
     }
 }
 
-
+inline
 int64_t mod_l (int64_t a, int64_t b) {
     int64_t m = a % b;
     if (b>0) {
@@ -1230,7 +1231,7 @@ int round_vector_i(KDVEC(v),IVEC(r)) {
 int mod_vector(int m, KIVEC(v), IVEC(r)) {
     int k;
     for(k=0; k<vn; k++) {
-        rp[k] = vp[k] % m;
+        rp[k] = mod(vp[k],m);
     }
     OK
 }
@@ -1266,7 +1267,7 @@ int round_vector_l(KDVEC(v),LVEC(r)) {
 int mod_vector_l(int64_t m, KLVEC(v), LVEC(r)) {
     int k;
     for(k=0; k<vn; k++) {
-        rp[k] = vp[k] % m;
+        rp[k] = mod_l(vp[k],m);
     }
     OK
 }