summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--packages/base/src/Internal/C/lapack-aux.c291
-rw-r--r--packages/base/src/Internal/C/vector-aux.c286
2 files changed, 291 insertions, 286 deletions
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 1402050..36d1e57 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -1300,6 +1300,35 @@ int multiplyI(KOIMAT(a), KOIMAT(b), OIMAT(r)) {
1300 OK 1300 OK
1301} 1301}
1302 1302
1303
1304////////////////// sparse matrix-product ///////////////////////////////////////
1305
1306
1307int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
1308 int r, c;
1309 for (r = 0; r < rowsn - 1; r++) {
1310 rp[r] = 0;
1311 for (c = rowsp[r]; c < rowsp[r+1]; c++) {
1312 rp[r] += valsp[c-1] * xp[colsp[c-1]-1];
1313 }
1314 }
1315 OK
1316}
1317
1318int smTXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
1319 int r,c;
1320 for (c = 0; c < rn; c++) {
1321 rp[c] = 0;
1322 }
1323 for (r = 0; r < rowsn - 1; r++) {
1324 for (c = rowsp[r]; c < rowsp[r+1]; c++) {
1325 rp[colsp[c-1]-1] += valsp[c-1] * xp[r];
1326 }
1327 }
1328 OK
1329}
1330
1331
1303//////////////////// transpose ///////////////////////// 1332//////////////////// transpose /////////////////////////
1304 1333
1305int transF(KFMAT(x),FMAT(t)) { 1334int transF(KFMAT(x),FMAT(t)) {
@@ -1376,248 +1405,6 @@ int transI(KIMAT(x),IMAT(t)) {
1376} 1405}
1377 1406
1378 1407
1379//////////////////// constant /////////////////////////
1380
1381int constantF(float * pval, FVEC(r)) {
1382 DEBUGMSG("constantF")
1383 int k;
1384 double val = *pval;
1385 for(k=0;k<rn;k++) {
1386 rp[k]=val;
1387 }
1388 OK
1389}
1390
1391int constantR(double * pval, DVEC(r)) {
1392 DEBUGMSG("constantR")
1393 int k;
1394 double val = *pval;
1395 for(k=0;k<rn;k++) {
1396 rp[k]=val;
1397 }
1398 OK
1399}
1400
1401int constantQ(complex* pval, QVEC(r)) {
1402 DEBUGMSG("constantQ")
1403 int k;
1404 complex val = *pval;
1405 for(k=0;k<rn;k++) {
1406 rp[k]=val;
1407 }
1408 OK
1409}
1410
1411int constantC(doublecomplex* pval, CVEC(r)) {
1412 DEBUGMSG("constantC")
1413 int k;
1414 doublecomplex val = *pval;
1415 for(k=0;k<rn;k++) {
1416 rp[k]=val;
1417 }
1418 OK
1419}
1420
1421int constantP(void* pval, PVEC(r)) {
1422 DEBUGMSG("constantP")
1423 int k;
1424 for(k=0;k<rn;k++) {
1425 memcpy(rp+k*rs,pval,rs);
1426 }
1427 OK
1428}
1429
1430
1431int constantI(int * pval, IVEC(r)) {
1432 DEBUGMSG("constantI")
1433 int k;
1434 int val = *pval;
1435 for(k=0;k<rn;k++) {
1436 rp[k]=val;
1437 }
1438 OK
1439}
1440
1441
1442//////////////////// float-double conversion /////////////////////////
1443
1444int float2double(FVEC(x),DVEC(y)) {
1445 DEBUGMSG("float2double")
1446 int k;
1447 for(k=0;k<xn;k++) {
1448 yp[k]=xp[k];
1449 }
1450 OK
1451}
1452
1453int float2int(KFVEC(x),IVEC(y)) {
1454 DEBUGMSG("float2int")
1455 int k;
1456 for(k=0;k<xn;k++) {
1457 yp[k]=xp[k];
1458 }
1459 OK
1460}
1461
1462
1463int double2float(DVEC(x),FVEC(y)) {
1464 DEBUGMSG("double2float")
1465 int k;
1466 for(k=0;k<xn;k++) {
1467 yp[k]=xp[k];
1468 }
1469 OK
1470}
1471
1472
1473int double2int(KDVEC(x),IVEC(y)) {
1474 DEBUGMSG("double2int")
1475 int k;
1476 for(k=0;k<xn;k++) {
1477 yp[k]=xp[k];
1478 }
1479 OK
1480}
1481
1482
1483int int2float(KIVEC(x),FVEC(y)) {
1484 DEBUGMSG("int2float")
1485 int k;
1486 for(k=0;k<xn;k++) {
1487 yp[k]=xp[k];
1488 }
1489 OK
1490}
1491
1492
1493int int2double(KIVEC(x),DVEC(y)) {
1494 DEBUGMSG("int2double")
1495 int k;
1496 for(k=0;k<xn;k++) {
1497 yp[k]=xp[k];
1498 }
1499 OK
1500}
1501
1502
1503//////////////////// conjugate /////////////////////////
1504
1505int conjugateQ(KQVEC(x),QVEC(t)) {
1506 REQUIRES(xn==tn,BAD_SIZE);
1507 DEBUGMSG("conjugateQ");
1508 int k;
1509 for(k=0;k<xn;k++) {
1510 tp[k].r = xp[k].r;
1511 tp[k].i = -xp[k].i;
1512 }
1513 OK
1514}
1515
1516int conjugateC(KCVEC(x),CVEC(t)) {
1517 REQUIRES(xn==tn,BAD_SIZE);
1518 DEBUGMSG("conjugateC");
1519 int k;
1520 for(k=0;k<xn;k++) {
1521 tp[k].r = xp[k].r;
1522 tp[k].i = -xp[k].i;
1523 }
1524 OK
1525}
1526
1527//////////////////// step /////////////////////////
1528
1529#define STEP_IMP \
1530 int k; \
1531 for(k=0;k<xn;k++) { \
1532 yp[k]=xp[k]>0; \
1533 } \
1534 OK
1535
1536int stepF(KFVEC(x),FVEC(y)) {
1537 STEP_IMP
1538}
1539
1540int stepD(KDVEC(x),DVEC(y)) {
1541 STEP_IMP
1542}
1543
1544int stepI(KIVEC(x),IVEC(y)) {
1545 STEP_IMP
1546}
1547
1548//////////////////// cond /////////////////////////
1549
1550#define COMPARE_IMP \
1551 REQUIRES(xn==yn && xn==rn ,BAD_SIZE); \
1552 int k; \
1553 for(k=0;k<xn;k++) { \
1554 rp[k] = xp[k]<yp[k]?-1:(xp[k]>yp[k]?1:0); \
1555 } \
1556 OK
1557
1558
1559int compareF(KFVEC(x),KFVEC(y),IVEC(r)) {
1560 COMPARE_IMP
1561}
1562
1563int compareD(KDVEC(x),KDVEC(y),IVEC(r)) {
1564 COMPARE_IMP
1565}
1566
1567int compareI(KIVEC(x),KIVEC(y),IVEC(r)) {
1568 COMPARE_IMP
1569}
1570
1571
1572#define COND_IMP \
1573 REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); \
1574 int k; \
1575 for(k=0;k<xn;k++) { \
1576 rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]); \
1577 } \
1578 OK
1579
1580int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) {
1581 COND_IMP
1582}
1583
1584int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {
1585 COND_IMP
1586}
1587
1588int condI(KIVEC(x),KIVEC(y),KIVEC(lt),KIVEC(eq),KIVEC(gt),IVEC(r)) {
1589 COND_IMP
1590}
1591
1592
1593#define CHOOSE_IMP \
1594 REQUIRES(condn==ltn && ltn==eqn && ltn==gtn && ltn==rn ,BAD_SIZE); \
1595 int k; \
1596 for(k=0;k<condn;k++) { \
1597 rp[k] = condp[k]<0?ltp[k]:(condp[k]>0?gtp[k]:eqp[k]); \
1598 } \
1599 OK
1600
1601int chooseF(KIVEC(cond),KFVEC(lt),KFVEC(eq),KFVEC(gt),FVEC(r)) {
1602 CHOOSE_IMP
1603}
1604
1605int chooseD(KIVEC(cond),KDVEC(lt),KDVEC(eq),KDVEC(gt),DVEC(r)) {
1606 CHOOSE_IMP
1607}
1608
1609int chooseI(KIVEC(cond),KIVEC(lt),KIVEC(eq),KIVEC(gt),IVEC(r)) {
1610 CHOOSE_IMP
1611}
1612
1613int chooseC(KIVEC(cond),KCVEC(lt),KCVEC(eq),KCVEC(gt),CVEC(r)) {
1614 CHOOSE_IMP
1615}
1616
1617int chooseQ(KIVEC(cond),KQVEC(lt),KQVEC(eq),KQVEC(gt),QVEC(r)) {
1618 CHOOSE_IMP
1619}
1620
1621//////////////////////// extract ///////////////////////////////// 1408//////////////////////// extract /////////////////////////////////
1622 1409
1623#define EXTRACT_IMP \ 1410#define EXTRACT_IMP \
@@ -1684,3 +1471,23 @@ int remapQ(KOIMAT(i), KOIMAT(j), KOQMAT(m), OQMAT(r)) {
1684 REMAP_IMP 1471 REMAP_IMP
1685} 1472}
1686 1473
1474////////////////////////////////////////////////////////////////////////////////
1475
1476int saveMatrix(char * file, char * format, KDMAT(a)){
1477 FILE * fp;
1478 fp = fopen (file, "w");
1479 int r, c;
1480 for (r=0;r<ar; r++) {
1481 for (c=0; c<ac; c++) {
1482 fprintf(fp,format,ap[r*ac+c]);
1483 if (c<ac-1) {
1484 fprintf(fp," ");
1485 } else {
1486 fprintf(fp,"\n");
1487 }
1488 }
1489 }
1490 fclose(fp);
1491 OK
1492}
1493
diff --git a/packages/base/src/Internal/C/vector-aux.c b/packages/base/src/Internal/C/vector-aux.c
index 5662697..a260f1e 100644
--- a/packages/base/src/Internal/C/vector-aux.c
+++ b/packages/base/src/Internal/C/vector-aux.c
@@ -809,24 +809,6 @@ int vectorScan(char * file, int* n, double**pp){
809 OK 809 OK
810} 810}
811 811
812int saveMatrix(char * file, char * format, KDMAT(a)){
813 FILE * fp;
814 fp = fopen (file, "w");
815 int r, c;
816 for (r=0;r<ar; r++) {
817 for (c=0; c<ac; c++) {
818 fprintf(fp,format,ap[r*ac+c]);
819 if (c<ac-1) {
820 fprintf(fp," ");
821 } else {
822 fprintf(fp,"\n");
823 }
824 }
825 }
826 fclose(fp);
827 OK
828}
829
830//////////////////////////////////////////////////////////////////////////////// 812////////////////////////////////////////////////////////////////////////////////
831 813
832#if defined (__APPLE__) || (__FreeBSD__) 814#if defined (__APPLE__) || (__FreeBSD__)
@@ -975,32 +957,6 @@ int random_vector(unsigned int seed, int code, DVEC(r)) {
975 957
976//////////////////////////////////////////////////////////////////////////////// 958////////////////////////////////////////////////////////////////////////////////
977 959
978int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
979 int r, c;
980 for (r = 0; r < rowsn - 1; r++) {
981 rp[r] = 0;
982 for (c = rowsp[r]; c < rowsp[r+1]; c++) {
983 rp[r] += valsp[c-1] * xp[colsp[c-1]-1];
984 }
985 }
986 OK
987}
988
989int smTXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
990 int r,c;
991 for (c = 0; c < rn; c++) {
992 rp[c] = 0;
993 }
994 for (r = 0; r < rowsn - 1; r++) {
995 for (c = rowsp[r]; c < rowsp[r+1]; c++) {
996 rp[colsp[c-1]-1] += valsp[c-1] * xp[r];
997 }
998 }
999 OK
1000}
1001
1002////////////////////////////////////////////////////////////////////////////////
1003
1004int 960int
1005compare_doubles (const void *a, const void *b) { 961compare_doubles (const void *a, const void *b) {
1006 return *(double*)a > *(double*)b; 962 return *(double*)a > *(double*)b;
@@ -1132,3 +1088,245 @@ int range_vector(IVEC(r)) {
1132 OK 1088 OK
1133} 1089}
1134 1090
1091//////////////////// constant /////////////////////////
1092
1093int constantF(float * pval, FVEC(r)) {
1094 DEBUGMSG("constantF")
1095 int k;
1096 double val = *pval;
1097 for(k=0;k<rn;k++) {
1098 rp[k]=val;
1099 }
1100 OK
1101}
1102
1103int constantR(double * pval, DVEC(r)) {
1104 DEBUGMSG("constantR")
1105 int k;
1106 double val = *pval;
1107 for(k=0;k<rn;k++) {
1108 rp[k]=val;
1109 }
1110 OK
1111}
1112
1113int constantQ(complex* pval, QVEC(r)) {
1114 DEBUGMSG("constantQ")
1115 int k;
1116 complex val = *pval;
1117 for(k=0;k<rn;k++) {
1118 rp[k]=val;
1119 }
1120 OK
1121}
1122
1123int constantC(doublecomplex* pval, CVEC(r)) {
1124 DEBUGMSG("constantC")
1125 int k;
1126 doublecomplex val = *pval;
1127 for(k=0;k<rn;k++) {
1128 rp[k]=val;
1129 }
1130 OK
1131}
1132
1133int constantP(void* pval, PVEC(r)) {
1134 DEBUGMSG("constantP")
1135 int k;
1136 for(k=0;k<rn;k++) {
1137 memcpy(rp+k*rs,pval,rs);
1138 }
1139 OK
1140}
1141
1142
1143int constantI(int * pval, IVEC(r)) {
1144 DEBUGMSG("constantI")
1145 int k;
1146 int val = *pval;
1147 for(k=0;k<rn;k++) {
1148 rp[k]=val;
1149 }
1150 OK
1151}
1152
1153
1154//////////////////// float-double conversion /////////////////////////
1155
1156int float2double(FVEC(x),DVEC(y)) {
1157 DEBUGMSG("float2double")
1158 int k;
1159 for(k=0;k<xn;k++) {
1160 yp[k]=xp[k];
1161 }
1162 OK
1163}
1164
1165int float2int(KFVEC(x),IVEC(y)) {
1166 DEBUGMSG("float2int")
1167 int k;
1168 for(k=0;k<xn;k++) {
1169 yp[k]=xp[k];
1170 }
1171 OK
1172}
1173
1174
1175int double2float(DVEC(x),FVEC(y)) {
1176 DEBUGMSG("double2float")
1177 int k;
1178 for(k=0;k<xn;k++) {
1179 yp[k]=xp[k];
1180 }
1181 OK
1182}
1183
1184
1185int double2int(KDVEC(x),IVEC(y)) {
1186 DEBUGMSG("double2int")
1187 int k;
1188 for(k=0;k<xn;k++) {
1189 yp[k]=xp[k];
1190 }
1191 OK
1192}
1193
1194
1195int int2float(KIVEC(x),FVEC(y)) {
1196 DEBUGMSG("int2float")
1197 int k;
1198 for(k=0;k<xn;k++) {
1199 yp[k]=xp[k];
1200 }
1201 OK
1202}
1203
1204
1205int int2double(KIVEC(x),DVEC(y)) {
1206 DEBUGMSG("int2double")
1207 int k;
1208 for(k=0;k<xn;k++) {
1209 yp[k]=xp[k];
1210 }
1211 OK
1212}
1213
1214
1215//////////////////// conjugate /////////////////////////
1216
1217int conjugateQ(KQVEC(x),QVEC(t)) {
1218 REQUIRES(xn==tn,BAD_SIZE);
1219 DEBUGMSG("conjugateQ");
1220 int k;
1221 for(k=0;k<xn;k++) {
1222 tp[k].r = xp[k].r;
1223 tp[k].i = -xp[k].i;
1224 }
1225 OK
1226}
1227
1228int conjugateC(KCVEC(x),CVEC(t)) {
1229 REQUIRES(xn==tn,BAD_SIZE);
1230 DEBUGMSG("conjugateC");
1231 int k;
1232 for(k=0;k<xn;k++) {
1233 tp[k].r = xp[k].r;
1234 tp[k].i = -xp[k].i;
1235 }
1236 OK
1237}
1238
1239//////////////////// step /////////////////////////
1240
1241#define STEP_IMP \
1242 int k; \
1243 for(k=0;k<xn;k++) { \
1244 yp[k]=xp[k]>0; \
1245 } \
1246 OK
1247
1248int stepF(KFVEC(x),FVEC(y)) {
1249 STEP_IMP
1250}
1251
1252int stepD(KDVEC(x),DVEC(y)) {
1253 STEP_IMP
1254}
1255
1256int stepI(KIVEC(x),IVEC(y)) {
1257 STEP_IMP
1258}
1259
1260//////////////////// cond /////////////////////////
1261
1262#define COMPARE_IMP \
1263 REQUIRES(xn==yn && xn==rn ,BAD_SIZE); \
1264 int k; \
1265 for(k=0;k<xn;k++) { \
1266 rp[k] = xp[k]<yp[k]?-1:(xp[k]>yp[k]?1:0); \
1267 } \
1268 OK
1269
1270
1271int compareF(KFVEC(x),KFVEC(y),IVEC(r)) {
1272 COMPARE_IMP
1273}
1274
1275int compareD(KDVEC(x),KDVEC(y),IVEC(r)) {
1276 COMPARE_IMP
1277}
1278
1279int compareI(KIVEC(x),KIVEC(y),IVEC(r)) {
1280 COMPARE_IMP
1281}
1282
1283
1284#define COND_IMP \
1285 REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); \
1286 int k; \
1287 for(k=0;k<xn;k++) { \
1288 rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]); \
1289 } \
1290 OK
1291
1292int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) {
1293 COND_IMP
1294}
1295
1296int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {
1297 COND_IMP
1298}
1299
1300int condI(KIVEC(x),KIVEC(y),KIVEC(lt),KIVEC(eq),KIVEC(gt),IVEC(r)) {
1301 COND_IMP
1302}
1303
1304
1305#define CHOOSE_IMP \
1306 REQUIRES(condn==ltn && ltn==eqn && ltn==gtn && ltn==rn ,BAD_SIZE); \
1307 int k; \
1308 for(k=0;k<condn;k++) { \
1309 rp[k] = condp[k]<0?ltp[k]:(condp[k]>0?gtp[k]:eqp[k]); \
1310 } \
1311 OK
1312
1313int chooseF(KIVEC(cond),KFVEC(lt),KFVEC(eq),KFVEC(gt),FVEC(r)) {
1314 CHOOSE_IMP
1315}
1316
1317int chooseD(KIVEC(cond),KDVEC(lt),KDVEC(eq),KDVEC(gt),DVEC(r)) {
1318 CHOOSE_IMP
1319}
1320
1321int chooseI(KIVEC(cond),KIVEC(lt),KIVEC(eq),KIVEC(gt),IVEC(r)) {
1322 CHOOSE_IMP
1323}
1324
1325int chooseC(KIVEC(cond),KCVEC(lt),KCVEC(eq),KCVEC(gt),CVEC(r)) {
1326 CHOOSE_IMP
1327}
1328
1329int chooseQ(KIVEC(cond),KQVEC(lt),KQVEC(eq),KQVEC(gt),QVEC(r)) {
1330 CHOOSE_IMP
1331}
1332
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 08:26:06 +0000