diff options
Diffstat (limited to 'packages/base/src/Internal/C')
-rw-r--r-- | packages/base/src/Internal/C/lapack-aux.c | 87 | ||||
-rw-r--r-- | packages/base/src/Internal/C/lapack-aux.h | 1 | ||||
-rw-r--r-- | packages/base/src/Internal/C/vector-aux.c | 7 |
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 @@ | |||
4 | #include <math.h> | 4 | #include <math.h> |
5 | #include <time.h> | 5 | #include <time.h> |
6 | #include <inttypes.h> | 6 | #include <inttypes.h> |
7 | #include <complex.h> | ||
8 | |||
9 | typedef double complex TCD; | ||
10 | typedef float complex TCF; | ||
11 | |||
12 | #undef complex | ||
7 | 13 | ||
8 | #include "lapack-aux.h" | 14 | #include "lapack-aux.h" |
9 | 15 | ||
@@ -46,6 +52,10 @@ | |||
46 | #define NODEFPOS 2006 | 52 | #define NODEFPOS 2006 |
47 | #define NOSPRTD 2007 | 53 | #define NOSPRTD 2007 |
48 | 54 | ||
55 | inline int mod (int a, int b); | ||
56 | |||
57 | inline int64_t mod_l (int64_t a, int64_t b); | ||
58 | |||
49 | //--------------------------------------- | 59 | //--------------------------------------- |
50 | void asm_finit() { | 60 | void asm_finit() { |
51 | #ifdef i386 | 61 | #ifdef i386 |
@@ -1310,6 +1320,83 @@ int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) { | |||
1310 | int multiplyI(int m, KOIMAT(a), KOIMAT(b), OIMAT(r)) MULT_IMP | 1320 | int multiplyI(int m, KOIMAT(a), KOIMAT(b), OIMAT(r)) MULT_IMP |
1311 | int multiplyL(int64_t m, KOLMAT(a), KOLMAT(b), OLMAT(r)) MULT_IMP | 1321 | int multiplyL(int64_t m, KOLMAT(a), KOLMAT(b), OLMAT(r)) MULT_IMP |
1312 | 1322 | ||
1323 | /////////////////////////////// inplace row ops //////////////////////////////// | ||
1324 | |||
1325 | #define AXPY_IMP { \ | ||
1326 | int j; \ | ||
1327 | for(j=j1; j<=j2; j++) { \ | ||
1328 | AT(r,i2,j) += a*AT(r,i1,j); \ | ||
1329 | } OK } | ||
1330 | |||
1331 | #define AXPY_MOD_IMP(M) { \ | ||
1332 | int j; \ | ||
1333 | for(j=j1; j<=j2; j++) { \ | ||
1334 | AT(r,i2,j) = M(AT(r,i2,j) + M(a*AT(r,i1,j), m) , m); \ | ||
1335 | } OK } | ||
1336 | |||
1337 | |||
1338 | #define SCAL_IMP { \ | ||
1339 | int i,j; \ | ||
1340 | for(i=i1; i<=i2; i++) { \ | ||
1341 | for(j=j1; j<=j2; j++) { \ | ||
1342 | AT(r,i,j) = a*AT(r,i,j); \ | ||
1343 | } \ | ||
1344 | } OK } | ||
1345 | |||
1346 | #define SCAL_MOD_IMP(M) { \ | ||
1347 | int i,j; \ | ||
1348 | for(i=i1; i<=i2; i++) { \ | ||
1349 | for(j=j1; j<=j2; j++) { \ | ||
1350 | AT(r,i,j) = M(a*AT(r,i,j) , m); \ | ||
1351 | } \ | ||
1352 | } OK } | ||
1353 | |||
1354 | |||
1355 | #define SWAP_IMP(T) { \ | ||
1356 | T aux; \ | ||
1357 | int k; \ | ||
1358 | if (i1 != i2) { \ | ||
1359 | for (k=j1; k<=j2; k++) { \ | ||
1360 | aux = AT(r,i1,k); \ | ||
1361 | AT(r,i1,k) = AT(r,i2,k); \ | ||
1362 | AT(r,i2,k) = aux; \ | ||
1363 | } \ | ||
1364 | } OK } | ||
1365 | |||
1366 | |||
1367 | #define ROWOP_IMP(T) { \ | ||
1368 | T a = *pa; \ | ||
1369 | switch(code) { \ | ||
1370 | case 0: AXPY_IMP \ | ||
1371 | case 1: SCAL_IMP \ | ||
1372 | case 2: SWAP_IMP(T) \ | ||
1373 | default: ERROR(BAD_CODE); \ | ||
1374 | } \ | ||
1375 | } | ||
1376 | |||
1377 | #define ROWOP_MOD_IMP(T,M) { \ | ||
1378 | T a = *pa; \ | ||
1379 | switch(code) { \ | ||
1380 | case 0: AXPY_MOD_IMP(M) \ | ||
1381 | case 1: SCAL_MOD_IMP(M) \ | ||
1382 | case 2: SWAP_IMP(T) \ | ||
1383 | default: ERROR(BAD_CODE); \ | ||
1384 | } \ | ||
1385 | } | ||
1386 | |||
1387 | |||
1388 | #define ROWOP(T) int rowop_##T(int code, T* pa, int i1, int i2, int j1, int j2, MATG(T,r)) ROWOP_IMP(T) | ||
1389 | |||
1390 | #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) | ||
1391 | |||
1392 | ROWOP(double) | ||
1393 | ROWOP(float) | ||
1394 | ROWOP(TCD) | ||
1395 | ROWOP(TCF) | ||
1396 | ROWOP(int32_t) | ||
1397 | ROWOP(int64_t) | ||
1398 | ROWOP_MOD(int32_t,mod) | ||
1399 | ROWOP_MOD(int64_t,mod_l) | ||
1313 | 1400 | ||
1314 | ////////////////// sparse matrix-product /////////////////////////////////////// | 1401 | ////////////////// sparse matrix-product /////////////////////////////////////// |
1315 | 1402 | ||
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; | |||
59 | #define OQMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, complex* A##p | 59 | #define OQMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, complex* A##p |
60 | #define OCMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, doublecomplex* A##p | 60 | #define OCMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, doublecomplex* A##p |
61 | 61 | ||
62 | #define MATG(T,A) int A##r, int A##c, int A##Xr, int A##Xc, T* A##p | ||
62 | 63 | ||
63 | #define KIVEC(A) int A##n, const int*A##p | 64 | #define KIVEC(A) int A##n, const int*A##p |
64 | #define KLVEC(A) int A##n, const int64_t*A##p | 65 | #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)) { | |||
716 | } | 716 | } |
717 | } | 717 | } |
718 | 718 | ||
719 | inline | ||
719 | int mod (int a, int b) { | 720 | int mod (int a, int b) { |
720 | int m = a % b; | 721 | int m = a % b; |
721 | if (b>0) { | 722 | if (b>0) { |
@@ -741,7 +742,7 @@ int mapValI(int code, int* pval, KIVEC(x), IVEC(r)) { | |||
741 | } | 742 | } |
742 | } | 743 | } |
743 | 744 | ||
744 | 745 | inline | |
745 | int64_t mod_l (int64_t a, int64_t b) { | 746 | int64_t mod_l (int64_t a, int64_t b) { |
746 | int64_t m = a % b; | 747 | int64_t m = a % b; |
747 | if (b>0) { | 748 | if (b>0) { |
@@ -1230,7 +1231,7 @@ int round_vector_i(KDVEC(v),IVEC(r)) { | |||
1230 | int mod_vector(int m, KIVEC(v), IVEC(r)) { | 1231 | int mod_vector(int m, KIVEC(v), IVEC(r)) { |
1231 | int k; | 1232 | int k; |
1232 | for(k=0; k<vn; k++) { | 1233 | for(k=0; k<vn; k++) { |
1233 | rp[k] = vp[k] % m; | 1234 | rp[k] = mod(vp[k],m); |
1234 | } | 1235 | } |
1235 | OK | 1236 | OK |
1236 | } | 1237 | } |
@@ -1266,7 +1267,7 @@ int round_vector_l(KDVEC(v),LVEC(r)) { | |||
1266 | int mod_vector_l(int64_t m, KLVEC(v), LVEC(r)) { | 1267 | int mod_vector_l(int64_t m, KLVEC(v), LVEC(r)) { |
1267 | int k; | 1268 | int k; |
1268 | for(k=0; k<vn; k++) { | 1269 | for(k=0; k<vn; k++) { |
1269 | rp[k] = vp[k] % m; | 1270 | rp[k] = mod_l(vp[k],m); |
1270 | } | 1271 | } |
1271 | OK | 1272 | OK |
1272 | } | 1273 | } |