summaryrefslogtreecommitdiff
path: root/packages/base/src/Internal/C
diff options
context:
space:
mode:
authormaxc01 <xingchen92@gmail.com>2015-10-07 13:48:26 +0800
committermaxc01 <xingchen92@gmail.com>2015-10-07 13:48:26 +0800
commita61af756ddca4544de5e4969edc73131f4fccdd1 (patch)
tree2ac1755695a42d3964208e0029e74d446f5c3bd8 /packages/base/src/Internal/C
parent0840304af1564fa86a6006d648450372f301a6c8 (diff)
parentc84a485f148063f6d0c23f016fe348ec94fb6b19 (diff)
Merge pull request #1 from albertoruiz/master
sync from albetoruiz/hmatrix
Diffstat (limited to 'packages/base/src/Internal/C')
-rw-r--r--packages/base/src/Internal/C/lapack-aux.c1544
-rw-r--r--packages/base/src/Internal/C/lapack-aux.h111
-rw-r--r--packages/base/src/Internal/C/vector-aux.c1486
3 files changed, 3141 insertions, 0 deletions
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
new file mode 100644
index 0000000..ff7ad92
--- /dev/null
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -0,0 +1,1544 @@
1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4#include <math.h>
5#include <time.h>
6#include <inttypes.h>
7#include <complex.h>
8
9typedef double complex TCD;
10typedef float complex TCF;
11
12#undef complex
13
14#include "lapack-aux.h"
15
16#define MACRO(B) do {B} while (0)
17#define ERROR(CODE) MACRO(return CODE;)
18#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
19
20#define MIN(A,B) ((A)<(B)?(A):(B))
21#define MAX(A,B) ((A)>(B)?(A):(B))
22
23// #define DBGL
24
25#ifdef DBGL
26#define DEBUGMSG(M) printf("\nLAPACK "M"\n");
27#else
28#define DEBUGMSG(M)
29#endif
30
31#define OK return 0;
32
33// #ifdef DBGL
34// #define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL);
35// #define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;);
36// #else
37// #define DEBUGMSG(M)
38// #define OK return 0;
39// #endif
40
41
42#define INFOMAT(M) printf("%dx%d %d:%d\n",M##r,M##c,M##Xr,M##Xc);
43
44#define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \
45 for(q=0;q<M##r*M##c;q++) printf("%.1f ",M##p[q]); printf("\n");}
46
47#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
48
49#define BAD_SIZE 2000
50#define BAD_CODE 2001
51#define MEM 2002
52#define BAD_FILE 2003
53#define SINGULAR 2004
54#define NOCONVER 2005
55#define NODEFPOS 2006
56#define NOSPRTD 2007
57
58////////////////////////////////////////////////////////////////////////////////
59void asm_finit() {
60#ifdef i386
61
62// asm("finit");
63
64 static unsigned char buf[108];
65 asm("FSAVE %0":"=m" (buf));
66
67 #if FPUDEBUG
68 if(buf[8]!=255 || buf[9]!=255) { // print warning in red
69 printf("%c[;31mWarning: FPU TAG = %x %x\%c[0m\n",0x1B,buf[8],buf[9],0x1B);
70 }
71 #endif
72
73 #if NANDEBUG
74 asm("FRSTOR %0":"=m" (buf));
75 #endif
76
77#endif
78}
79
80#if NANDEBUG
81
82#define CHECKNANR(M,msg) \
83{ int k; \
84for(k=0; k<(M##r * M##c); k++) { \
85 if(M##p[k] != M##p[k]) { \
86 printf(msg); \
87 TRACEMAT(M) \
88 /*exit(1);*/ \
89 } \
90} \
91}
92
93#define CHECKNANC(M,msg) \
94{ int k; \
95for(k=0; k<(M##r * M##c); k++) { \
96 if( M##p[k].r != M##p[k].r \
97 || M##p[k].i != M##p[k].i) { \
98 printf(msg); \
99 /*exit(1);*/ \
100 } \
101} \
102}
103
104#else
105#define CHECKNANC(M,msg)
106#define CHECKNANR(M,msg)
107#endif
108
109
110////////////////////////////////////////////////////////////////////////////////
111//////////////////// real svd ///////////////////////////////////////////////////
112
113int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
114 doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
115 ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
116 integer *info);
117
118int svd_l_R(ODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
119 integer m = ar;
120 integer n = ac;
121 integer q = MIN(m,n);
122 REQUIRES(sn==q,BAD_SIZE);
123 REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE);
124 char* jobu = "A";
125 if (up==NULL) {
126 jobu = "N";
127 } else {
128 if (uc==q) {
129 jobu = "S";
130 }
131 }
132 REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE);
133 char* jobvt = "A";
134 integer ldvt = n;
135 if (vp==NULL) {
136 jobvt = "N";
137 } else {
138 if (vr==q) {
139 jobvt = "S";
140 ldvt = q;
141 }
142 }
143 DEBUGMSG("svd_l_R");
144 integer lwork = -1;
145 integer res;
146 // ask for optimal lwork
147 double ans;
148 dgesvd_ (jobu,jobvt,
149 &m,&n,ap,&m,
150 sp,
151 up,&m,
152 vp,&ldvt,
153 &ans, &lwork,
154 &res);
155 lwork = ceil(ans);
156 double * work = (double*)malloc(lwork*sizeof(double));
157 CHECK(!work,MEM);
158 dgesvd_ (jobu,jobvt,
159 &m,&n,ap,&m,
160 sp,
161 up,&m,
162 vp,&ldvt,
163 work, &lwork,
164 &res);
165 CHECK(res,res);
166 free(work);
167 OK
168}
169
170// (alternative version)
171
172int dgesdd_(char *jobz, integer *m, integer *n, doublereal *
173 a, integer *lda, doublereal *s, doublereal *u, integer *ldu,
174 doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
175 integer *iwork, integer *info);
176
177int svd_l_Rdd(ODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
178 integer m = ar;
179 integer n = ac;
180 integer q = MIN(m,n);
181 REQUIRES(sn==q,BAD_SIZE);
182 REQUIRES((up == NULL && vp == NULL)
183 || (ur==m && vc==n
184 && ((uc == q && vr == q)
185 || (uc == m && vc==n))),BAD_SIZE);
186 char* jobz = "A";
187 integer ldvt = n;
188 if (up==NULL) {
189 jobz = "N";
190 } else {
191 if (uc==q && vr == q) {
192 jobz = "S";
193 ldvt = q;
194 }
195 }
196 DEBUGMSG("svd_l_Rdd");
197 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
198 CHECK(!iwk,MEM);
199 integer lwk = -1;
200 integer res;
201 // ask for optimal lwk
202 double ans;
203 dgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res);
204 lwk = ans;
205 double * workv = (double*)malloc(lwk*sizeof(double));
206 CHECK(!workv,MEM);
207 dgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res);
208 CHECK(res,res);
209 free(iwk);
210 free(workv);
211 OK
212}
213
214//////////////////// complex svd ////////////////////////////////////
215
216int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
217 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
218 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
219 integer *lwork, doublereal *rwork, integer *info);
220
221int svd_l_C(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
222 integer m = ar;
223 integer n = ac;
224 integer q = MIN(m,n);
225 REQUIRES(sn==q,BAD_SIZE);
226 REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE);
227 char* jobu = "A";
228 if (up==NULL) {
229 jobu = "N";
230 } else {
231 if (uc==q) {
232 jobu = "S";
233 }
234 }
235 REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE);
236 char* jobvt = "A";
237 integer ldvt = n;
238 if (vp==NULL) {
239 jobvt = "N";
240 } else {
241 if (vr==q) {
242 jobvt = "S";
243 ldvt = q;
244 }
245 }DEBUGMSG("svd_l_C");
246
247 double *rwork = (double*) malloc(5*q*sizeof(double));
248 CHECK(!rwork,MEM);
249 integer lwork = -1;
250 integer res;
251 // ask for optimal lwork
252 doublecomplex ans;
253 zgesvd_ (jobu,jobvt,
254 &m,&n,ap,&m,
255 sp,
256 up,&m,
257 vp,&ldvt,
258 &ans, &lwork,
259 rwork,
260 &res);
261 lwork = ceil(ans.r);
262 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
263 CHECK(!work,MEM);
264 zgesvd_ (jobu,jobvt,
265 &m,&n,ap,&m,
266 sp,
267 up,&m,
268 vp,&ldvt,
269 work, &lwork,
270 rwork,
271 &res);
272 CHECK(res,res);
273 free(work);
274 free(rwork);
275 OK
276}
277
278int zgesdd_ (char *jobz, integer *m, integer *n,
279 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
280 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
281 integer *lwork, doublereal *rwork, integer* iwork, integer *info);
282
283int svd_l_Cdd(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
284 integer m = ar;
285 integer n = ac;
286 integer q = MIN(m,n);
287 REQUIRES(sn==q,BAD_SIZE);
288 REQUIRES((up == NULL && vp == NULL)
289 || (ur==m && vc==n
290 && ((uc == q && vr == q)
291 || (uc == m && vc==n))),BAD_SIZE);
292 char* jobz = "A";
293 integer ldvt = n;
294 if (up==NULL) {
295 jobz = "N";
296 } else {
297 if (uc==q && vr == q) {
298 jobz = "S";
299 ldvt = q;
300 }
301 }
302 DEBUGMSG("svd_l_Cdd");
303 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
304 CHECK(!iwk,MEM);
305 int lrwk;
306 if (0 && *jobz == 'N') {
307 lrwk = 5*q; // does not work, crash at free below
308 } else {
309 lrwk = 5*q*q + 7*q;
310 }
311 double *rwk = (double*)malloc(lrwk*sizeof(double));;
312 CHECK(!rwk,MEM);
313 integer lwk = -1;
314 integer res;
315 // ask for optimal lwk
316 doublecomplex ans;
317 zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
318 lwk = ans.r;
319 doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));
320 CHECK(!workv,MEM);
321 zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res);
322 CHECK(res,res);
323 free(workv);
324 free(rwk);
325 free(iwk);
326 OK
327}
328
329//////////////////// general complex eigensystem ////////////
330
331int zgeev_(char *jobvl, char *jobvr, integer *n,
332 doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl,
333 integer *ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work,
334 integer *lwork, doublereal *rwork, integer *info);
335
336int eig_l_C(OCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
337 integer n = ar;
338 REQUIRES(ac==n && sn==n, BAD_SIZE);
339 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
340 char jobvl = up==NULL?'N':'V';
341 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
342 char jobvr = vp==NULL?'N':'V';
343 DEBUGMSG("eig_l_C");
344 double *rwork = (double*) malloc(2*n*sizeof(double));
345 CHECK(!rwork,MEM);
346 integer lwork = -1;
347 integer res;
348 // ask for optimal lwork
349 doublecomplex ans;
350 zgeev_ (&jobvl,&jobvr,
351 &n,ap,&n,
352 sp,
353 up,&n,
354 vp,&n,
355 &ans, &lwork,
356 rwork,
357 &res);
358 lwork = ceil(ans.r);
359 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
360 CHECK(!work,MEM);
361 zgeev_ (&jobvl,&jobvr,
362 &n,ap,&n,
363 sp,
364 up,&n,
365 vp,&n,
366 work, &lwork,
367 rwork,
368 &res);
369 CHECK(res,res);
370 free(work);
371 free(rwork);
372 OK
373}
374
375
376
377//////////////////// general real eigensystem ////////////
378
379int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *
380 a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl,
381 integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
382 integer *lwork, integer *info);
383
384int eig_l_R(ODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
385 integer n = ar;
386 REQUIRES(ac==n && sn==n, BAD_SIZE);
387 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
388 char jobvl = up==NULL?'N':'V';
389 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
390 char jobvr = vp==NULL?'N':'V';
391 DEBUGMSG("eig_l_R");
392 integer lwork = -1;
393 integer res;
394 // ask for optimal lwork
395 double ans;
396 dgeev_ (&jobvl,&jobvr,
397 &n,ap,&n,
398 (double*)sp, (double*)sp+n,
399 up,&n,
400 vp,&n,
401 &ans, &lwork,
402 &res);
403 lwork = ceil(ans);
404 double * work = (double*)malloc(lwork*sizeof(double));
405 CHECK(!work,MEM);
406 dgeev_ (&jobvl,&jobvr,
407 &n,ap,&n,
408 (double*)sp, (double*)sp+n,
409 up,&n,
410 vp,&n,
411 work, &lwork,
412 &res);
413 CHECK(res,res);
414 free(work);
415 OK
416}
417
418
419//////////////////// symmetric real eigensystem ////////////
420
421int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a,
422 integer *lda, doublereal *w, doublereal *work, integer *lwork,
423 integer *info);
424
425int eig_l_S(int wantV,DVEC(s),ODMAT(v)) {
426 integer n = sn;
427 REQUIRES(vr==n && vc==n, BAD_SIZE);
428 char jobz = wantV?'V':'N';
429 DEBUGMSG("eig_l_S");
430 integer lwork = -1;
431 char uplo = 'U';
432 integer res;
433 // ask for optimal lwork
434 double ans;
435 dsyev_ (&jobz,&uplo,
436 &n,vp,&n,
437 sp,
438 &ans, &lwork,
439 &res);
440 lwork = ceil(ans);
441 double * work = (double*)malloc(lwork*sizeof(double));
442 CHECK(!work,MEM);
443 dsyev_ (&jobz,&uplo,
444 &n,vp,&n,
445 sp,
446 work, &lwork,
447 &res);
448 CHECK(res,res);
449 free(work);
450 OK
451}
452
453//////////////////// hermitian complex eigensystem ////////////
454
455int zheev_(char *jobz, char *uplo, integer *n, doublecomplex
456 *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork,
457 doublereal *rwork, integer *info);
458
459int eig_l_H(int wantV,DVEC(s),OCMAT(v)) {
460 integer n = sn;
461 REQUIRES(vr==n && vc==n, BAD_SIZE);
462 char jobz = wantV?'V':'N';
463 DEBUGMSG("eig_l_H");
464 double *rwork = (double*) malloc((3*n-2)*sizeof(double));
465 CHECK(!rwork,MEM);
466 integer lwork = -1;
467 char uplo = 'U';
468 integer res;
469 // ask for optimal lwork
470 doublecomplex ans;
471 zheev_ (&jobz,&uplo,
472 &n,vp,&n,
473 sp,
474 &ans, &lwork,
475 rwork,
476 &res);
477 lwork = ceil(ans.r);
478 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
479 CHECK(!work,MEM);
480 zheev_ (&jobz,&uplo,
481 &n,vp,&n,
482 sp,
483 work, &lwork,
484 rwork,
485 &res);
486 CHECK(res,res);
487 free(work);
488 free(rwork);
489 OK
490}
491
492//////////////////// general real linear system ////////////
493
494int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
495 *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info);
496
497int linearSolveR_l(ODMAT(a),ODMAT(b)) {
498 integer n = ar;
499 integer nhrs = bc;
500 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
501 DEBUGMSG("linearSolveR_l");
502 integer * ipiv = (integer*)malloc(n*sizeof(integer));
503 integer res;
504 dgesv_ (&n,&nhrs,
505 ap, &n,
506 ipiv,
507 bp, &n,
508 &res);
509 if(res>0) {
510 return SINGULAR;
511 }
512 CHECK(res,res);
513 free(ipiv);
514 OK
515}
516
517//////////////////// general complex linear system ////////////
518
519int zgesv_(integer *n, integer *nrhs, doublecomplex *a,
520 integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer *
521 info);
522
523int linearSolveC_l(OCMAT(a),OCMAT(b)) {
524 integer n = ar;
525 integer nhrs = bc;
526 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
527 DEBUGMSG("linearSolveC_l");
528 integer * ipiv = (integer*)malloc(n*sizeof(integer));
529 integer res;
530 zgesv_ (&n,&nhrs,
531 ap, &n,
532 ipiv,
533 bp, &n,
534 &res);
535 if(res>0) {
536 return SINGULAR;
537 }
538 CHECK(res,res);
539 free(ipiv);
540 OK
541}
542
543//////// symmetric positive definite real linear system using Cholesky ////////////
544
545int dpotrs_(char *uplo, integer *n, integer *nrhs,
546 doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *
547 info);
548
549int cholSolveR_l(KODMAT(a),ODMAT(b)) {
550 integer n = ar;
551 integer lda = aXc;
552 integer nhrs = bc;
553 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
554 DEBUGMSG("cholSolveR_l");
555 integer res;
556 dpotrs_ ("U",
557 &n,&nhrs,
558 (double*)ap, &lda,
559 bp, &n,
560 &res);
561 CHECK(res,res);
562 OK
563}
564
565//////// Hermitian positive definite real linear system using Cholesky ////////////
566
567int zpotrs_(char *uplo, integer *n, integer *nrhs,
568 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
569 integer *info);
570
571int cholSolveC_l(KOCMAT(a),OCMAT(b)) {
572 integer n = ar;
573 integer lda = aXc;
574 integer nhrs = bc;
575 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
576 DEBUGMSG("cholSolveC_l");
577 integer res;
578 zpotrs_ ("U",
579 &n,&nhrs,
580 (doublecomplex*)ap, &lda,
581 bp, &n,
582 &res);
583 CHECK(res,res);
584 OK
585}
586
587//////////////////// least squares real linear system ////////////
588
589int dgels_(char *trans, integer *m, integer *n, integer *
590 nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb,
591 doublereal *work, integer *lwork, integer *info);
592
593int linearSolveLSR_l(ODMAT(a),ODMAT(b)) {
594 integer m = ar;
595 integer n = ac;
596 integer nrhs = bc;
597 integer ldb = bXc;
598 REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);
599 DEBUGMSG("linearSolveLSR_l");
600 integer res;
601 integer lwork = -1;
602 double ans;
603 dgels_ ("N",&m,&n,&nrhs,
604 ap,&m,
605 bp,&ldb,
606 &ans,&lwork,
607 &res);
608 lwork = ceil(ans);
609 double * work = (double*)malloc(lwork*sizeof(double));
610 dgels_ ("N",&m,&n,&nrhs,
611 ap,&m,
612 bp,&ldb,
613 work,&lwork,
614 &res);
615 if(res>0) {
616 return SINGULAR;
617 }
618 CHECK(res,res);
619 free(work);
620 OK
621}
622
623//////////////////// least squares complex linear system ////////////
624
625int zgels_(char *trans, integer *m, integer *n, integer *
626 nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
627 doublecomplex *work, integer *lwork, integer *info);
628
629int linearSolveLSC_l(OCMAT(a),OCMAT(b)) {
630 integer m = ar;
631 integer n = ac;
632 integer nrhs = bc;
633 integer ldb = bXc;
634 REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);
635 DEBUGMSG("linearSolveLSC_l");
636 integer res;
637 integer lwork = -1;
638 doublecomplex ans;
639 zgels_ ("N",&m,&n,&nrhs,
640 ap,&m,
641 bp,&ldb,
642 &ans,&lwork,
643 &res);
644 lwork = ceil(ans.r);
645 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
646 zgels_ ("N",&m,&n,&nrhs,
647 ap,&m,
648 bp,&ldb,
649 work,&lwork,
650 &res);
651 if(res>0) {
652 return SINGULAR;
653 }
654 CHECK(res,res);
655 free(work);
656 OK
657}
658
659//////////////////// least squares real linear system using SVD ////////////
660
661int dgelss_(integer *m, integer *n, integer *nrhs,
662 doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *
663 s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork,
664 integer *info);
665
666int linearSolveSVDR_l(double rcond,ODMAT(a),ODMAT(b)) {
667 integer m = ar;
668 integer n = ac;
669 integer nrhs = bc;
670 integer ldb = bXc;
671 REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);
672 DEBUGMSG("linearSolveSVDR_l");
673 double*S = (double*)malloc(MIN(m,n)*sizeof(double));
674 integer res;
675 integer lwork = -1;
676 integer rank;
677 double ans;
678 dgelss_ (&m,&n,&nrhs,
679 ap,&m,
680 bp,&ldb,
681 S,
682 &rcond,&rank,
683 &ans,&lwork,
684 &res);
685 lwork = ceil(ans);
686 double * work = (double*)malloc(lwork*sizeof(double));
687 dgelss_ (&m,&n,&nrhs,
688 ap,&m,
689 bp,&ldb,
690 S,
691 &rcond,&rank,
692 work,&lwork,
693 &res);
694 if(res>0) {
695 return NOCONVER;
696 }
697 CHECK(res,res);
698 free(work);
699 free(S);
700 OK
701}
702
703//////////////////// least squares complex linear system using SVD ////////////
704
705int zgelss_(integer *m, integer *n, integer *nhrs,
706 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublereal *s,
707 doublereal *rcond, integer* rank,
708 doublecomplex *work, integer* lwork, doublereal* rwork,
709 integer *info);
710
711int linearSolveSVDC_l(double rcond, OCMAT(a),OCMAT(b)) {
712 integer m = ar;
713 integer n = ac;
714 integer nrhs = bc;
715 integer ldb = bXc;
716 REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);
717 DEBUGMSG("linearSolveSVDC_l");
718 double*S = (double*)malloc(MIN(m,n)*sizeof(double));
719 double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double));
720 integer res;
721 integer lwork = -1;
722 integer rank;
723 doublecomplex ans;
724 zgelss_ (&m,&n,&nrhs,
725 ap,&m,
726 bp,&ldb,
727 S,
728 &rcond,&rank,
729 &ans,&lwork,
730 RWORK,
731 &res);
732 lwork = ceil(ans.r);
733 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
734 zgelss_ (&m,&n,&nrhs,
735 ap,&m,
736 bp,&ldb,
737 S,
738 &rcond,&rank,
739 work,&lwork,
740 RWORK,
741 &res);
742 if(res>0) {
743 return NOCONVER;
744 }
745 CHECK(res,res);
746 free(work);
747 free(RWORK);
748 free(S);
749 OK
750}
751
752//////////////////// Cholesky factorization /////////////////////////
753
754int zpotrf_(char *uplo, integer *n, doublecomplex *a, integer *lda, integer *info);
755
756int chol_l_H(OCMAT(l)) {
757 integer n = lr;
758 REQUIRES(n>=1 && lc == n,BAD_SIZE);
759 DEBUGMSG("chol_l_H");
760 char uplo = 'U';
761 integer res;
762 zpotrf_ (&uplo,&n,lp,&n,&res);
763 CHECK(res>0,NODEFPOS);
764 CHECK(res,res);
765 doublecomplex zero = {0.,0.};
766 int r,c;
767 for (r=0; r<lr; r++) {
768 for(c=0; c<r; c++) {
769 AT(l,r,c) = zero;
770 }
771 }
772 OK
773}
774
775
776int dpotrf_(char *uplo, integer *n, doublereal *a, integer * lda, integer *info);
777
778int chol_l_S(ODMAT(l)) {
779 integer n = lr;
780 REQUIRES(n>=1 && lc == n,BAD_SIZE);
781 DEBUGMSG("chol_l_S");
782 char uplo = 'U';
783 integer res;
784 dpotrf_ (&uplo,&n,lp,&n,&res);
785 CHECK(res>0,NODEFPOS);
786 CHECK(res,res);
787 int r,c;
788 for (r=0; r<lr; r++) {
789 for(c=0; c<r; c++) {
790 AT(l,r,c) = 0.;
791 }
792 }
793 OK
794}
795
796//////////////////// QR factorization /////////////////////////
797
798int dgeqr2_(integer *m, integer *n, doublereal *a, integer *
799 lda, doublereal *tau, doublereal *work, integer *info);
800
801int qr_l_R(DVEC(tau), ODMAT(r)) {
802 integer m = rr;
803 integer n = rc;
804 integer mn = MIN(m,n);
805 REQUIRES(m>=1 && n >=1 && taun == mn, BAD_SIZE);
806 DEBUGMSG("qr_l_R");
807 double *WORK = (double*)malloc(n*sizeof(double));
808 CHECK(!WORK,MEM);
809 integer res;
810 dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);
811 CHECK(res,res);
812 free(WORK);
813 OK
814}
815
816int zgeqr2_(integer *m, integer *n, doublecomplex *a,
817 integer *lda, doublecomplex *tau, doublecomplex *work, integer *info);
818
819int qr_l_C(CVEC(tau), OCMAT(r)) {
820 integer m = rr;
821 integer n = rc;
822 integer mn = MIN(m,n);
823 REQUIRES(m>=1 && n >=1 && taun == mn, BAD_SIZE);
824 DEBUGMSG("qr_l_C");
825 doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex));
826 CHECK(!WORK,MEM);
827 integer res;
828 zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);
829 CHECK(res,res);
830 free(WORK);
831 OK
832}
833
834int dorgqr_(integer *m, integer *n, integer *k, doublereal *
835 a, integer *lda, doublereal *tau, doublereal *work, integer *lwork,
836 integer *info);
837
838int c_dorgqr(KDVEC(tau), ODMAT(r)) {
839 integer m = rr;
840 integer n = MIN(rc,rr);
841 integer k = taun;
842 DEBUGMSG("c_dorgqr");
843 integer lwork = 8*n; // FIXME
844 double *WORK = (double*)malloc(lwork*sizeof(double));
845 CHECK(!WORK,MEM);
846 integer res;
847 dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res);
848 CHECK(res,res);
849 free(WORK);
850 OK
851}
852
853int zungqr_(integer *m, integer *n, integer *k,
854 doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
855 work, integer *lwork, integer *info);
856
857int c_zungqr(KCVEC(tau), OCMAT(r)) {
858 integer m = rr;
859 integer n = MIN(rc,rr);
860 integer k = taun;
861 DEBUGMSG("z_ungqr");
862 integer lwork = 8*n; // FIXME
863 doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
864 CHECK(!WORK,MEM);
865 integer res;
866 zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res);
867 CHECK(res,res);
868 free(WORK);
869 OK
870}
871
872
873//////////////////// Hessenberg factorization /////////////////////////
874
875int dgehrd_(integer *n, integer *ilo, integer *ihi,
876 doublereal *a, integer *lda, doublereal *tau, doublereal *work,
877 integer *lwork, integer *info);
878
879int hess_l_R(DVEC(tau), ODMAT(r)) {
880 integer m = rr;
881 integer n = rc;
882 integer mn = MIN(m,n);
883 REQUIRES(m>=1 && n == m && taun == mn-1, BAD_SIZE);
884 DEBUGMSG("hess_l_R");
885 integer lwork = 5*n; // FIXME
886 double *WORK = (double*)malloc(lwork*sizeof(double));
887 CHECK(!WORK,MEM);
888 integer res;
889 integer one = 1;
890 dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);
891 CHECK(res,res);
892 free(WORK);
893 OK
894}
895
896
897int zgehrd_(integer *n, integer *ilo, integer *ihi,
898 doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
899 work, integer *lwork, integer *info);
900
901int hess_l_C(CVEC(tau), OCMAT(r)) {
902 integer m = rr;
903 integer n = rc;
904 integer mn = MIN(m,n);
905 REQUIRES(m>=1 && n == m && taun == mn-1, BAD_SIZE);
906 DEBUGMSG("hess_l_C");
907 integer lwork = 5*n; // FIXME
908 doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
909 CHECK(!WORK,MEM);
910 integer res;
911 integer one = 1;
912 zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);
913 CHECK(res,res);
914 free(WORK);
915 OK
916}
917
918//////////////////// Schur factorization /////////////////////////
919
920int dgees_(char *jobvs, char *sort, L_fp select, integer *n,
921 doublereal *a, integer *lda, integer *sdim, doublereal *wr,
922 doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work,
923 integer *lwork, logical *bwork, integer *info);
924
925int schur_l_R(ODMAT(u), ODMAT(s)) {
926 integer m = sr;
927 integer n = sc;
928 REQUIRES(m>=1 && n==m && ur==n && uc==n, BAD_SIZE);
929 DEBUGMSG("schur_l_R");
930 integer lwork = 6*n; // FIXME
931 double *WORK = (double*)malloc(lwork*sizeof(double));
932 double *WR = (double*)malloc(n*sizeof(double));
933 double *WI = (double*)malloc(n*sizeof(double));
934 // WR and WI not really required in this call
935 logical *BWORK = (logical*)malloc(n*sizeof(logical));
936 integer res;
937 integer sdim;
938 dgees_ ("V","N",NULL,&n,sp,&n,&sdim,WR,WI,up,&n,WORK,&lwork,BWORK,&res);
939 if(res>0) {
940 return NOCONVER;
941 }
942 CHECK(res,res);
943 free(WR);
944 free(WI);
945 free(BWORK);
946 free(WORK);
947 OK
948}
949
950
951int zgees_(char *jobvs, char *sort, L_fp select, integer *n,
952 doublecomplex *a, integer *lda, integer *sdim, doublecomplex *w,
953 doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork,
954 doublereal *rwork, logical *bwork, integer *info);
955
956int schur_l_C(OCMAT(u), OCMAT(s)) {
957 integer m = sr;
958 integer n = sc;
959 REQUIRES(m>=1 && n==m && ur==n && uc==n, BAD_SIZE);
960 DEBUGMSG("schur_l_C");
961 integer lwork = 6*n; // FIXME
962 doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
963 doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex));
964 // W not really required in this call
965 logical *BWORK = (logical*)malloc(n*sizeof(logical));
966 double *RWORK = (double*)malloc(n*sizeof(double));
967 integer res;
968 integer sdim;
969 zgees_ ("V","N",NULL,&n,sp,&n,&sdim,W,
970 up,&n,
971 WORK,&lwork,RWORK,BWORK,&res);
972 if(res>0) {
973 return NOCONVER;
974 }
975 CHECK(res,res);
976 free(W);
977 free(BWORK);
978 free(WORK);
979 OK
980}
981
982//////////////////// LU factorization /////////////////////////
983
984int dgetrf_(integer *m, integer *n, doublereal *a, integer *
985 lda, integer *ipiv, integer *info);
986
987int lu_l_R(DVEC(ipiv), ODMAT(r)) {
988 integer m = rr;
989 integer n = rc;
990 integer mn = MIN(m,n);
991 REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
992 DEBUGMSG("lu_l_R");
993 integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
994 integer res;
995 dgetrf_ (&m,&n,rp,&m,auxipiv,&res);
996 if(res>0) {
997 res = 0; // FIXME
998 }
999 CHECK(res,res);
1000 int k;
1001 for (k=0; k<mn; k++) {
1002 ipivp[k] = auxipiv[k];
1003 }
1004 free(auxipiv);
1005 OK
1006}
1007
1008
1009int zgetrf_(integer *m, integer *n, doublecomplex *a,
1010 integer *lda, integer *ipiv, integer *info);
1011
1012int lu_l_C(DVEC(ipiv), OCMAT(r)) {
1013 integer m = rr;
1014 integer n = rc;
1015 integer mn = MIN(m,n);
1016 REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
1017 DEBUGMSG("lu_l_C");
1018 integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
1019 integer res;
1020 zgetrf_ (&m,&n,rp,&m,auxipiv,&res);
1021 if(res>0) {
1022 res = 0; // FIXME
1023 }
1024 CHECK(res,res);
1025 int k;
1026 for (k=0; k<mn; k++) {
1027 ipivp[k] = auxipiv[k];
1028 }
1029 free(auxipiv);
1030 OK
1031}
1032
1033
1034//////////////////// LU substitution /////////////////////////
1035
1036int dgetrs_(char *trans, integer *n, integer *nrhs,
1037 doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *
1038 ldb, integer *info);
1039
1040int luS_l_R(KODMAT(a), KDVEC(ipiv), ODMAT(b)) {
1041 integer m = ar;
1042 integer n = ac;
1043 integer lda = aXc;
1044 integer mrhs = br;
1045 integer nrhs = bc;
1046
1047 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1048 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1049 int k;
1050 for (k=0; k<n; k++) {
1051 auxipiv[k] = (integer)ipivp[k];
1052 }
1053 integer res;
1054 dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&lda,auxipiv,bp,&mrhs,&res);
1055 CHECK(res,res);
1056 free(auxipiv);
1057 OK
1058}
1059
1060
1061int zgetrs_(char *trans, integer *n, integer *nrhs,
1062 doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b,
1063 integer *ldb, integer *info);
1064
1065int luS_l_C(KOCMAT(a), KDVEC(ipiv), OCMAT(b)) {
1066 integer m = ar;
1067 integer n = ac;
1068 integer lda = aXc;
1069 integer mrhs = br;
1070 integer nrhs = bc;
1071
1072 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1073 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1074 int k;
1075 for (k=0; k<n; k++) {
1076 auxipiv[k] = (integer)ipivp[k];
1077 }
1078 integer res;
1079 zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&lda,auxipiv,bp,&mrhs,&res);
1080 CHECK(res,res);
1081 free(auxipiv);
1082 OK
1083}
1084
1085
1086//////////////////// LDL factorization /////////////////////////
1087
1088int dsytrf_(char *uplo, integer *n, doublereal *a, integer *lda, integer *ipiv,
1089 doublereal *work, integer *lwork, integer *info);
1090
1091int ldl_R(DVEC(ipiv), ODMAT(r)) {
1092 integer n = rr;
1093 REQUIRES(n>=1 && rc==n && ipivn == n, BAD_SIZE);
1094 DEBUGMSG("ldl_R");
1095 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1096 integer res;
1097 integer lda = rXc;
1098 integer lwork = -1;
1099 doublereal ans;
1100 dsytrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res);
1101 lwork = ceil(ans);
1102 doublereal* work = (doublereal*)malloc(lwork*sizeof(doublereal));
1103 dsytrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res);
1104 CHECK(res,res);
1105 int k;
1106 for (k=0; k<n; k++) {
1107 ipivp[k] = auxipiv[k];
1108 }
1109 free(auxipiv);
1110 free(work);
1111 OK
1112}
1113
1114
1115int zhetrf_(char *uplo, integer *n, doublecomplex *a, integer *lda, integer *ipiv,
1116 doublecomplex *work, integer *lwork, integer *info);
1117
1118int ldl_C(DVEC(ipiv), OCMAT(r)) {
1119 integer n = rr;
1120 REQUIRES(n>=1 && rc==n && ipivn == n, BAD_SIZE);
1121 DEBUGMSG("ldl_R");
1122 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1123 integer res;
1124 integer lda = rXc;
1125 integer lwork = -1;
1126 doublecomplex ans;
1127 zhetrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res);
1128 lwork = ceil(ans.r);
1129 doublecomplex* work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
1130 zhetrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res);
1131 CHECK(res,res);
1132 int k;
1133 for (k=0; k<n; k++) {
1134 ipivp[k] = auxipiv[k];
1135 }
1136 free(auxipiv);
1137 free(work);
1138 OK
1139
1140}
1141
1142//////////////////// LDL solve /////////////////////////
1143
1144int dsytrs_(char *uplo, integer *n, integer *nrhs, doublereal *a, integer *lda,
1145 integer *ipiv, doublereal *b, integer *ldb, integer *info);
1146
1147int ldl_S_R(KODMAT(a), KDVEC(ipiv), ODMAT(b)) {
1148 integer m = ar;
1149 integer n = ac;
1150 integer lda = aXc;
1151 integer mrhs = br;
1152 integer nrhs = bc;
1153
1154 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1155 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1156 int k;
1157 for (k=0; k<n; k++) {
1158 auxipiv[k] = (integer)ipivp[k];
1159 }
1160 integer res;
1161 dsytrs_ ("L",&n,&nrhs,(/*no const (!?)*/ double*)ap,&lda,auxipiv,bp,&mrhs,&res);
1162 CHECK(res,res);
1163 free(auxipiv);
1164 OK
1165}
1166
1167
1168int zhetrs_(char *uplo, integer *n, integer *nrhs, doublecomplex *a, integer *lda,
1169 integer *ipiv, doublecomplex *b, integer *ldb, integer *info);
1170
1171int ldl_S_C(KOCMAT(a), KDVEC(ipiv), OCMAT(b)) {
1172 integer m = ar;
1173 integer n = ac;
1174 integer lda = aXc;
1175 integer mrhs = br;
1176 integer nrhs = bc;
1177
1178 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1179 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1180 int k;
1181 for (k=0; k<n; k++) {
1182 auxipiv[k] = (integer)ipivp[k];
1183 }
1184 integer res;
1185 zhetrs_ ("L",&n,&nrhs,(doublecomplex*)ap,&lda,auxipiv,bp,&mrhs,&res);
1186 CHECK(res,res);
1187 free(auxipiv);
1188 OK
1189}
1190
1191
1192//////////////////// Matrix Product /////////////////////////
1193
1194void dgemm_(char *, char *, integer *, integer *, integer *,
1195 double *, const double *, integer *, const double *,
1196 integer *, double *, double *, integer *);
1197
1198int multiplyR(int ta, int tb, KODMAT(a),KODMAT(b),ODMAT(r)) {
1199 DEBUGMSG("dgemm_");
1200 CHECKNANR(a,"NaN multR Input\n")
1201 CHECKNANR(b,"NaN multR Input\n")
1202 integer m = ta?ac:ar;
1203 integer n = tb?br:bc;
1204 integer k = ta?ar:ac;
1205 integer lda = aXc;
1206 integer ldb = bXc;
1207 integer ldc = rXc;
1208 double alpha = 1;
1209 double beta = 0;
1210 dgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);
1211 CHECKNANR(r,"NaN multR Output\n")
1212 OK
1213}
1214
1215void zgemm_(char *, char *, integer *, integer *, integer *,
1216 doublecomplex *, const doublecomplex *, integer *, const doublecomplex *,
1217 integer *, doublecomplex *, doublecomplex *, integer *);
1218
1219int multiplyC(int ta, int tb, KOCMAT(a),KOCMAT(b),OCMAT(r)) {
1220 DEBUGMSG("zgemm_");
1221 CHECKNANC(a,"NaN multC Input\n")
1222 CHECKNANC(b,"NaN multC Input\n")
1223 integer m = ta?ac:ar;
1224 integer n = tb?br:bc;
1225 integer k = ta?ar:ac;
1226 integer lda = aXc;
1227 integer ldb = bXc;
1228 integer ldc = rXc;
1229 doublecomplex alpha = {1,0};
1230 doublecomplex beta = {0,0};
1231 zgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,
1232 ap,&lda,
1233 bp,&ldb,&beta,
1234 rp,&ldc);
1235 CHECKNANC(r,"NaN multC Output\n")
1236 OK
1237}
1238
1239void sgemm_(char *, char *, integer *, integer *, integer *,
1240 float *, const float *, integer *, const float *,
1241 integer *, float *, float *, integer *);
1242
1243int multiplyF(int ta, int tb, KOFMAT(a),KOFMAT(b),OFMAT(r)) {
1244 DEBUGMSG("sgemm_");
1245 integer m = ta?ac:ar;
1246 integer n = tb?br:bc;
1247 integer k = ta?ar:ac;
1248 integer lda = aXc;
1249 integer ldb = bXc;
1250 integer ldc = rXc;
1251 float alpha = 1;
1252 float beta = 0;
1253 sgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);
1254 OK
1255}
1256
1257void cgemm_(char *, char *, integer *, integer *, integer *,
1258 complex *, const complex *, integer *, const complex *,
1259 integer *, complex *, complex *, integer *);
1260
1261int multiplyQ(int ta, int tb, KOQMAT(a),KOQMAT(b),OQMAT(r)) {
1262 DEBUGMSG("cgemm_");
1263 integer m = ta?ac:ar;
1264 integer n = tb?br:bc;
1265 integer k = ta?ar:ac;
1266 integer lda = aXc;
1267 integer ldb = bXc;
1268 integer ldc = rXc;
1269 complex alpha = {1,0};
1270 complex beta = {0,0};
1271 cgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,
1272 ap,&lda,
1273 bp,&ldb,&beta,
1274 rp,&ldc);
1275 OK
1276}
1277
1278
1279#define MULT_IMP_VER(OP) \
1280 { TRAV(r,i,j) { \
1281 int k; \
1282 AT(r,i,j) = 0; \
1283 for (k=0;k<ac;k++) { \
1284 OP \
1285 } \
1286 } \
1287 }
1288
1289#define MULT_IMP(M) { \
1290 if (m==1) { \
1291 MULT_IMP_VER( AT(r,i,j) += AT(a,i,k) * AT(b,k,j); ) \
1292 } else { \
1293 MULT_IMP_VER( AT(r,i,j) = M(AT(r,i,j) + M(AT(a,i,k) * AT(b,k,j), m) , m) ; ) \
1294 } OK }
1295
1296int multiplyI(int m, KOIMAT(a), KOIMAT(b), OIMAT(r)) MULT_IMP(mod)
1297int multiplyL(int64_t m, KOLMAT(a), KOLMAT(b), OLMAT(r)) MULT_IMP(mod_l)
1298
1299/////////////////////////////// inplace row ops ////////////////////////////////
1300
1301#define AXPY_IMP { \
1302 int j; \
1303 for(j=j1; j<=j2; j++) { \
1304 AT(r,i2,j) += a*AT(r,i1,j); \
1305 } OK }
1306
1307#define AXPY_MOD_IMP(M) { \
1308 int j; \
1309 for(j=j1; j<=j2; j++) { \
1310 AT(r,i2,j) = M(AT(r,i2,j) + M(a*AT(r,i1,j), m) , m); \
1311 } OK }
1312
1313
1314#define SCAL_IMP { \
1315 int i,j; \
1316 for(i=i1; i<=i2; i++) { \
1317 for(j=j1; j<=j2; j++) { \
1318 AT(r,i,j) = a*AT(r,i,j); \
1319 } \
1320 } OK }
1321
1322#define SCAL_MOD_IMP(M) { \
1323 int i,j; \
1324 for(i=i1; i<=i2; i++) { \
1325 for(j=j1; j<=j2; j++) { \
1326 AT(r,i,j) = M(a*AT(r,i,j) , m); \
1327 } \
1328 } OK }
1329
1330
1331#define SWAP_IMP(T) { \
1332 T aux; \
1333 int k; \
1334 if (i1 != i2) { \
1335 for (k=j1; k<=j2; k++) { \
1336 aux = AT(r,i1,k); \
1337 AT(r,i1,k) = AT(r,i2,k); \
1338 AT(r,i2,k) = aux; \
1339 } \
1340 } OK }
1341
1342
1343#define ROWOP_IMP(T) { \
1344 T a = *pa; \
1345 switch(code) { \
1346 case 0: AXPY_IMP \
1347 case 1: SCAL_IMP \
1348 case 2: SWAP_IMP(T) \
1349 default: ERROR(BAD_CODE); \
1350 } \
1351}
1352
1353#define ROWOP_MOD_IMP(T,M) { \
1354 T a = *pa; \
1355 switch(code) { \
1356 case 0: AXPY_MOD_IMP(M) \
1357 case 1: SCAL_MOD_IMP(M) \
1358 case 2: SWAP_IMP(T) \
1359 default: ERROR(BAD_CODE); \
1360 } \
1361}
1362
1363
1364#define ROWOP(T) int rowop_##T(int code, T* pa, int i1, int i2, int j1, int j2, MATG(T,r)) ROWOP_IMP(T)
1365
1366#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)
1367
1368ROWOP(double)
1369ROWOP(float)
1370ROWOP(TCD)
1371ROWOP(TCF)
1372ROWOP(int32_t)
1373ROWOP(int64_t)
1374ROWOP_MOD(int32_t,mod)
1375ROWOP_MOD(int64_t,mod_l)
1376
1377/////////////////////////////// inplace GEMM ////////////////////////////////
1378
1379#define GEMM(T) int gemm_##T(VECG(T,c),MATG(T,a),MATG(T,b),MATG(T,r)) { \
1380 T a = cp[0], b = cp[1]; \
1381 T t; \
1382 int k; \
1383 { TRAV(r,i,j) { \
1384 t = 0; \
1385 for(k=0; k<ac; k++) { \
1386 t += AT(a,i,k) * AT(b,k,j); \
1387 } \
1388 AT(r,i,j) = b*AT(r,i,j) + a*t; \
1389 } \
1390 } OK }
1391
1392
1393GEMM(double)
1394GEMM(float)
1395GEMM(TCD)
1396GEMM(TCF)
1397GEMM(int32_t)
1398GEMM(int64_t)
1399
1400#define GEMM_MOD(T,M) int gemm_mod_##T(T m, VECG(T,c),MATG(T,a),MATG(T,b),MATG(T,r)) { \
1401 T a = cp[0], b = cp[1]; \
1402 int k; \
1403 T t; \
1404 { TRAV(r,i,j) { \
1405 t = 0; \
1406 for(k=0; k<ac; k++) { \
1407 t = M(t+M(AT(a,i,k) * AT(b,k,j))); \
1408 } \
1409 AT(r,i,j) = M(M(b*AT(r,i,j)) + M(a*t)); \
1410 } \
1411 } OK }
1412
1413
1414#define MOD32(X) mod(X,m)
1415#define MOD64(X) mod_l(X,m)
1416
1417GEMM_MOD(int32_t,MOD32)
1418GEMM_MOD(int64_t,MOD64)
1419
1420////////////////// sparse matrix-product ///////////////////////////////////////
1421
1422
1423int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
1424 int r, c;
1425 for (r = 0; r < rowsn - 1; r++) {
1426 rp[r] = 0;
1427 for (c = rowsp[r]; c < rowsp[r+1]; c++) {
1428 rp[r] += valsp[c-1] * xp[colsp[c-1]-1];
1429 }
1430 }
1431 OK
1432}
1433
1434int smTXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
1435 int r,c;
1436 for (c = 0; c < rn; c++) {
1437 rp[c] = 0;
1438 }
1439 for (r = 0; r < rowsn - 1; r++) {
1440 for (c = rowsp[r]; c < rowsp[r+1]; c++) {
1441 rp[colsp[c-1]-1] += valsp[c-1] * xp[r];
1442 }
1443 }
1444 OK
1445}
1446
1447
1448//////////////////////// extract /////////////////////////////////
1449
1450#define EXTRACT_IMP { \
1451 int i,j,si,sj,ni,nj; \
1452 ni = modei ? in : ip[1]-ip[0]+1; \
1453 nj = modej ? jn : jp[1]-jp[0]+1; \
1454 \
1455 for (i=0; i<ni; i++) { \
1456 si = modei ? ip[i] : i+ip[0]; \
1457 \
1458 for (j=0; j<nj; j++) { \
1459 sj = modej ? jp[j] : j+jp[0]; \
1460 \
1461 AT(r,i,j) = AT(m,si,sj); \
1462 } \
1463 } OK }
1464
1465#define EXTRACT(T) int extract##T(int modei, int modej, KIVEC(i), KIVEC(j), KO##T##MAT(m), O##T##MAT(r)) EXTRACT_IMP
1466
1467EXTRACT(D)
1468EXTRACT(F)
1469EXTRACT(C)
1470EXTRACT(Q)
1471EXTRACT(I)
1472EXTRACT(L)
1473
1474//////////////////////// setRect /////////////////////////////////
1475
1476#define SETRECT(T) \
1477int setRect##T(int i, int j, KO##T##MAT(m), O##T##MAT(r)) { \
1478 { TRAV(m,a,b) { \
1479 int x = a+i, y = b+j; \
1480 if(x>=0 && x<rr && y>=0 && y<rc) { \
1481 AT(r,x,y) = AT(m,a,b); \
1482 } \
1483 } \
1484 } OK }
1485
1486SETRECT(D)
1487SETRECT(F)
1488SETRECT(C)
1489SETRECT(Q)
1490SETRECT(I)
1491SETRECT(L)
1492
1493//////////////////////// remap /////////////////////////////////
1494
1495#define REMAP_IMP \
1496 REQUIRES(ir==jr && ic==jc && ir==rr && ic==rc ,BAD_SIZE); \
1497 { TRAV(r,a,b) { AT(r,a,b) = AT(m,AT(i,a,b),AT(j,a,b)); } \
1498 } \
1499 OK
1500
1501int remapD(KOIMAT(i), KOIMAT(j), KODMAT(m), ODMAT(r)) {
1502 REMAP_IMP
1503}
1504
1505int remapF(KOIMAT(i), KOIMAT(j), KOFMAT(m), OFMAT(r)) {
1506 REMAP_IMP
1507}
1508
1509int remapI(KOIMAT(i), KOIMAT(j), KOIMAT(m), OIMAT(r)) {
1510 REMAP_IMP
1511}
1512
1513int remapL(KOIMAT(i), KOIMAT(j), KOLMAT(m), OLMAT(r)) {
1514 REMAP_IMP
1515}
1516
1517int remapC(KOIMAT(i), KOIMAT(j), KOCMAT(m), OCMAT(r)) {
1518 REMAP_IMP
1519}
1520
1521int remapQ(KOIMAT(i), KOIMAT(j), KOQMAT(m), OQMAT(r)) {
1522 REMAP_IMP
1523}
1524
1525////////////////////////////////////////////////////////////////////////////////
1526
1527int saveMatrix(char * file, char * format, KODMAT(a)){
1528 FILE * fp;
1529 fp = fopen (file, "w");
1530 int r, c;
1531 for (r=0;r<ar; r++) {
1532 for (c=0; c<ac; c++) {
1533 fprintf(fp,format,AT(a,r,c));
1534 if (c<ac-1) {
1535 fprintf(fp," ");
1536 } else {
1537 fprintf(fp,"\n");
1538 }
1539 }
1540 }
1541 fclose(fp);
1542 OK
1543}
1544
diff --git a/packages/base/src/Internal/C/lapack-aux.h b/packages/base/src/Internal/C/lapack-aux.h
new file mode 100644
index 0000000..7a6fcbf
--- /dev/null
+++ b/packages/base/src/Internal/C/lapack-aux.h
@@ -0,0 +1,111 @@
1/*
2 * We have copied the definitions in f2c.h required
3 * to compile clapack.h, modified to support both
4 * 32 and 64 bit
5
6 http://opengrok.creo.hu/dragonfly/xref/src/contrib/gcc-3.4/libf2c/readme.netlib
7 http://www.ibm.com/developerworks/library/l-port64.html
8 */
9
10#ifdef _LP64
11typedef int integer;
12typedef unsigned int uinteger;
13typedef int logical;
14typedef long longint; /* system-dependent */
15typedef unsigned long ulongint; /* system-dependent */
16#else
17typedef long int integer;
18typedef unsigned long int uinteger;
19typedef long int logical;
20typedef long long longint; /* system-dependent */
21typedef unsigned long long ulongint; /* system-dependent */
22#endif
23
24typedef char *address;
25typedef short int shortint;
26typedef float real;
27typedef double doublereal;
28typedef struct { real r, i; } complex;
29typedef struct { doublereal r, i; } doublecomplex;
30typedef short int shortlogical;
31typedef char logical1;
32typedef char integer1;
33
34typedef logical (*L_fp)();
35typedef short ftnlen;
36
37/********************************************************/
38
39#define IVEC(A) int A##n, int*A##p
40#define LVEC(A) int A##n, int64_t*A##p
41#define FVEC(A) int A##n, float*A##p
42#define DVEC(A) int A##n, double*A##p
43#define QVEC(A) int A##n, complex*A##p
44#define CVEC(A) int A##n, doublecomplex*A##p
45#define PVEC(A) int A##n, void* A##p, int A##s
46
47#define IMAT(A) int A##r, int A##c, int* A##p
48#define LMAT(A) int A##r, int A##c, int64_t* A##p
49#define FMAT(A) int A##r, int A##c, float* A##p
50#define DMAT(A) int A##r, int A##c, double* A##p
51#define QMAT(A) int A##r, int A##c, complex* A##p
52#define CMAT(A) int A##r, int A##c, doublecomplex* A##p
53#define PMAT(A) int A##r, int A##c, void* A##p, int A##s
54
55#define KIVEC(A) int A##n, const int*A##p
56#define KLVEC(A) int A##n, const int64_t*A##p
57#define KFVEC(A) int A##n, const float*A##p
58#define KDVEC(A) int A##n, const double*A##p
59#define KQVEC(A) int A##n, const complex*A##p
60#define KCVEC(A) int A##n, const doublecomplex*A##p
61#define KPVEC(A) int A##n, const void* A##p, int A##s
62
63#define KIMAT(A) int A##r, int A##c, const int* A##p
64#define KLMAT(A) int A##r, int A##c, const int64_t* A##p
65#define KFMAT(A) int A##r, int A##c, const float* A##p
66#define KDMAT(A) int A##r, int A##c, const double* A##p
67#define KQMAT(A) int A##r, int A##c, const complex* A##p
68#define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p
69#define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s
70
71#define VECG(T,A) int A##n, T* A##p
72#define MATG(T,A) int A##r, int A##c, int A##Xr, int A##Xc, T* A##p
73
74#define OIMAT(A) MATG(int,A)
75#define OLMAT(A) MATG(int64_t,A)
76#define OFMAT(A) MATG(float,A)
77#define ODMAT(A) MATG(double,A)
78#define OQMAT(A) MATG(complex,A)
79#define OCMAT(A) MATG(doublecomplex,A)
80
81#define KOIMAT(A) MATG(const int,A)
82#define KOLMAT(A) MATG(const int64_t,A)
83#define KOFMAT(A) MATG(const float,A)
84#define KODMAT(A) MATG(const double,A)
85#define KOQMAT(A) MATG(const complex,A)
86#define KOCMAT(A) MATG(const doublecomplex,A)
87
88#define AT(m,i,j) (m##p[(i)*m##Xr + (j)*m##Xc])
89#define TRAV(m,i,j) int i,j; for (i=0;i<m##r;i++) for (j=0;j<m##c;j++)
90
91/********************************************************/
92
93static inline
94int mod (int a, int b) {
95 int m = a % b;
96 if (b>0) {
97 return m >=0 ? m : m+b;
98 } else {
99 return m <=0 ? m : m+b;
100 }
101}
102
103static inline
104int64_t mod_l (int64_t a, int64_t b) {
105 int64_t m = a % b;
106 if (b>0) {
107 return m >=0 ? m : m+b;
108 } else {
109 return m <=0 ? m : m+b;
110 }
111}
diff --git a/packages/base/src/Internal/C/vector-aux.c b/packages/base/src/Internal/C/vector-aux.c
new file mode 100644
index 0000000..9dbf536
--- /dev/null
+++ b/packages/base/src/Internal/C/vector-aux.c
@@ -0,0 +1,1486 @@
1#include <complex.h>
2#include <inttypes.h>
3
4typedef double complex TCD;
5typedef float complex TCF;
6
7#undef complex
8
9#include "lapack-aux.h"
10
11#define V(x) x##n,x##p
12
13#include <string.h>
14#include <math.h>
15#include <stdio.h>
16#include <stdlib.h>
17#include <stdint.h>
18
19#define MACRO(B) do {B} while (0)
20#define ERROR(CODE) MACRO(return CODE;)
21#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
22#define OK return 0;
23
24#define MIN(A,B) ((A)<(B)?(A):(B))
25#define MAX(A,B) ((A)>(B)?(A):(B))
26
27#ifdef DBG
28#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
29#else
30#define DEBUGMSG(M)
31#endif
32
33#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
34
35#define BAD_SIZE 2000
36#define BAD_CODE 2001
37#define MEM 2002
38#define BAD_FILE 2003
39
40
41int sumF(KFVEC(x),FVEC(r)) {
42 DEBUGMSG("sumF");
43 REQUIRES(rn==1,BAD_SIZE);
44 int i;
45 float res = 0;
46 for (i = 0; i < xn; i++) res += xp[i];
47 rp[0] = res;
48 OK
49}
50
51int sumR(KDVEC(x),DVEC(r)) {
52 DEBUGMSG("sumR");
53 REQUIRES(rn==1,BAD_SIZE);
54 int i;
55 double res = 0;
56 for (i = 0; i < xn; i++) res += xp[i];
57 rp[0] = res;
58 OK
59}
60
61int sumI(int m, KIVEC(x),IVEC(r)) {
62 REQUIRES(rn==1,BAD_SIZE);
63 int i;
64 int res = 0;
65 if (m==1) {
66 for (i = 0; i < xn; i++) res += xp[i];
67 } else {
68 for (i = 0; i < xn; i++) res = (res + xp[i]) % m;
69 }
70 rp[0] = res;
71 OK
72}
73
74int sumL(int64_t m, KLVEC(x),LVEC(r)) {
75 REQUIRES(rn==1,BAD_SIZE);
76 int i;
77 int res = 0;
78 if (m==1) {
79 for (i = 0; i < xn; i++) res += xp[i];
80 } else {
81 for (i = 0; i < xn; i++) res = (res + xp[i]) % m;
82 }
83 rp[0] = res;
84 OK
85}
86
87int sumQ(KQVEC(x),QVEC(r)) {
88 DEBUGMSG("sumQ");
89 REQUIRES(rn==1,BAD_SIZE);
90 int i;
91 complex res;
92 res.r = 0;
93 res.i = 0;
94 for (i = 0; i < xn; i++) {
95 res.r += xp[i].r;
96 res.i += xp[i].i;
97 }
98 rp[0] = res;
99 OK
100}
101
102int sumC(KCVEC(x),CVEC(r)) {
103 DEBUGMSG("sumC");
104 REQUIRES(rn==1,BAD_SIZE);
105 int i;
106 doublecomplex res;
107 res.r = 0;
108 res.i = 0;
109 for (i = 0; i < xn; i++) {
110 res.r += xp[i].r;
111 res.i += xp[i].i;
112 }
113 rp[0] = res;
114 OK
115}
116
117
118int prodF(KFVEC(x),FVEC(r)) {
119 DEBUGMSG("prodF");
120 REQUIRES(rn==1,BAD_SIZE);
121 int i;
122 float res = 1;
123 for (i = 0; i < xn; i++) res *= xp[i];
124 rp[0] = res;
125 OK
126}
127
128int prodR(KDVEC(x),DVEC(r)) {
129 DEBUGMSG("prodR");
130 REQUIRES(rn==1,BAD_SIZE);
131 int i;
132 double res = 1;
133 for (i = 0; i < xn; i++) res *= xp[i];
134 rp[0] = res;
135 OK
136}
137
138int prodI(int m, KIVEC(x),IVEC(r)) {
139 REQUIRES(rn==1,BAD_SIZE);
140 int i;
141 int res = 1;
142 if (m==1) {
143 for (i = 0; i < xn; i++) res *= xp[i];
144 } else {
145 for (i = 0; i < xn; i++) res = (res * xp[i]) % m;
146 }
147 rp[0] = res;
148 OK
149}
150
151int prodL(int64_t m, KLVEC(x),LVEC(r)) {
152 REQUIRES(rn==1,BAD_SIZE);
153 int i;
154 int res = 1;
155 if (m==1) {
156 for (i = 0; i < xn; i++) res *= xp[i];
157 } else {
158 for (i = 0; i < xn; i++) res = (res * xp[i]) % m;
159 }
160 rp[0] = res;
161 OK
162}
163
164int prodQ(KQVEC(x),QVEC(r)) {
165 DEBUGMSG("prodQ");
166 REQUIRES(rn==1,BAD_SIZE);
167 int i;
168 complex res;
169 float temp;
170 res.r = 1;
171 res.i = 0;
172 for (i = 0; i < xn; i++) {
173 temp = res.r * xp[i].r - res.i * xp[i].i;
174 res.i = res.r * xp[i].i + res.i * xp[i].r;
175 res.r = temp;
176 }
177 rp[0] = res;
178 OK
179}
180
181int prodC(KCVEC(x),CVEC(r)) {
182 DEBUGMSG("prodC");
183 REQUIRES(rn==1,BAD_SIZE);
184 int i;
185 doublecomplex res;
186 double temp;
187 res.r = 1;
188 res.i = 0;
189 for (i = 0; i < xn; i++) {
190 temp = res.r * xp[i].r - res.i * xp[i].i;
191 res.i = res.r * xp[i].i + res.i * xp[i].r;
192 res.r = temp;
193 }
194 rp[0] = res;
195 OK
196}
197
198
199double dnrm2_(integer*, const double*, integer*);
200double dasum_(integer*, const double*, integer*);
201
202double vector_max(KDVEC(x)) {
203 double r = xp[0];
204 int k;
205 for (k = 1; k<xn; k++) {
206 if(xp[k]>r) {
207 r = xp[k];
208 }
209 }
210 return r;
211}
212
213double vector_min(KDVEC(x)) {
214 double r = xp[0];
215 int k;
216 for (k = 1; k<xn; k++) {
217 if(xp[k]<r) {
218 r = xp[k];
219 }
220 }
221 return r;
222}
223
224int vector_max_index(KDVEC(x)) {
225 int k, r = 0;
226 for (k = 1; k<xn; k++) {
227 if(xp[k]>xp[r]) {
228 r = k;
229 }
230 }
231 return r;
232}
233
234int vector_min_index(KDVEC(x)) {
235 int k, r = 0;
236 for (k = 1; k<xn; k++) {
237 if(xp[k]<xp[r]) {
238 r = k;
239 }
240 }
241 return r;
242}
243
244int toScalarR(int code, KDVEC(x), DVEC(r)) {
245 REQUIRES(rn==1,BAD_SIZE);
246 DEBUGMSG("toScalarR");
247 double res;
248 integer one = 1;
249 integer n = xn;
250 switch(code) {
251 case 0: { res = dnrm2_(&n,xp,&one); break; }
252 case 1: { res = dasum_(&n,xp,&one); break; }
253 case 2: { res = vector_max_index(V(x)); break; }
254 case 3: { res = vector_max(V(x)); break; }
255 case 4: { res = vector_min_index(V(x)); break; }
256 case 5: { res = vector_min(V(x)); break; }
257 default: ERROR(BAD_CODE);
258 }
259 rp[0] = res;
260 OK
261}
262
263
264float snrm2_(integer*, const float*, integer*);
265float sasum_(integer*, const float*, integer*);
266
267float vector_max_f(KFVEC(x)) {
268 float r = xp[0];
269 int k;
270 for (k = 1; k<xn; k++) {
271 if(xp[k]>r) {
272 r = xp[k];
273 }
274 }
275 return r;
276}
277
278float vector_min_f(KFVEC(x)) {
279 float r = xp[0];
280 int k;
281 for (k = 1; k<xn; k++) {
282 if(xp[k]<r) {
283 r = xp[k];
284 }
285 }
286 return r;
287}
288
289int vector_max_index_f(KFVEC(x)) {
290 int k, r = 0;
291 for (k = 1; k<xn; k++) {
292 if(xp[k]>xp[r]) {
293 r = k;
294 }
295 }
296 return r;
297}
298
299int vector_min_index_f(KFVEC(x)) {
300 int k, r = 0;
301 for (k = 1; k<xn; k++) {
302 if(xp[k]<xp[r]) {
303 r = k;
304 }
305 }
306 return r;
307}
308
309
310int toScalarF(int code, KFVEC(x), FVEC(r)) {
311 REQUIRES(rn==1,BAD_SIZE);
312 DEBUGMSG("toScalarF");
313 float res;
314 integer one = 1;
315 integer n = xn;
316 switch(code) {
317 case 0: { res = snrm2_(&n,xp,&one); break; }
318 case 1: { res = sasum_(&n,xp,&one); break; }
319 case 2: { res = vector_max_index_f(V(x)); break; }
320 case 3: { res = vector_max_f(V(x)); break; }
321 case 4: { res = vector_min_index_f(V(x)); break; }
322 case 5: { res = vector_min_f(V(x)); break; }
323 default: ERROR(BAD_CODE);
324 }
325 rp[0] = res;
326 OK
327}
328
329int vector_max_i(KIVEC(x)) {
330 int r = xp[0];
331 int k;
332 for (k = 1; k<xn; k++) {
333 if(xp[k]>r) {
334 r = xp[k];
335 }
336 }
337 return r;
338}
339
340int vector_min_i(KIVEC(x)) {
341 int r = xp[0];
342 int k;
343 for (k = 1; k<xn; k++) {
344 if(xp[k]<r) {
345 r = xp[k];
346 }
347 }
348 return r;
349}
350
351int vector_max_index_i(KIVEC(x)) {
352 int k, r = 0;
353 for (k = 1; k<xn; k++) {
354 if(xp[k]>xp[r]) {
355 r = k;
356 }
357 }
358 return r;
359}
360
361int vector_min_index_i(KIVEC(x)) {
362 int k, r = 0;
363 for (k = 1; k<xn; k++) {
364 if(xp[k]<xp[r]) {
365 r = k;
366 }
367 }
368 return r;
369}
370
371
372int toScalarI(int code, KIVEC(x), IVEC(r)) {
373 REQUIRES(rn==1,BAD_SIZE);
374 int res;
375 switch(code) {
376 case 2: { res = vector_max_index_i(V(x)); break; }
377 case 3: { res = vector_max_i(V(x)); break; }
378 case 4: { res = vector_min_index_i(V(x)); break; }
379 case 5: { res = vector_min_i(V(x)); break; }
380 default: ERROR(BAD_CODE);
381 }
382 rp[0] = res;
383 OK
384}
385
386
387int64_t vector_max_l(KLVEC(x)) {
388 int64_t r = xp[0];
389 int k;
390 for (k = 1; k<xn; k++) {
391 if(xp[k]>r) {
392 r = xp[k];
393 }
394 }
395 return r;
396}
397
398int64_t vector_min_l(KLVEC(x)) {
399 int64_t r = xp[0];
400 int k;
401 for (k = 1; k<xn; k++) {
402 if(xp[k]<r) {
403 r = xp[k];
404 }
405 }
406 return r;
407}
408
409int vector_max_index_l(KLVEC(x)) {
410 int k, r = 0;
411 for (k = 1; k<xn; k++) {
412 if(xp[k]>xp[r]) {
413 r = k;
414 }
415 }
416 return r;
417}
418
419int vector_min_index_l(KLVEC(x)) {
420 int k, r = 0;
421 for (k = 1; k<xn; k++) {
422 if(xp[k]<xp[r]) {
423 r = k;
424 }
425 }
426 return r;
427}
428
429
430int toScalarL(int code, KLVEC(x), LVEC(r)) {
431 REQUIRES(rn==1,BAD_SIZE);
432 int64_t res;
433 switch(code) {
434 case 2: { res = vector_max_index_l(V(x)); break; }
435 case 3: { res = vector_max_l(V(x)); break; }
436 case 4: { res = vector_min_index_l(V(x)); break; }
437 case 5: { res = vector_min_l(V(x)); break; }
438 default: ERROR(BAD_CODE);
439 }
440 rp[0] = res;
441 OK
442}
443
444
445double dznrm2_(integer*, const doublecomplex*, integer*);
446double dzasum_(integer*, const doublecomplex*, integer*);
447
448int toScalarC(int code, KCVEC(x), DVEC(r)) {
449 REQUIRES(rn==1,BAD_SIZE);
450 DEBUGMSG("toScalarC");
451 double res;
452 integer one = 1;
453 integer n = xn;
454 switch(code) {
455 case 0: { res = dznrm2_(&n,xp,&one); break; }
456 case 1: { res = dzasum_(&n,xp,&one); break; }
457 default: ERROR(BAD_CODE);
458 }
459 rp[0] = res;
460 OK
461}
462
463
464double scnrm2_(integer*, const complex*, integer*);
465double scasum_(integer*, const complex*, integer*);
466
467int toScalarQ(int code, KQVEC(x), FVEC(r)) {
468 REQUIRES(rn==1,BAD_SIZE);
469 DEBUGMSG("toScalarQ");
470 float res;
471 integer one = 1;
472 integer n = xn;
473 switch(code) {
474 case 0: { res = scnrm2_(&n,xp,&one); break; }
475 case 1: { res = scasum_(&n,xp,&one); break; }
476 default: ERROR(BAD_CODE);
477 }
478 rp[0] = res;
479 OK
480}
481
482
483inline double sign(double x) {
484 if(x>0) {
485 return +1.0;
486 } else if (x<0) {
487 return -1.0;
488 } else {
489 return 0.0;
490 }
491}
492
493inline float float_sign(float x) {
494 if(x>0) {
495 return +1.0;
496 } else if (x<0) {
497 return -1.0;
498 } else {
499 return 0.0;
500 }
501}
502
503
504#define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
505#define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
506int mapR(int code, KDVEC(x), DVEC(r)) {
507 int k;
508 REQUIRES(xn == rn,BAD_SIZE);
509 DEBUGMSG("mapR");
510 switch (code) {
511 OP(0,sin)
512 OP(1,cos)
513 OP(2,tan)
514 OP(3,fabs)
515 OP(4,asin)
516 OP(5,acos)
517 OP(6,atan)
518 OP(7,sinh)
519 OP(8,cosh)
520 OP(9,tanh)
521 OP(10,asinh)
522 OP(11,acosh)
523 OP(12,atanh)
524 OP(13,exp)
525 OP(14,log)
526 OP(15,sign)
527 OP(16,sqrt)
528 default: ERROR(BAD_CODE);
529 }
530}
531
532int mapF(int code, KFVEC(x), FVEC(r)) {
533 int k;
534 REQUIRES(xn == rn,BAD_SIZE);
535 DEBUGMSG("mapF");
536 switch (code) {
537 OP(0,sin)
538 OP(1,cos)
539 OP(2,tan)
540 OP(3,fabs)
541 OP(4,asin)
542 OP(5,acos)
543 OP(6,atan)
544 OP(7,sinh)
545 OP(8,cosh)
546 OP(9,tanh)
547 OP(10,asinh)
548 OP(11,acosh)
549 OP(12,atanh)
550 OP(13,exp)
551 OP(14,log)
552 OP(15,sign)
553 OP(16,sqrt)
554 default: ERROR(BAD_CODE);
555 }
556}
557
558
559int mapI(int code, KIVEC(x), IVEC(r)) {
560 int k;
561 REQUIRES(xn == rn,BAD_SIZE);
562 switch (code) {
563 OP(3,abs)
564 OP(15,sign)
565 default: ERROR(BAD_CODE);
566 }
567}
568
569
570int mapL(int code, KLVEC(x), LVEC(r)) {
571 int k;
572 REQUIRES(xn == rn,BAD_SIZE);
573 switch (code) {
574 OP(3,abs)
575 OP(15,sign)
576 default: ERROR(BAD_CODE);
577 }
578}
579
580
581
582inline double abs_complex(doublecomplex z) {
583 return sqrt(z.r*z.r + z.i*z.i);
584}
585
586inline doublecomplex complex_abs_complex(doublecomplex z) {
587 doublecomplex r;
588 r.r = abs_complex(z);
589 r.i = 0;
590 return r;
591}
592
593inline doublecomplex complex_signum_complex(doublecomplex z) {
594 doublecomplex r;
595 double mag;
596 if (z.r == 0 && z.i == 0) {
597 r.r = 0;
598 r.i = 0;
599 } else {
600 mag = abs_complex(z);
601 r.r = z.r/mag;
602 r.i = z.i/mag;
603 }
604 return r;
605}
606
607#define OPb(C,F) case C: { for(k=0;k<xn;k++) r2p[k] = F(x2p[k]); OK }
608int mapC(int code, KCVEC(x), CVEC(r)) {
609 TCD* x2p = (TCD*)xp;
610 TCD* r2p = (TCD*)rp;
611 int k;
612 REQUIRES(xn == rn,BAD_SIZE);
613 DEBUGMSG("mapC");
614 switch (code) {
615 OPb(0,csin)
616 OPb(1,ccos)
617 OPb(2,ctan)
618 OP(3,complex_abs_complex)
619 OPb(4,casin)
620 OPb(5,cacos)
621 OPb(6,catan)
622 OPb(7,csinh)
623 OPb(8,ccosh)
624 OPb(9,ctanh)
625 OPb(10,casinh)
626 OPb(11,cacosh)
627 OPb(12,catanh)
628 OPb(13,cexp)
629 OPb(14,clog)
630 OP(15,complex_signum_complex)
631 OPb(16,csqrt)
632 default: ERROR(BAD_CODE);
633 }
634}
635
636
637
638inline complex complex_f_math_fun(doublecomplex (*cf)(doublecomplex), complex a)
639{
640 doublecomplex c;
641 doublecomplex r;
642
643 complex float_r;
644
645 c.r = a.r;
646 c.i = a.i;
647
648 r = (*cf)(c);
649
650 float_r.r = r.r;
651 float_r.i = r.i;
652
653 return float_r;
654}
655
656
657#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_f_math_fun(&F,xp[k]); OK }
658int mapQ(int code, KQVEC(x), QVEC(r)) {
659 TCF* x2p = (TCF*)xp;
660 TCF* r2p = (TCF*)rp;
661 int k;
662 REQUIRES(xn == rn,BAD_SIZE);
663 DEBUGMSG("mapQ");
664 switch (code) {
665 OPb(0,csinf)
666 OPb(1,ccosf)
667 OPb(2,ctanf)
668 OPC(3,complex_abs_complex)
669 OPb(4,casinf)
670 OPb(5,cacosf)
671 OPb(6,catanf)
672 OPb(7,csinhf)
673 OPb(8,ccoshf)
674 OPb(9,ctanhf)
675 OPb(10,casinhf)
676 OPb(11,cacoshf)
677 OPb(12,catanhf)
678 OPb(13,cexpf)
679 OPb(14,clogf)
680 OPC(15,complex_signum_complex)
681 OPb(16,csqrtf)
682 default: ERROR(BAD_CODE);
683 }
684}
685
686
687int mapValR(int code, double* pval, KDVEC(x), DVEC(r)) {
688 int k;
689 double val = *pval;
690 REQUIRES(xn == rn,BAD_SIZE);
691 DEBUGMSG("mapValR");
692 switch (code) {
693 OPV(0,val*xp[k])
694 OPV(1,val/xp[k])
695 OPV(2,val+xp[k])
696 OPV(3,val-xp[k])
697 OPV(4,pow(val,xp[k]))
698 OPV(5,pow(xp[k],val))
699 default: ERROR(BAD_CODE);
700 }
701}
702
703int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) {
704 int k;
705 float val = *pval;
706 REQUIRES(xn == rn,BAD_SIZE);
707 DEBUGMSG("mapValF");
708 switch (code) {
709 OPV(0,val*xp[k])
710 OPV(1,val/xp[k])
711 OPV(2,val+xp[k])
712 OPV(3,val-xp[k])
713 OPV(4,pow(val,xp[k]))
714 OPV(5,pow(xp[k],val))
715 default: ERROR(BAD_CODE);
716 }
717}
718
719int mapValI(int code, int* pval, KIVEC(x), IVEC(r)) {
720 int k;
721 int val = *pval;
722 REQUIRES(xn == rn,BAD_SIZE);
723 DEBUGMSG("mapValI");
724 switch (code) {
725 OPV(0,val*xp[k])
726 OPV(1,val/xp[k])
727 OPV(2,val+xp[k])
728 OPV(3,val-xp[k])
729 OPV(6,mod(val,xp[k]))
730 OPV(7,mod(xp[k],val))
731 default: ERROR(BAD_CODE);
732 }
733}
734
735int mapValL(int code, int64_t* pval, KLVEC(x), LVEC(r)) {
736 int k;
737 int64_t val = *pval;
738 REQUIRES(xn == rn,BAD_SIZE);
739 DEBUGMSG("mapValL");
740 switch (code) {
741 OPV(0,val*xp[k])
742 OPV(1,val/xp[k])
743 OPV(2,val+xp[k])
744 OPV(3,val-xp[k])
745 OPV(6,mod_l(val,xp[k]))
746 OPV(7,mod_l(xp[k],val))
747 default: ERROR(BAD_CODE);
748 }
749}
750
751
752
753inline doublecomplex complex_add(doublecomplex a, doublecomplex b) {
754 doublecomplex r;
755 r.r = a.r+b.r;
756 r.i = a.i+b.i;
757 return r;
758}
759
760#define OPVb(C,E) case C: { for(k=0;k<xn;k++) r2p[k] = E; OK }
761int mapValC(int code, doublecomplex* pval, KCVEC(x), CVEC(r)) {
762 TCD* x2p = (TCD*)xp;
763 TCD* r2p = (TCD*)rp;
764 int k;
765 TCD val = * (TCD*)pval;
766 REQUIRES(xn == rn,BAD_SIZE);
767 DEBUGMSG("mapValC");
768 switch (code) {
769 OPVb(0,val*x2p[k])
770 OPVb(1,val/x2p[k])
771 OPVb(2,val+x2p[k])
772 OPVb(3,val-x2p[k])
773 OPVb(4,cpow(val,x2p[k]))
774 OPVb(5,cpow(x2p[k],val))
775 default: ERROR(BAD_CODE);
776 }
777}
778
779
780int mapValQ(int code, complex* pval, KQVEC(x), QVEC(r)) {
781 TCF* x2p = (TCF*)xp;
782 TCF* r2p = (TCF*)rp;
783 int k;
784 TCF val = *(TCF*)pval;
785 REQUIRES(xn == rn,BAD_SIZE);
786 DEBUGMSG("mapValQ");
787 switch (code) {
788 OPVb(0,val*x2p[k])
789 OPVb(1,val/x2p[k])
790 OPVb(2,val+x2p[k])
791 OPVb(3,val-x2p[k])
792 OPVb(4,cpow(val,x2p[k]))
793 OPVb(5,cpow(x2p[k],val))
794 default: ERROR(BAD_CODE);
795 }
796}
797
798
799
800#define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
801#define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
802#define OPZO(C,msg,O) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = ap[k] O bp[k]; OK }
803
804int zipR(int code, KDVEC(a), KDVEC(b), DVEC(r)) {
805REQUIRES(an == bn && an == rn, BAD_SIZE);
806 int k;
807 switch(code) {
808 OPZO(0,"zipR Add",+)
809 OPZO(1,"zipR Sub",-)
810 OPZO(2,"zipR Mul",*)
811 OPZO(3,"zipR Div",/)
812 OPZE(4,"zipR Pow", pow)
813 OPZE(5,"zipR ATan2",atan2)
814 default: ERROR(BAD_CODE);
815 }
816}
817
818int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) {
819REQUIRES(an == bn && an == rn, BAD_SIZE);
820 int k;
821 switch(code) {
822 OPZO(0,"zipR Add",+)
823 OPZO(1,"zipR Sub",-)
824 OPZO(2,"zipR Mul",*)
825 OPZO(3,"zipR Div",/)
826 OPZE(4,"zipR Pow", pow)
827 OPZE(5,"zipR ATan2",atan2)
828 default: ERROR(BAD_CODE);
829 }
830}
831
832
833int zipI(int code, KIVEC(a), KIVEC(b), IVEC(r)) {
834REQUIRES(an == bn && an == rn, BAD_SIZE);
835 int k;
836 switch(code) {
837 OPZO(0,"zipI Add",+)
838 OPZO(1,"zipI Sub",-)
839 OPZO(2,"zipI Mul",*)
840 OPZO(3,"zipI Div",/)
841 OPZO(6,"zipI Mod",%)
842 default: ERROR(BAD_CODE);
843 }
844}
845
846
847int zipL(int code, KLVEC(a), KLVEC(b), LVEC(r)) {
848REQUIRES(an == bn && an == rn, BAD_SIZE);
849 int k;
850 switch(code) {
851 OPZO(0,"zipI Add",+)
852 OPZO(1,"zipI Sub",-)
853 OPZO(2,"zipI Mul",*)
854 OPZO(3,"zipI Div",/)
855 OPZO(6,"zipI Mod",%)
856 default: ERROR(BAD_CODE);
857 }
858}
859
860
861#define OPZOb(C,msg,O) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) r2p[k] = a2p[k] O b2p[k]; OK }
862#define OPZEb(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) r2p[k] = E(a2p[k],b2p[k]); OK }
863int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
864 REQUIRES(an == bn && an == rn, BAD_SIZE);
865 TCD* a2p = (TCD*)ap;
866 TCD* b2p = (TCD*)bp;
867 TCD* r2p = (TCD*)rp;
868 int k;
869 switch(code) {
870 OPZOb(0,"zipC Add",+)
871 OPZOb(1,"zipC Sub",-)
872 OPZOb(2,"zipC Mul",*)
873 OPZOb(3,"zipC Div",/)
874 OPZEb(4,"zipC Pow",cpow)
875 default: ERROR(BAD_CODE);
876 }
877}
878
879
880
881
882
883#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_f_math_op(&E,ap[k],bp[k]); OK }
884
885int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
886 REQUIRES(an == bn && an == rn, BAD_SIZE);
887 TCF* a2p = (TCF*)ap;
888 TCF* b2p = (TCF*)bp;
889 TCF* r2p = (TCF*)rp;
890
891 int k;
892 switch(code) {
893 OPZOb(0,"zipC Add",+)
894 OPZOb(1,"zipC Sub",-)
895 OPZOb(2,"zipC Mul",*)
896 OPZOb(3,"zipC Div",/)
897 OPZEb(4,"zipC Pow",cpowf)
898 default: ERROR(BAD_CODE);
899 }
900}
901
902////////////////////////////////////////////////////////////////////////////////
903
904int vectorScan(char * file, int* n, double**pp){
905 FILE * fp;
906 fp = fopen (file, "r");
907 if(!fp) {
908 ERROR(BAD_FILE);
909 }
910 int nbuf = 100*100;
911 double * p = (double*)malloc(nbuf*sizeof(double));
912 int k=0;
913 double d;
914 int ok;
915 for (;;) {
916 ok = fscanf(fp,"%lf",&d);
917 if (ok<1) {
918 break;
919 }
920 if (k==nbuf) {
921 nbuf = nbuf * 2;
922 p = (double*)realloc(p,nbuf*sizeof(double));
923 // printf("R\n");
924 }
925 p[k++] = d;
926 }
927 *n = k;
928 *pp = p;
929 fclose(fp);
930 OK
931}
932
933////////////////////////////////////////////////////////////////////////////////
934
935#if defined (__APPLE__) || (__FreeBSD__)
936/* FreeBSD and Mac OS X do not provide random_r(), thread safety cannot be
937 guaranteed.
938 For FreeBSD and Mac OS X, nrand48() is much better than random().
939 See: http://www.evanjones.ca/random-thread-safe.html
940*/
941#pragma message "randomVector is not thread-safe in OSX and FreeBSD"
942#endif
943
944#if defined (__APPLE__) || (__FreeBSD__) || defined(_WIN32) || defined(WIN32)
945/* Windows use thread-safe random
946 See: http://stackoverflow.com/questions/143108/is-windows-rand-s-thread-safe
947*/
948inline double urandom() {
949 /* the probalility of matching will be theoretically p^3(in fact, it is not)
950 p is matching probalility of random().
951 using the test there, only 3 matches, using random(), 13783 matches
952 */
953 unsigned short state[3];
954 state[0] = random();
955 state[1] = random();
956 state[2] = random();
957
958 const long max_random = 2147483647; // 2**31 - 1
959 return (double)nrand48(state) / (double)max_random;
960}
961
962double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
963{
964 double V1=*pV1, V2=*pV2, S=*pS;
965 double X;
966
967 if(*phase == 0) {
968 do {
969 double U1 = urandom();
970 double U2 = urandom();
971
972 V1 = 2 * U1 - 1;
973 V2 = 2 * U2 - 1;
974 S = V1 * V1 + V2 * V2;
975 } while(S >= 1 || S == 0);
976
977 X = V1 * sqrt(-2 * log(S) / S);
978 } else
979 X = V2 * sqrt(-2 * log(S) / S);
980
981 *phase = 1 - *phase;
982 *pV1=V1; *pV2=V2; *pS=S;
983
984 return X;
985
986}
987
988int random_vector(unsigned int seed, int code, DVEC(r)) {
989 int phase = 0;
990 double V1,V2,S;
991
992 srandom(seed);
993
994 int k;
995 switch (code) {
996 case 0: { // uniform
997 for (k=0; k<rn; k++) {
998 rp[k] = urandom();
999 }
1000 OK
1001 }
1002 case 1: { // gaussian
1003 for (k=0; k<rn; k++) {
1004 rp[k] = gaussrand(&phase,&V1,&V2,&S);
1005 }
1006 OK
1007 }
1008
1009 default: ERROR(BAD_CODE);
1010 }
1011}
1012
1013#else
1014
1015inline double urandom(struct random_data * buffer) {
1016 int32_t res;
1017 random_r(buffer,&res);
1018 return (double)res/RAND_MAX;
1019}
1020
1021
1022// http://c-faq.com/lib/gaussian.html
1023double gaussrand(struct random_data *buffer,
1024 int *phase, double *pV1, double *pV2, double *pS)
1025{
1026 double V1=*pV1, V2=*pV2, S=*pS;
1027 double X;
1028
1029 if(*phase == 0) {
1030 do {
1031 double U1 = urandom(buffer);
1032 double U2 = urandom(buffer);
1033
1034 V1 = 2 * U1 - 1;
1035 V2 = 2 * U2 - 1;
1036 S = V1 * V1 + V2 * V2;
1037 } while(S >= 1 || S == 0);
1038
1039 X = V1 * sqrt(-2 * log(S) / S);
1040 } else
1041 X = V2 * sqrt(-2 * log(S) / S);
1042
1043 *phase = 1 - *phase;
1044 *pV1=V1; *pV2=V2; *pS=S;
1045
1046 return X;
1047
1048}
1049
1050int random_vector(unsigned int seed, int code, DVEC(r)) {
1051 struct random_data buffer;
1052 char random_state[128];
1053 memset(&buffer, 0, sizeof(struct random_data));
1054 memset(random_state, 0, sizeof(random_state));
1055
1056 initstate_r(seed,random_state,sizeof(random_state),&buffer);
1057 // setstate_r(random_state,&buffer);
1058 // srandom_r(seed,&buffer);
1059
1060 int phase = 0;
1061 double V1,V2,S;
1062
1063 int k;
1064 switch (code) {
1065 case 0: { // uniform
1066 for (k=0; k<rn; k++) {
1067 rp[k] = urandom(&buffer);
1068 }
1069 OK
1070 }
1071 case 1: { // gaussian
1072 for (k=0; k<rn; k++) {
1073 rp[k] = gaussrand(&buffer,&phase,&V1,&V2,&S);
1074 }
1075 OK
1076 }
1077
1078 default: ERROR(BAD_CODE);
1079 }
1080}
1081
1082#endif
1083
1084////////////////////////////////////////////////////////////////////////////////
1085
1086int
1087compare_doubles (const void *a, const void *b) {
1088 return *(double*)a > *(double*)b;
1089}
1090
1091int sort_valuesD(KDVEC(v),DVEC(r)) {
1092 memcpy(rp,vp,vn*sizeof(double));
1093 qsort(rp,rn,sizeof(double),compare_doubles);
1094 OK
1095}
1096
1097int
1098compare_floats (const void *a, const void *b) {
1099 return *(float*)a > *(float*)b;
1100}
1101
1102int sort_valuesF(KFVEC(v),FVEC(r)) {
1103 memcpy(rp,vp,vn*sizeof(float));
1104 qsort(rp,rn,sizeof(float),compare_floats);
1105 OK
1106}
1107
1108int
1109compare_ints(const void *a, const void *b) {
1110 return *(int*)a > *(int*)b;
1111}
1112
1113int sort_valuesI(KIVEC(v),IVEC(r)) {
1114 memcpy(rp,vp,vn*sizeof(int));
1115 qsort(rp,rn,sizeof(int),compare_ints);
1116 OK
1117}
1118
1119int
1120compare_longs(const void *a, const void *b) {
1121 return *(int64_t*)a > *(int64_t*)b;
1122}
1123
1124int sort_valuesL(KLVEC(v),LVEC(r)) {
1125 memcpy(rp,vp,vn*sizeof(int64_t));
1126 qsort(rp,rn,sizeof(int64_t),compare_ints);
1127 OK
1128}
1129
1130
1131////////////////////////////////////////
1132
1133
1134#define SORTIDX_IMP(T,C) \
1135 T* x = (T*)malloc(sizeof(T)*vn); \
1136 int k; \
1137 for (k=0;k<vn;k++) { \
1138 x[k].pos = k; \
1139 x[k].val = vp[k]; \
1140 } \
1141 \
1142 qsort(x,vn,sizeof(T),C); \
1143 \
1144 for (k=0;k<vn;k++) { \
1145 rp[k] = x[k].pos; \
1146 } \
1147 free(x); \
1148 OK
1149
1150
1151typedef struct DI { int pos; double val;} DI;
1152
1153int compare_doubles_i (const void *a, const void *b) {
1154 return ((DI*)a)->val > ((DI*)b)->val;
1155}
1156
1157int sort_indexD(KDVEC(v),IVEC(r)) {
1158 SORTIDX_IMP(DI,compare_doubles_i)
1159}
1160
1161
1162typedef struct FI { int pos; float val;} FI;
1163
1164int compare_floats_i (const void *a, const void *b) {
1165 return ((FI*)a)->val > ((FI*)b)->val;
1166}
1167
1168int sort_indexF(KFVEC(v),IVEC(r)) {
1169 SORTIDX_IMP(FI,compare_floats_i)
1170}
1171
1172
1173typedef struct II { int pos; int val;} II;
1174
1175int compare_ints_i (const void *a, const void *b) {
1176 return ((II*)a)->val > ((II*)b)->val;
1177}
1178
1179int sort_indexI(KIVEC(v),IVEC(r)) {
1180 SORTIDX_IMP(II,compare_ints_i)
1181}
1182
1183
1184typedef struct LI { int pos; int64_t val;} LI;
1185
1186int compare_longs_i (const void *a, const void *b) {
1187 return ((II*)a)->val > ((II*)b)->val;
1188}
1189
1190int sort_indexL(KLVEC(v),LVEC(r)) {
1191 SORTIDX_IMP(II,compare_longs_i)
1192}
1193
1194
1195////////////////////////////////////////////////////////////////////////////////
1196
1197int round_vector(KDVEC(v),DVEC(r)) {
1198 int k;
1199 for(k=0; k<vn; k++) {
1200 rp[k] = round(vp[k]);
1201 }
1202 OK
1203}
1204
1205////////////////////////////////////////////////////////////////////////////////
1206
1207int round_vector_i(KDVEC(v),IVEC(r)) {
1208 int k;
1209 for(k=0; k<vn; k++) {
1210 rp[k] = round(vp[k]);
1211 }
1212 OK
1213}
1214
1215
1216int mod_vector(int m, KIVEC(v), IVEC(r)) {
1217 int k;
1218 for(k=0; k<vn; k++) {
1219 rp[k] = mod(vp[k],m);
1220 }
1221 OK
1222}
1223
1224int div_vector(int m, KIVEC(v), IVEC(r)) {
1225 int k;
1226 for(k=0; k<vn; k++) {
1227 rp[k] = vp[k] / m;
1228 }
1229 OK
1230}
1231
1232int range_vector(IVEC(r)) {
1233 int k;
1234 for(k=0; k<rn; k++) {
1235 rp[k] = k;
1236 }
1237 OK
1238}
1239
1240///////////////////////////
1241
1242
1243int round_vector_l(KDVEC(v),LVEC(r)) {
1244 int k;
1245 for(k=0; k<vn; k++) {
1246 rp[k] = round(vp[k]);
1247 }
1248 OK
1249}
1250
1251
1252int mod_vector_l(int64_t m, KLVEC(v), LVEC(r)) {
1253 int k;
1254 for(k=0; k<vn; k++) {
1255 rp[k] = mod_l(vp[k],m);
1256 }
1257 OK
1258}
1259
1260int div_vector_l(int64_t m, KLVEC(v), LVEC(r)) {
1261 int k;
1262 for(k=0; k<vn; k++) {
1263 rp[k] = vp[k] / m;
1264 }
1265 OK
1266}
1267
1268int range_vector_l(LVEC(r)) {
1269 int k;
1270 for(k=0; k<rn; k++) {
1271 rp[k] = k;
1272 }
1273 OK
1274}
1275
1276
1277
1278//////////////////// constant /////////////////////////
1279
1280int constantF(float * pval, FVEC(r)) {
1281 DEBUGMSG("constantF")
1282 int k;
1283 double val = *pval;
1284 for(k=0;k<rn;k++) {
1285 rp[k]=val;
1286 }
1287 OK
1288}
1289
1290int constantR(double * pval, DVEC(r)) {
1291 DEBUGMSG("constantR")
1292 int k;
1293 double val = *pval;
1294 for(k=0;k<rn;k++) {
1295 rp[k]=val;
1296 }
1297 OK
1298}
1299
1300int constantQ(complex* pval, QVEC(r)) {
1301 DEBUGMSG("constantQ")
1302 int k;
1303 complex val = *pval;
1304 for(k=0;k<rn;k++) {
1305 rp[k]=val;
1306 }
1307 OK
1308}
1309
1310int constantC(doublecomplex* pval, CVEC(r)) {
1311 DEBUGMSG("constantC")
1312 int k;
1313 doublecomplex val = *pval;
1314 for(k=0;k<rn;k++) {
1315 rp[k]=val;
1316 }
1317 OK
1318}
1319
1320
1321
1322int constantI(int * pval, IVEC(r)) {
1323 DEBUGMSG("constantI")
1324 int k;
1325 int val = *pval;
1326 for(k=0;k<rn;k++) {
1327 rp[k]=val;
1328 }
1329 OK
1330}
1331
1332
1333
1334int constantL(int64_t * pval, LVEC(r)) {
1335 DEBUGMSG("constantL")
1336 int k;
1337 int64_t val = *pval;
1338 for(k=0;k<rn;k++) {
1339 rp[k]=val;
1340 }
1341 OK
1342}
1343
1344
1345//////////////////// type conversions /////////////////////////
1346
1347#define CONVERT_IMP { \
1348 int k; \
1349 for(k=0;k<xn;k++) { \
1350 yp[k]=xp[k]; \
1351 } \
1352 OK }
1353
1354int float2double(FVEC(x),DVEC(y)) CONVERT_IMP
1355
1356int float2int(KFVEC(x),IVEC(y)) CONVERT_IMP
1357
1358int double2float(DVEC(x),FVEC(y)) CONVERT_IMP
1359
1360int double2int(KDVEC(x),IVEC(y)) CONVERT_IMP
1361
1362int double2long(KDVEC(x),LVEC(y)) CONVERT_IMP
1363
1364int int2float(KIVEC(x),FVEC(y)) CONVERT_IMP
1365
1366int int2double(KIVEC(x),DVEC(y)) CONVERT_IMP
1367
1368int int2long(KIVEC(x),LVEC(y)) CONVERT_IMP
1369
1370int long2int(KLVEC(x),IVEC(y)) CONVERT_IMP
1371
1372int long2double(KLVEC(x),DVEC(y)) CONVERT_IMP
1373
1374
1375//////////////////// conjugate /////////////////////////
1376
1377int conjugateQ(KQVEC(x),QVEC(t)) {
1378 REQUIRES(xn==tn,BAD_SIZE);
1379 DEBUGMSG("conjugateQ");
1380 int k;
1381 for(k=0;k<xn;k++) {
1382 tp[k].r = xp[k].r;
1383 tp[k].i = -xp[k].i;
1384 }
1385 OK
1386}
1387
1388int conjugateC(KCVEC(x),CVEC(t)) {
1389 REQUIRES(xn==tn,BAD_SIZE);
1390 DEBUGMSG("conjugateC");
1391 int k;
1392 for(k=0;k<xn;k++) {
1393 tp[k].r = xp[k].r;
1394 tp[k].i = -xp[k].i;
1395 }
1396 OK
1397}
1398
1399//////////////////// step /////////////////////////
1400
1401#define STEP_IMP \
1402 int k; \
1403 for(k=0;k<xn;k++) { \
1404 yp[k]=xp[k]>0; \
1405 } \
1406 OK
1407
1408int stepF(KFVEC(x),FVEC(y)) {
1409 STEP_IMP
1410}
1411
1412int stepD(KDVEC(x),DVEC(y)) {
1413 STEP_IMP
1414}
1415
1416int stepI(KIVEC(x),IVEC(y)) {
1417 STEP_IMP
1418}
1419
1420int stepL(KLVEC(x),LVEC(y)) {
1421 STEP_IMP
1422}
1423
1424
1425//////////////////// cond /////////////////////////
1426
1427#define COMPARE_IMP \
1428 REQUIRES(xn==yn && xn==rn ,BAD_SIZE); \
1429 int k; \
1430 for(k=0;k<xn;k++) { \
1431 rp[k] = xp[k]<yp[k]?-1:(xp[k]>yp[k]?1:0); \
1432 } \
1433 OK
1434
1435
1436int compareF(KFVEC(x),KFVEC(y),IVEC(r)) {
1437 COMPARE_IMP
1438}
1439
1440int compareD(KDVEC(x),KDVEC(y),IVEC(r)) {
1441 COMPARE_IMP
1442}
1443
1444int compareI(KIVEC(x),KIVEC(y),IVEC(r)) {
1445 COMPARE_IMP
1446}
1447
1448int compareL(KLVEC(x),KLVEC(y),IVEC(r)) {
1449 COMPARE_IMP
1450}
1451
1452
1453
1454#define CHOOSE_IMP \
1455 REQUIRES(condn==ltn && ltn==eqn && ltn==gtn && ltn==rn ,BAD_SIZE); \
1456 int k; \
1457 for(k=0;k<condn;k++) { \
1458 rp[k] = condp[k]<0?ltp[k]:(condp[k]>0?gtp[k]:eqp[k]); \
1459 } \
1460 OK
1461
1462int chooseF(KIVEC(cond),KFVEC(lt),KFVEC(eq),KFVEC(gt),FVEC(r)) {
1463 CHOOSE_IMP
1464}
1465
1466int chooseD(KIVEC(cond),KDVEC(lt),KDVEC(eq),KDVEC(gt),DVEC(r)) {
1467 CHOOSE_IMP
1468}
1469
1470int chooseI(KIVEC(cond),KIVEC(lt),KIVEC(eq),KIVEC(gt),IVEC(r)) {
1471 CHOOSE_IMP
1472}
1473
1474int chooseL(KIVEC(cond),KLVEC(lt),KLVEC(eq),KLVEC(gt),LVEC(r)) {
1475 CHOOSE_IMP
1476}
1477
1478
1479int chooseC(KIVEC(cond),KCVEC(lt),KCVEC(eq),KCVEC(gt),CVEC(r)) {
1480 CHOOSE_IMP
1481}
1482
1483int chooseQ(KIVEC(cond),KQVEC(lt),KQVEC(eq),KQVEC(gt),QVEC(r)) {
1484 CHOOSE_IMP
1485}
1486
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-28 20:55:47 +0000