summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
Diffstat (limited to 'packages')
-rw-r--r--packages/base/hmatrix-base.cabal30
-rw-r--r--packages/base/src/C/lapack-aux.c1489
-rw-r--r--packages/base/src/C/lapack-aux.h60
-rw-r--r--packages/base/src/Data/Packed.hs (renamed from packages/hmatrix/src/Data/Packed.hs)20
-rw-r--r--packages/base/src/Data/Packed/Development.hs (renamed from packages/hmatrix/src/Data/Packed/Development.hs)7
-rw-r--r--packages/base/src/Data/Packed/Foreign.hs (renamed from packages/hmatrix/src/Data/Packed/Foreign.hs)1
-rw-r--r--packages/base/src/Data/Packed/Internal.hs (renamed from packages/hmatrix/src/Data/Packed/Internal.hs)0
-rw-r--r--packages/base/src/Data/Packed/Internal/Common.hs (renamed from packages/hmatrix/src/Data/Packed/Internal/Common.hs)27
-rw-r--r--packages/base/src/Data/Packed/Internal/Matrix.hs (renamed from packages/hmatrix/src/Data/Packed/Internal/Matrix.hs)75
-rw-r--r--packages/base/src/Data/Packed/Internal/Signatures.hs (renamed from packages/hmatrix/src/Data/Packed/Internal/Signatures.hs)10
-rw-r--r--packages/base/src/Data/Packed/Internal/Vector.hs (renamed from packages/hmatrix/src/Data/Packed/Internal/Vector.hs)56
-rw-r--r--packages/base/src/Data/Packed/Matrix.hs (renamed from packages/hmatrix/src/Data/Packed/Matrix.hs)0
-rw-r--r--packages/base/src/Data/Packed/ST.hs (renamed from packages/hmatrix/src/Data/Packed/ST.hs)7
-rw-r--r--packages/base/src/Data/Packed/Vector.hs (renamed from packages/hmatrix/src/Data/Packed/Vector.hs)0
14 files changed, 1606 insertions, 176 deletions
diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal
index 148c8f3..96f90e2 100644
--- a/packages/base/hmatrix-base.cabal
+++ b/packages/base/hmatrix-base.cabal
@@ -15,22 +15,42 @@ cabal-version: >=1.8
15 15
16build-type: Simple 16build-type: Simple
17 17
18extra-source-files: 18extra-source-files: src/C/lapack-aux.h
19 19
20library 20library
21 21
22 Build-Depends: base >= 4 && < 5 22 Build-Depends: base,
23 binary,
24 array,
25 deepseq,
26 storable-complex,
27 vector >= 0.8
23 28
24 hs-source-dirs: src 29 hs-source-dirs: src
25 30
26 exposed-modules: 31 exposed-modules: Data.Packed,
27 32 Data.Packed.Vector,
28 other-modules: 33 Data.Packed.Matrix,
34 Data.Packed.Foreign,
35 Data.Packed.ST,
36 Data.Packed.Development
37
38 other-modules: Data.Packed.Internal,
39 Data.Packed.Internal.Common,
40 Data.Packed.Internal.Signatures,
41 Data.Packed.Internal.Vector,
42 Data.Packed.Internal.Matrix
43
44 C-sources: src/C/lapack-aux.c
45
29 46
30 extensions: ForeignFunctionInterface, 47 extensions: ForeignFunctionInterface,
31 CPP 48 CPP
32 49
33 ghc-options: -Wall 50 ghc-options: -Wall
51 -fno-warn-missing-signatures
52 -fno-warn-orphans
53-- -fno-warn-unused-binds
34 54
35 cc-options: -O4 -msse2 -Wall 55 cc-options: -O4 -msse2 -Wall
36 56
diff --git a/packages/base/src/C/lapack-aux.c b/packages/base/src/C/lapack-aux.c
new file mode 100644
index 0000000..e5e45ef
--- /dev/null
+++ b/packages/base/src/C/lapack-aux.c
@@ -0,0 +1,1489 @@
1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4#include <math.h>
5#include <time.h>
6#include "lapack-aux.h"
7
8#define MACRO(B) do {B} while (0)
9#define ERROR(CODE) MACRO(return CODE;)
10#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
11
12#define MIN(A,B) ((A)<(B)?(A):(B))
13#define MAX(A,B) ((A)>(B)?(A):(B))
14
15// #define DBGL
16
17#ifdef DBGL
18#define DEBUGMSG(M) printf("\nLAPACK "M"\n");
19#else
20#define DEBUGMSG(M)
21#endif
22
23#define OK return 0;
24
25// #ifdef DBGL
26// #define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL);
27// #define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;);
28// #else
29// #define DEBUGMSG(M)
30// #define OK return 0;
31// #endif
32
33#define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \
34 for(q=0;q<M##r*M##c;q++) printf("%.1f ",M##p[q]); printf("\n");}
35
36#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
37
38#define BAD_SIZE 2000
39#define BAD_CODE 2001
40#define MEM 2002
41#define BAD_FILE 2003
42#define SINGULAR 2004
43#define NOCONVER 2005
44#define NODEFPOS 2006
45#define NOSPRTD 2007
46
47//---------------------------------------
48void asm_finit() {
49#ifdef i386
50
51// asm("finit");
52
53 static unsigned char buf[108];
54 asm("FSAVE %0":"=m" (buf));
55
56 #if FPUDEBUG
57 if(buf[8]!=255 || buf[9]!=255) { // print warning in red
58 printf("%c[;31mWarning: FPU TAG = %x %x\%c[0m\n",0x1B,buf[8],buf[9],0x1B);
59 }
60 #endif
61
62 #if NANDEBUG
63 asm("FRSTOR %0":"=m" (buf));
64 #endif
65
66#endif
67}
68
69//---------------------------------------
70
71#if NANDEBUG
72
73#define CHECKNANR(M,msg) \
74{ int k; \
75for(k=0; k<(M##r * M##c); k++) { \
76 if(M##p[k] != M##p[k]) { \
77 printf(msg); \
78 TRACEMAT(M) \
79 /*exit(1);*/ \
80 } \
81} \
82}
83
84#define CHECKNANC(M,msg) \
85{ int k; \
86for(k=0; k<(M##r * M##c); k++) { \
87 if( M##p[k].r != M##p[k].r \
88 || M##p[k].i != M##p[k].i) { \
89 printf(msg); \
90 /*exit(1);*/ \
91 } \
92} \
93}
94
95#else
96#define CHECKNANC(M,msg)
97#define CHECKNANR(M,msg)
98#endif
99
100//---------------------------------------
101
102//////////////////// real svd ////////////////////////////////////
103
104/* Subroutine */ int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
105 doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
106 ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
107 integer *info);
108
109int svd_l_R(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
110 integer m = ar;
111 integer n = ac;
112 integer q = MIN(m,n);
113 REQUIRES(sn==q,BAD_SIZE);
114 REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE);
115 char* jobu = "A";
116 if (up==NULL) {
117 jobu = "N";
118 } else {
119 if (uc==q) {
120 jobu = "S";
121 }
122 }
123 REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE);
124 char* jobvt = "A";
125 integer ldvt = n;
126 if (vp==NULL) {
127 jobvt = "N";
128 } else {
129 if (vr==q) {
130 jobvt = "S";
131 ldvt = q;
132 }
133 }
134 DEBUGMSG("svd_l_R");
135 double *B = (double*)malloc(m*n*sizeof(double));
136 CHECK(!B,MEM);
137 memcpy(B,ap,m*n*sizeof(double));
138 integer lwork = -1;
139 integer res;
140 // ask for optimal lwork
141 double ans;
142 dgesvd_ (jobu,jobvt,
143 &m,&n,B,&m,
144 sp,
145 up,&m,
146 vp,&ldvt,
147 &ans, &lwork,
148 &res);
149 lwork = ceil(ans);
150 double * work = (double*)malloc(lwork*sizeof(double));
151 CHECK(!work,MEM);
152 dgesvd_ (jobu,jobvt,
153 &m,&n,B,&m,
154 sp,
155 up,&m,
156 vp,&ldvt,
157 work, &lwork,
158 &res);
159 CHECK(res,res);
160 free(work);
161 free(B);
162 OK
163}
164
165// (alternative version)
166
167/* Subroutine */ int dgesdd_(char *jobz, integer *m, integer *n, doublereal *
168 a, integer *lda, doublereal *s, doublereal *u, integer *ldu,
169 doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
170 integer *iwork, integer *info);
171
172int svd_l_Rdd(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
173 integer m = ar;
174 integer n = ac;
175 integer q = MIN(m,n);
176 REQUIRES(sn==q,BAD_SIZE);
177 REQUIRES((up == NULL && vp == NULL)
178 || (ur==m && vc==n
179 && ((uc == q && vr == q)
180 || (uc == m && vc==n))),BAD_SIZE);
181 char* jobz = "A";
182 integer ldvt = n;
183 if (up==NULL) {
184 jobz = "N";
185 } else {
186 if (uc==q && vr == q) {
187 jobz = "S";
188 ldvt = q;
189 }
190 }
191 DEBUGMSG("svd_l_Rdd");
192 double *B = (double*)malloc(m*n*sizeof(double));
193 CHECK(!B,MEM);
194 memcpy(B,ap,m*n*sizeof(double));
195 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
196 CHECK(!iwk,MEM);
197 integer lwk = -1;
198 integer res;
199 // ask for optimal lwk
200 double ans;
201 dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res);
202 lwk = ans;
203 double * workv = (double*)malloc(lwk*sizeof(double));
204 CHECK(!workv,MEM);
205 dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res);
206 CHECK(res,res);
207 free(iwk);
208 free(workv);
209 free(B);
210 OK
211}
212
213//////////////////// complex svd ////////////////////////////////////
214
215// not in clapack.h
216
217int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
218 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
219 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
220 integer *lwork, doublereal *rwork, integer *info);
221
222int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
223 integer m = ar;
224 integer n = ac;
225 integer q = MIN(m,n);
226 REQUIRES(sn==q,BAD_SIZE);
227 REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE);
228 char* jobu = "A";
229 if (up==NULL) {
230 jobu = "N";
231 } else {
232 if (uc==q) {
233 jobu = "S";
234 }
235 }
236 REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE);
237 char* jobvt = "A";
238 integer ldvt = n;
239 if (vp==NULL) {
240 jobvt = "N";
241 } else {
242 if (vr==q) {
243 jobvt = "S";
244 ldvt = q;
245 }
246 }DEBUGMSG("svd_l_C");
247 doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
248 CHECK(!B,MEM);
249 memcpy(B,ap,m*n*sizeof(doublecomplex));
250
251 double *rwork = (double*) malloc(5*q*sizeof(double));
252 CHECK(!rwork,MEM);
253 integer lwork = -1;
254 integer res;
255 // ask for optimal lwork
256 doublecomplex ans;
257 zgesvd_ (jobu,jobvt,
258 &m,&n,B,&m,
259 sp,
260 up,&m,
261 vp,&ldvt,
262 &ans, &lwork,
263 rwork,
264 &res);
265 lwork = ceil(ans.r);
266 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
267 CHECK(!work,MEM);
268 zgesvd_ (jobu,jobvt,
269 &m,&n,B,&m,
270 sp,
271 up,&m,
272 vp,&ldvt,
273 work, &lwork,
274 rwork,
275 &res);
276 CHECK(res,res);
277 free(work);
278 free(rwork);
279 free(B);
280 OK
281}
282
283int zgesdd_ (char *jobz, integer *m, integer *n,
284 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
285 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
286 integer *lwork, doublereal *rwork, integer* iwork, integer *info);
287
288int svd_l_Cdd(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
289 //printf("entro\n");
290 integer m = ar;
291 integer n = ac;
292 integer q = MIN(m,n);
293 REQUIRES(sn==q,BAD_SIZE);
294 REQUIRES((up == NULL && vp == NULL)
295 || (ur==m && vc==n
296 && ((uc == q && vr == q)
297 || (uc == m && vc==n))),BAD_SIZE);
298 char* jobz = "A";
299 integer ldvt = n;
300 if (up==NULL) {
301 jobz = "N";
302 } else {
303 if (uc==q && vr == q) {
304 jobz = "S";
305 ldvt = q;
306 }
307 }
308 DEBUGMSG("svd_l_Cdd");
309 doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
310 CHECK(!B,MEM);
311 memcpy(B,ap,m*n*sizeof(doublecomplex));
312 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
313 CHECK(!iwk,MEM);
314 int lrwk;
315 if (0 && *jobz == 'N') {
316 lrwk = 5*q; // does not work, crash at free below
317 } else {
318 lrwk = 5*q*q + 7*q;
319 }
320 double *rwk = (double*)malloc(lrwk*sizeof(double));;
321 CHECK(!rwk,MEM);
322 //printf("%s %ld %d\n",jobz,q,lrwk);
323 integer lwk = -1;
324 integer res;
325 // ask for optimal lwk
326 doublecomplex ans;
327 zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
328 lwk = ans.r;
329 //printf("lwk = %ld\n",lwk);
330 doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));
331 CHECK(!workv,MEM);
332 zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res);
333 //printf("res = %ld\n",res);
334 CHECK(res,res);
335 free(workv); // printf("freed workv\n");
336 free(rwk); // printf("freed rwk\n");
337 free(iwk); // printf("freed iwk\n");
338 free(B); // printf("freed B, salgo\n");
339 OK
340}
341
342//////////////////// general complex eigensystem ////////////
343
344/* Subroutine */ int zgeev_(char *jobvl, char *jobvr, integer *n,
345 doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl,
346 integer *ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work,
347 integer *lwork, doublereal *rwork, integer *info);
348
349int eig_l_C(KCMAT(a), CMAT(u), CVEC(s),CMAT(v)) {
350 integer n = ar;
351 REQUIRES(ac==n && sn==n, BAD_SIZE);
352 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
353 char jobvl = up==NULL?'N':'V';
354 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
355 char jobvr = vp==NULL?'N':'V';
356 DEBUGMSG("eig_l_C");
357 doublecomplex *B = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
358 CHECK(!B,MEM);
359 memcpy(B,ap,n*n*sizeof(doublecomplex));
360 double *rwork = (double*) malloc(2*n*sizeof(double));
361 CHECK(!rwork,MEM);
362 integer lwork = -1;
363 integer res;
364 // ask for optimal lwork
365 doublecomplex ans;
366 //printf("ask zgeev\n");
367 zgeev_ (&jobvl,&jobvr,
368 &n,B,&n,
369 sp,
370 up,&n,
371 vp,&n,
372 &ans, &lwork,
373 rwork,
374 &res);
375 lwork = ceil(ans.r);
376 //printf("ans = %d\n",lwork);
377 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
378 CHECK(!work,MEM);
379 //printf("zgeev\n");
380 zgeev_ (&jobvl,&jobvr,
381 &n,B,&n,
382 sp,
383 up,&n,
384 vp,&n,
385 work, &lwork,
386 rwork,
387 &res);
388 CHECK(res,res);
389 free(work);
390 free(rwork);
391 free(B);
392 OK
393}
394
395
396
397//////////////////// general real eigensystem ////////////
398
399/* Subroutine */ int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *
400 a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl,
401 integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
402 integer *lwork, integer *info);
403
404int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) {
405 integer n = ar;
406 REQUIRES(ac==n && sn==n, BAD_SIZE);
407 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
408 char jobvl = up==NULL?'N':'V';
409 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
410 char jobvr = vp==NULL?'N':'V';
411 DEBUGMSG("eig_l_R");
412 double *B = (double*)malloc(n*n*sizeof(double));
413 CHECK(!B,MEM);
414 memcpy(B,ap,n*n*sizeof(double));
415 integer lwork = -1;
416 integer res;
417 // ask for optimal lwork
418 double ans;
419 //printf("ask dgeev\n");
420 dgeev_ (&jobvl,&jobvr,
421 &n,B,&n,
422 (double*)sp, (double*)sp+n,
423 up,&n,
424 vp,&n,
425 &ans, &lwork,
426 &res);
427 lwork = ceil(ans);
428 //printf("ans = %d\n",lwork);
429 double * work = (double*)malloc(lwork*sizeof(double));
430 CHECK(!work,MEM);
431 //printf("dgeev\n");
432 dgeev_ (&jobvl,&jobvr,
433 &n,B,&n,
434 (double*)sp, (double*)sp+n,
435 up,&n,
436 vp,&n,
437 work, &lwork,
438 &res);
439 CHECK(res,res);
440 free(work);
441 free(B);
442 OK
443}
444
445
446//////////////////// symmetric real eigensystem ////////////
447
448/* Subroutine */ int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a,
449 integer *lda, doublereal *w, doublereal *work, integer *lwork,
450 integer *info);
451
452int eig_l_S(int wantV,KDMAT(a),DVEC(s),DMAT(v)) {
453 integer n = ar;
454 REQUIRES(ac==n && sn==n, BAD_SIZE);
455 REQUIRES(vr==n && vc==n, BAD_SIZE);
456 char jobz = wantV?'V':'N';
457 DEBUGMSG("eig_l_S");
458 memcpy(vp,ap,n*n*sizeof(double));
459 integer lwork = -1;
460 char uplo = 'U';
461 integer res;
462 // ask for optimal lwork
463 double ans;
464 //printf("ask dsyev\n");
465 dsyev_ (&jobz,&uplo,
466 &n,vp,&n,
467 sp,
468 &ans, &lwork,
469 &res);
470 lwork = ceil(ans);
471 //printf("ans = %d\n",lwork);
472 double * work = (double*)malloc(lwork*sizeof(double));
473 CHECK(!work,MEM);
474 dsyev_ (&jobz,&uplo,
475 &n,vp,&n,
476 sp,
477 work, &lwork,
478 &res);
479 CHECK(res,res);
480 free(work);
481 OK
482}
483
484//////////////////// hermitian complex eigensystem ////////////
485
486/* Subroutine */ int zheev_(char *jobz, char *uplo, integer *n, doublecomplex
487 *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork,
488 doublereal *rwork, integer *info);
489
490int eig_l_H(int wantV,KCMAT(a),DVEC(s),CMAT(v)) {
491 integer n = ar;
492 REQUIRES(ac==n && sn==n, BAD_SIZE);
493 REQUIRES(vr==n && vc==n, BAD_SIZE);
494 char jobz = wantV?'V':'N';
495 DEBUGMSG("eig_l_H");
496 memcpy(vp,ap,2*n*n*sizeof(double));
497 double *rwork = (double*) malloc((3*n-2)*sizeof(double));
498 CHECK(!rwork,MEM);
499 integer lwork = -1;
500 char uplo = 'U';
501 integer res;
502 // ask for optimal lwork
503 doublecomplex ans;
504 //printf("ask zheev\n");
505 zheev_ (&jobz,&uplo,
506 &n,vp,&n,
507 sp,
508 &ans, &lwork,
509 rwork,
510 &res);
511 lwork = ceil(ans.r);
512 //printf("ans = %d\n",lwork);
513 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
514 CHECK(!work,MEM);
515 zheev_ (&jobz,&uplo,
516 &n,vp,&n,
517 sp,
518 work, &lwork,
519 rwork,
520 &res);
521 CHECK(res,res);
522 free(work);
523 free(rwork);
524 OK
525}
526
527//////////////////// general real linear system ////////////
528
529/* Subroutine */ int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
530 *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info);
531
532int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
533 integer n = ar;
534 integer nhrs = bc;
535 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
536 DEBUGMSG("linearSolveR_l");
537 double*AC = (double*)malloc(n*n*sizeof(double));
538 memcpy(AC,ap,n*n*sizeof(double));
539 memcpy(xp,bp,n*nhrs*sizeof(double));
540 integer * ipiv = (integer*)malloc(n*sizeof(integer));
541 integer res;
542 dgesv_ (&n,&nhrs,
543 AC, &n,
544 ipiv,
545 xp, &n,
546 &res);
547 if(res>0) {
548 return SINGULAR;
549 }
550 CHECK(res,res);
551 free(ipiv);
552 free(AC);
553 OK
554}
555
556//////////////////// general complex linear system ////////////
557
558/* Subroutine */ int zgesv_(integer *n, integer *nrhs, doublecomplex *a,
559 integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer *
560 info);
561
562int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
563 integer n = ar;
564 integer nhrs = bc;
565 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
566 DEBUGMSG("linearSolveC_l");
567 doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
568 memcpy(AC,ap,n*n*sizeof(doublecomplex));
569 memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));
570 integer * ipiv = (integer*)malloc(n*sizeof(integer));
571 integer res;
572 zgesv_ (&n,&nhrs,
573 AC, &n,
574 ipiv,
575 xp, &n,
576 &res);
577 if(res>0) {
578 return SINGULAR;
579 }
580 CHECK(res,res);
581 free(ipiv);
582 free(AC);
583 OK
584}
585
586//////// symmetric positive definite real linear system using Cholesky ////////////
587
588/* Subroutine */ int dpotrs_(char *uplo, integer *n, integer *nrhs,
589 doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *
590 info);
591
592int cholSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
593 integer n = ar;
594 integer nhrs = bc;
595 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
596 DEBUGMSG("cholSolveR_l");
597 memcpy(xp,bp,n*nhrs*sizeof(double));
598 integer res;
599 dpotrs_ ("U",
600 &n,&nhrs,
601 (double*)ap, &n,
602 xp, &n,
603 &res);
604 CHECK(res,res);
605 OK
606}
607
608//////// Hermitian positive definite real linear system using Cholesky ////////////
609
610/* Subroutine */ int zpotrs_(char *uplo, integer *n, integer *nrhs,
611 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
612 integer *info);
613
614int cholSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
615 integer n = ar;
616 integer nhrs = bc;
617 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
618 DEBUGMSG("cholSolveC_l");
619 memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));
620 integer res;
621 zpotrs_ ("U",
622 &n,&nhrs,
623 (doublecomplex*)ap, &n,
624 xp, &n,
625 &res);
626 CHECK(res,res);
627 OK
628}
629
630//////////////////// least squares real linear system ////////////
631
632/* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer *
633 nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb,
634 doublereal *work, integer *lwork, integer *info);
635
636int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
637 integer m = ar;
638 integer n = ac;
639 integer nrhs = bc;
640 integer ldb = xr;
641 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
642 DEBUGMSG("linearSolveLSR_l");
643 double*AC = (double*)malloc(m*n*sizeof(double));
644 memcpy(AC,ap,m*n*sizeof(double));
645 if (m>=n) {
646 memcpy(xp,bp,m*nrhs*sizeof(double));
647 } else {
648 int k;
649 for(k = 0; k<nrhs; k++) {
650 memcpy(xp+ldb*k,bp+m*k,m*sizeof(double));
651 }
652 }
653 integer res;
654 integer lwork = -1;
655 double ans;
656 dgels_ ("N",&m,&n,&nrhs,
657 AC,&m,
658 xp,&ldb,
659 &ans,&lwork,
660 &res);
661 lwork = ceil(ans);
662 //printf("ans = %d\n",lwork);
663 double * work = (double*)malloc(lwork*sizeof(double));
664 dgels_ ("N",&m,&n,&nrhs,
665 AC,&m,
666 xp,&ldb,
667 work,&lwork,
668 &res);
669 if(res>0) {
670 return SINGULAR;
671 }
672 CHECK(res,res);
673 free(work);
674 free(AC);
675 OK
676}
677
678//////////////////// least squares complex linear system ////////////
679
680/* Subroutine */ int zgels_(char *trans, integer *m, integer *n, integer *
681 nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
682 doublecomplex *work, integer *lwork, integer *info);
683
684int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
685 integer m = ar;
686 integer n = ac;
687 integer nrhs = bc;
688 integer ldb = xr;
689 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
690 DEBUGMSG("linearSolveLSC_l");
691 doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
692 memcpy(AC,ap,m*n*sizeof(doublecomplex));
693 if (m>=n) {
694 memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));
695 } else {
696 int k;
697 for(k = 0; k<nrhs; k++) {
698 memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));
699 }
700 }
701 integer res;
702 integer lwork = -1;
703 doublecomplex ans;
704 zgels_ ("N",&m,&n,&nrhs,
705 AC,&m,
706 xp,&ldb,
707 &ans,&lwork,
708 &res);
709 lwork = ceil(ans.r);
710 //printf("ans = %d\n",lwork);
711 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
712 zgels_ ("N",&m,&n,&nrhs,
713 AC,&m,
714 xp,&ldb,
715 work,&lwork,
716 &res);
717 if(res>0) {
718 return SINGULAR;
719 }
720 CHECK(res,res);
721 free(work);
722 free(AC);
723 OK
724}
725
726//////////////////// least squares real linear system using SVD ////////////
727
728/* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs,
729 doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *
730 s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork,
731 integer *info);
732
733int linearSolveSVDR_l(double rcond,KDMAT(a),KDMAT(b),DMAT(x)) {
734 integer m = ar;
735 integer n = ac;
736 integer nrhs = bc;
737 integer ldb = xr;
738 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
739 DEBUGMSG("linearSolveSVDR_l");
740 double*AC = (double*)malloc(m*n*sizeof(double));
741 double*S = (double*)malloc(MIN(m,n)*sizeof(double));
742 memcpy(AC,ap,m*n*sizeof(double));
743 if (m>=n) {
744 memcpy(xp,bp,m*nrhs*sizeof(double));
745 } else {
746 int k;
747 for(k = 0; k<nrhs; k++) {
748 memcpy(xp+ldb*k,bp+m*k,m*sizeof(double));
749 }
750 }
751 integer res;
752 integer lwork = -1;
753 integer rank;
754 double ans;
755 dgelss_ (&m,&n,&nrhs,
756 AC,&m,
757 xp,&ldb,
758 S,
759 &rcond,&rank,
760 &ans,&lwork,
761 &res);
762 lwork = ceil(ans);
763 //printf("ans = %d\n",lwork);
764 double * work = (double*)malloc(lwork*sizeof(double));
765 dgelss_ (&m,&n,&nrhs,
766 AC,&m,
767 xp,&ldb,
768 S,
769 &rcond,&rank,
770 work,&lwork,
771 &res);
772 if(res>0) {
773 return NOCONVER;
774 }
775 CHECK(res,res);
776 free(work);
777 free(S);
778 free(AC);
779 OK
780}
781
782//////////////////// least squares complex linear system using SVD ////////////
783
784// not in clapack.h
785
786int zgelss_(integer *m, integer *n, integer *nhrs,
787 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublereal *s,
788 doublereal *rcond, integer* rank,
789 doublecomplex *work, integer* lwork, doublereal* rwork,
790 integer *info);
791
792int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
793 integer m = ar;
794 integer n = ac;
795 integer nrhs = bc;
796 integer ldb = xr;
797 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
798 DEBUGMSG("linearSolveSVDC_l");
799 doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
800 double*S = (double*)malloc(MIN(m,n)*sizeof(double));
801 double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double));
802 memcpy(AC,ap,m*n*sizeof(doublecomplex));
803 if (m>=n) {
804 memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));
805 } else {
806 int k;
807 for(k = 0; k<nrhs; k++) {
808 memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));
809 }
810 }
811 integer res;
812 integer lwork = -1;
813 integer rank;
814 doublecomplex ans;
815 zgelss_ (&m,&n,&nrhs,
816 AC,&m,
817 xp,&ldb,
818 S,
819 &rcond,&rank,
820 &ans,&lwork,
821 RWORK,
822 &res);
823 lwork = ceil(ans.r);
824 //printf("ans = %d\n",lwork);
825 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
826 zgelss_ (&m,&n,&nrhs,
827 AC,&m,
828 xp,&ldb,
829 S,
830 &rcond,&rank,
831 work,&lwork,
832 RWORK,
833 &res);
834 if(res>0) {
835 return NOCONVER;
836 }
837 CHECK(res,res);
838 free(work);
839 free(RWORK);
840 free(S);
841 free(AC);
842 OK
843}
844
845//////////////////// Cholesky factorization /////////////////////////
846
847/* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a,
848 integer *lda, integer *info);
849
850int chol_l_H(KCMAT(a),CMAT(l)) {
851 integer n = ar;
852 REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
853 DEBUGMSG("chol_l_H");
854 memcpy(lp,ap,n*n*sizeof(doublecomplex));
855 char uplo = 'U';
856 integer res;
857 zpotrf_ (&uplo,&n,lp,&n,&res);
858 CHECK(res>0,NODEFPOS);
859 CHECK(res,res);
860 doublecomplex zero = {0.,0.};
861 int r,c;
862 for (r=0; r<lr-1; r++) {
863 for(c=r+1; c<lc; c++) {
864 lp[r*lc+c] = zero;
865 }
866 }
867 OK
868}
869
870
871/* Subroutine */ int dpotrf_(char *uplo, integer *n, doublereal *a, integer *
872 lda, integer *info);
873
874int chol_l_S(KDMAT(a),DMAT(l)) {
875 integer n = ar;
876 REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
877 DEBUGMSG("chol_l_S");
878 memcpy(lp,ap,n*n*sizeof(double));
879 char uplo = 'U';
880 integer res;
881 dpotrf_ (&uplo,&n,lp,&n,&res);
882 CHECK(res>0,NODEFPOS);
883 CHECK(res,res);
884 int r,c;
885 for (r=0; r<lr-1; r++) {
886 for(c=r+1; c<lc; c++) {
887 lp[r*lc+c] = 0.;
888 }
889 }
890 OK
891}
892
893//////////////////// QR factorization /////////////////////////
894
895/* Subroutine */ int dgeqr2_(integer *m, integer *n, doublereal *a, integer *
896 lda, doublereal *tau, doublereal *work, integer *info);
897
898int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
899 integer m = ar;
900 integer n = ac;
901 integer mn = MIN(m,n);
902 REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE);
903 DEBUGMSG("qr_l_R");
904 double *WORK = (double*)malloc(n*sizeof(double));
905 CHECK(!WORK,MEM);
906 memcpy(rp,ap,m*n*sizeof(double));
907 integer res;
908 dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);
909 CHECK(res,res);
910 free(WORK);
911 OK
912}
913
914/* Subroutine */ int zgeqr2_(integer *m, integer *n, doublecomplex *a,
915 integer *lda, doublecomplex *tau, doublecomplex *work, integer *info);
916
917int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
918 integer m = ar;
919 integer n = ac;
920 integer mn = MIN(m,n);
921 REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE);
922 DEBUGMSG("qr_l_C");
923 doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex));
924 CHECK(!WORK,MEM);
925 memcpy(rp,ap,m*n*sizeof(doublecomplex));
926 integer res;
927 zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);
928 CHECK(res,res);
929 free(WORK);
930 OK
931}
932
933/* Subroutine */ int dorgqr_(integer *m, integer *n, integer *k, doublereal *
934 a, integer *lda, doublereal *tau, doublereal *work, integer *lwork,
935 integer *info);
936
937int c_dorgqr(KDMAT(a), KDVEC(tau), DMAT(r)) {
938 integer m = ar;
939 integer n = MIN(ac,ar);
940 integer k = taun;
941 DEBUGMSG("c_dorgqr");
942 integer lwork = 8*n; // FIXME
943 double *WORK = (double*)malloc(lwork*sizeof(double));
944 CHECK(!WORK,MEM);
945 memcpy(rp,ap,m*k*sizeof(double));
946 integer res;
947 dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res);
948 CHECK(res,res);
949 free(WORK);
950 OK
951}
952
953/* Subroutine */ int zungqr_(integer *m, integer *n, integer *k,
954 doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
955 work, integer *lwork, integer *info);
956
957int c_zungqr(KCMAT(a), KCVEC(tau), CMAT(r)) {
958 integer m = ar;
959 integer n = MIN(ac,ar);
960 integer k = taun;
961 DEBUGMSG("z_ungqr");
962 integer lwork = 8*n; // FIXME
963 doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
964 CHECK(!WORK,MEM);
965 memcpy(rp,ap,m*k*sizeof(doublecomplex));
966 integer res;
967 zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res);
968 CHECK(res,res);
969 free(WORK);
970 OK
971}
972
973
974//////////////////// Hessenberg factorization /////////////////////////
975
976/* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi,
977 doublereal *a, integer *lda, doublereal *tau, doublereal *work,
978 integer *lwork, integer *info);
979
980int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
981 integer m = ar;
982 integer n = ac;
983 integer mn = MIN(m,n);
984 REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE);
985 DEBUGMSG("hess_l_R");
986 integer lwork = 5*n; // fixme
987 double *WORK = (double*)malloc(lwork*sizeof(double));
988 CHECK(!WORK,MEM);
989 memcpy(rp,ap,m*n*sizeof(double));
990 integer res;
991 integer one = 1;
992 dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);
993 CHECK(res,res);
994 free(WORK);
995 OK
996}
997
998
999/* Subroutine */ int zgehrd_(integer *n, integer *ilo, integer *ihi,
1000 doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
1001 work, integer *lwork, integer *info);
1002
1003int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
1004 integer m = ar;
1005 integer n = ac;
1006 integer mn = MIN(m,n);
1007 REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE);
1008 DEBUGMSG("hess_l_C");
1009 integer lwork = 5*n; // fixme
1010 doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
1011 CHECK(!WORK,MEM);
1012 memcpy(rp,ap,m*n*sizeof(doublecomplex));
1013 integer res;
1014 integer one = 1;
1015 zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);
1016 CHECK(res,res);
1017 free(WORK);
1018 OK
1019}
1020
1021//////////////////// Schur factorization /////////////////////////
1022
1023/* Subroutine */ int dgees_(char *jobvs, char *sort, L_fp select, integer *n,
1024 doublereal *a, integer *lda, integer *sdim, doublereal *wr,
1025 doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work,
1026 integer *lwork, logical *bwork, integer *info);
1027
1028int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) {
1029 integer m = ar;
1030 integer n = ac;
1031 REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
1032 DEBUGMSG("schur_l_R");
1033 //int k;
1034 //printf("---------------------------\n");
1035 //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n");
1036 //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n");
1037 //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n");
1038 memcpy(sp,ap,n*n*sizeof(double));
1039 integer lwork = 6*n; // fixme
1040 double *WORK = (double*)malloc(lwork*sizeof(double));
1041 double *WR = (double*)malloc(n*sizeof(double));
1042 double *WI = (double*)malloc(n*sizeof(double));
1043 // WR and WI not really required in this call
1044 logical *BWORK = (logical*)malloc(n*sizeof(logical));
1045 integer res;
1046 integer sdim;
1047 dgees_ ("V","N",NULL,&n,sp,&n,&sdim,WR,WI,up,&n,WORK,&lwork,BWORK,&res);
1048 //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n");
1049 //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n");
1050 //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n");
1051 if(res>0) {
1052 return NOCONVER;
1053 }
1054 CHECK(res,res);
1055 free(WR);
1056 free(WI);
1057 free(BWORK);
1058 free(WORK);
1059 OK
1060}
1061
1062
1063/* Subroutine */ int zgees_(char *jobvs, char *sort, L_fp select, integer *n,
1064 doublecomplex *a, integer *lda, integer *sdim, doublecomplex *w,
1065 doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork,
1066 doublereal *rwork, logical *bwork, integer *info);
1067
1068int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) {
1069 integer m = ar;
1070 integer n = ac;
1071 REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
1072 DEBUGMSG("schur_l_C");
1073 memcpy(sp,ap,n*n*sizeof(doublecomplex));
1074 integer lwork = 6*n; // fixme
1075 doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
1076 doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex));
1077 // W not really required in this call
1078 logical *BWORK = (logical*)malloc(n*sizeof(logical));
1079 double *RWORK = (double*)malloc(n*sizeof(double));
1080 integer res;
1081 integer sdim;
1082 zgees_ ("V","N",NULL,&n,sp,&n,&sdim,W,
1083 up,&n,
1084 WORK,&lwork,RWORK,BWORK,&res);
1085 if(res>0) {
1086 return NOCONVER;
1087 }
1088 CHECK(res,res);
1089 free(W);
1090 free(BWORK);
1091 free(WORK);
1092 OK
1093}
1094
1095//////////////////// LU factorization /////////////////////////
1096
1097/* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer *
1098 lda, integer *ipiv, integer *info);
1099
1100int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r)) {
1101 integer m = ar;
1102 integer n = ac;
1103 integer mn = MIN(m,n);
1104 REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
1105 DEBUGMSG("lu_l_R");
1106 integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
1107 memcpy(rp,ap,m*n*sizeof(double));
1108 integer res;
1109 dgetrf_ (&m,&n,rp,&m,auxipiv,&res);
1110 if(res>0) {
1111 res = 0; // fixme
1112 }
1113 CHECK(res,res);
1114 int k;
1115 for (k=0; k<mn; k++) {
1116 ipivp[k] = auxipiv[k];
1117 }
1118 free(auxipiv);
1119 OK
1120}
1121
1122
1123/* Subroutine */ int zgetrf_(integer *m, integer *n, doublecomplex *a,
1124 integer *lda, integer *ipiv, integer *info);
1125
1126int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {
1127 integer m = ar;
1128 integer n = ac;
1129 integer mn = MIN(m,n);
1130 REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
1131 DEBUGMSG("lu_l_C");
1132 integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
1133 memcpy(rp,ap,m*n*sizeof(doublecomplex));
1134 integer res;
1135 zgetrf_ (&m,&n,rp,&m,auxipiv,&res);
1136 if(res>0) {
1137 res = 0; // fixme
1138 }
1139 CHECK(res,res);
1140 int k;
1141 for (k=0; k<mn; k++) {
1142 ipivp[k] = auxipiv[k];
1143 }
1144 free(auxipiv);
1145 OK
1146}
1147
1148
1149//////////////////// LU substitution /////////////////////////
1150
1151/* Subroutine */ int dgetrs_(char *trans, integer *n, integer *nrhs,
1152 doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *
1153 ldb, integer *info);
1154
1155int luS_l_R(KDMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x)) {
1156 integer m = ar;
1157 integer n = ac;
1158 integer mrhs = br;
1159 integer nrhs = bc;
1160
1161 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1162 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1163 int k;
1164 for (k=0; k<n; k++) {
1165 auxipiv[k] = (integer)ipivp[k];
1166 }
1167 integer res;
1168 memcpy(xp,bp,mrhs*nrhs*sizeof(double));
1169 dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&m,auxipiv,xp,&mrhs,&res);
1170 CHECK(res,res);
1171 free(auxipiv);
1172 OK
1173}
1174
1175
1176/* Subroutine */ int zgetrs_(char *trans, integer *n, integer *nrhs,
1177 doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b,
1178 integer *ldb, integer *info);
1179
1180int luS_l_C(KCMAT(a), KDVEC(ipiv), KCMAT(b), CMAT(x)) {
1181 integer m = ar;
1182 integer n = ac;
1183 integer mrhs = br;
1184 integer nrhs = bc;
1185
1186 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1187 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1188 int k;
1189 for (k=0; k<n; k++) {
1190 auxipiv[k] = (integer)ipivp[k];
1191 }
1192 integer res;
1193 memcpy(xp,bp,mrhs*nrhs*sizeof(doublecomplex));
1194 zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,xp,&mrhs,&res);
1195 CHECK(res,res);
1196 free(auxipiv);
1197 OK
1198}
1199
1200//////////////////// Matrix Product /////////////////////////
1201
1202void dgemm_(char *, char *, integer *, integer *, integer *,
1203 double *, const double *, integer *, const double *,
1204 integer *, double *, double *, integer *);
1205
1206int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) {
1207 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
1208 DEBUGMSG("dgemm_");
1209 CHECKNANR(a,"NaN multR Input\n")
1210 CHECKNANR(b,"NaN multR Input\n")
1211 integer m = ta?ac:ar;
1212 integer n = tb?br:bc;
1213 integer k = ta?ar:ac;
1214 integer lda = ar;
1215 integer ldb = br;
1216 integer ldc = rr;
1217 double alpha = 1;
1218 double beta = 0;
1219 dgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);
1220 CHECKNANR(r,"NaN multR Output\n")
1221 OK
1222}
1223
1224void zgemm_(char *, char *, integer *, integer *, integer *,
1225 doublecomplex *, const doublecomplex *, integer *, const doublecomplex *,
1226 integer *, doublecomplex *, doublecomplex *, integer *);
1227
1228int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) {
1229 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
1230 DEBUGMSG("zgemm_");
1231 CHECKNANC(a,"NaN multC Input\n")
1232 CHECKNANC(b,"NaN multC Input\n")
1233 integer m = ta?ac:ar;
1234 integer n = tb?br:bc;
1235 integer k = ta?ar:ac;
1236 integer lda = ar;
1237 integer ldb = br;
1238 integer ldc = rr;
1239 doublecomplex alpha = {1,0};
1240 doublecomplex beta = {0,0};
1241 zgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,
1242 ap,&lda,
1243 bp,&ldb,&beta,
1244 rp,&ldc);
1245 CHECKNANC(r,"NaN multC Output\n")
1246 OK
1247}
1248
1249void sgemm_(char *, char *, integer *, integer *, integer *,
1250 float *, const float *, integer *, const float *,
1251 integer *, float *, float *, integer *);
1252
1253int multiplyF(int ta, int tb, KFMAT(a),KFMAT(b),FMAT(r)) {
1254 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
1255 DEBUGMSG("sgemm_");
1256 integer m = ta?ac:ar;
1257 integer n = tb?br:bc;
1258 integer k = ta?ar:ac;
1259 integer lda = ar;
1260 integer ldb = br;
1261 integer ldc = rr;
1262 float alpha = 1;
1263 float beta = 0;
1264 sgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);
1265 OK
1266}
1267
1268void cgemm_(char *, char *, integer *, integer *, integer *,
1269 complex *, const complex *, integer *, const complex *,
1270 integer *, complex *, complex *, integer *);
1271
1272int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) {
1273 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
1274 DEBUGMSG("cgemm_");
1275 integer m = ta?ac:ar;
1276 integer n = tb?br:bc;
1277 integer k = ta?ar:ac;
1278 integer lda = ar;
1279 integer ldb = br;
1280 integer ldc = rr;
1281 complex alpha = {1,0};
1282 complex beta = {0,0};
1283 cgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,
1284 ap,&lda,
1285 bp,&ldb,&beta,
1286 rp,&ldc);
1287 OK
1288}
1289
1290//////////////////// transpose /////////////////////////
1291
1292int transF(KFMAT(x),FMAT(t)) {
1293 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1294 DEBUGMSG("transF");
1295 int i,j;
1296 for (i=0; i<tr; i++) {
1297 for (j=0; j<tc; j++) {
1298 tp[i*tc+j] = xp[j*xc+i];
1299 }
1300 }
1301 OK
1302}
1303
1304int transR(KDMAT(x),DMAT(t)) {
1305 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1306 DEBUGMSG("transR");
1307 int i,j;
1308 for (i=0; i<tr; i++) {
1309 for (j=0; j<tc; j++) {
1310 tp[i*tc+j] = xp[j*xc+i];
1311 }
1312 }
1313 OK
1314}
1315
1316int transQ(KQMAT(x),QMAT(t)) {
1317 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1318 DEBUGMSG("transQ");
1319 int i,j;
1320 for (i=0; i<tr; i++) {
1321 for (j=0; j<tc; j++) {
1322 tp[i*tc+j] = xp[j*xc+i];
1323 }
1324 }
1325 OK
1326}
1327
1328int transC(KCMAT(x),CMAT(t)) {
1329 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1330 DEBUGMSG("transC");
1331 int i,j;
1332 for (i=0; i<tr; i++) {
1333 for (j=0; j<tc; j++) {
1334 tp[i*tc+j] = xp[j*xc+i];
1335 }
1336 }
1337 OK
1338}
1339
1340int transP(KPMAT(x), PMAT(t)) {
1341 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1342 REQUIRES(xs==ts,NOCONVER);
1343 DEBUGMSG("transP");
1344 int i,j;
1345 for (i=0; i<tr; i++) {
1346 for (j=0; j<tc; j++) {
1347 memcpy(tp+(i*tc+j)*xs,xp +(j*xc+i)*xs,xs);
1348 }
1349 }
1350 OK
1351}
1352
1353//////////////////// constant /////////////////////////
1354
1355int constantF(float * pval, FVEC(r)) {
1356 DEBUGMSG("constantF")
1357 int k;
1358 double val = *pval;
1359 for(k=0;k<rn;k++) {
1360 rp[k]=val;
1361 }
1362 OK
1363}
1364
1365int constantR(double * pval, DVEC(r)) {
1366 DEBUGMSG("constantR")
1367 int k;
1368 double val = *pval;
1369 for(k=0;k<rn;k++) {
1370 rp[k]=val;
1371 }
1372 OK
1373}
1374
1375int constantQ(complex* pval, QVEC(r)) {
1376 DEBUGMSG("constantQ")
1377 int k;
1378 complex val = *pval;
1379 for(k=0;k<rn;k++) {
1380 rp[k]=val;
1381 }
1382 OK
1383}
1384
1385int constantC(doublecomplex* pval, CVEC(r)) {
1386 DEBUGMSG("constantC")
1387 int k;
1388 doublecomplex val = *pval;
1389 for(k=0;k<rn;k++) {
1390 rp[k]=val;
1391 }
1392 OK
1393}
1394
1395int constantP(void* pval, PVEC(r)) {
1396 DEBUGMSG("constantP")
1397 int k;
1398 for(k=0;k<rn;k++) {
1399 memcpy(rp+k*rs,pval,rs);
1400 }
1401 OK
1402}
1403
1404//////////////////// float-double conversion /////////////////////////
1405
1406int float2double(FVEC(x),DVEC(y)) {
1407 DEBUGMSG("float2double")
1408 int k;
1409 for(k=0;k<xn;k++) {
1410 yp[k]=xp[k];
1411 }
1412 OK
1413}
1414
1415int double2float(DVEC(x),FVEC(y)) {
1416 DEBUGMSG("double2float")
1417 int k;
1418 for(k=0;k<xn;k++) {
1419 yp[k]=xp[k];
1420 }
1421 OK
1422}
1423
1424//////////////////// conjugate /////////////////////////
1425
1426int conjugateQ(KQVEC(x),QVEC(t)) {
1427 REQUIRES(xn==tn,BAD_SIZE);
1428 DEBUGMSG("conjugateQ");
1429 int k;
1430 for(k=0;k<xn;k++) {
1431 tp[k].r = xp[k].r;
1432 tp[k].i = -xp[k].i;
1433 }
1434 OK
1435}
1436
1437int conjugateC(KCVEC(x),CVEC(t)) {
1438 REQUIRES(xn==tn,BAD_SIZE);
1439 DEBUGMSG("conjugateC");
1440 int k;
1441 for(k=0;k<xn;k++) {
1442 tp[k].r = xp[k].r;
1443 tp[k].i = -xp[k].i;
1444 }
1445 OK
1446}
1447
1448//////////////////// step /////////////////////////
1449
1450int stepF(FVEC(x),FVEC(y)) {
1451 DEBUGMSG("stepF")
1452 int k;
1453 for(k=0;k<xn;k++) {
1454 yp[k]=xp[k]>0;
1455 }
1456 OK
1457}
1458
1459int stepD(DVEC(x),DVEC(y)) {
1460 DEBUGMSG("stepD")
1461 int k;
1462 for(k=0;k<xn;k++) {
1463 yp[k]=xp[k]>0;
1464 }
1465 OK
1466}
1467
1468//////////////////// cond /////////////////////////
1469
1470int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) {
1471 REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);
1472 DEBUGMSG("condF")
1473 int k;
1474 for(k=0;k<xn;k++) {
1475 rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]);
1476 }
1477 OK
1478}
1479
1480int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {
1481 REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);
1482 DEBUGMSG("condD")
1483 int k;
1484 for(k=0;k<xn;k++) {
1485 rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]);
1486 }
1487 OK
1488}
1489
diff --git a/packages/base/src/C/lapack-aux.h b/packages/base/src/C/lapack-aux.h
new file mode 100644
index 0000000..a3f1899
--- /dev/null
+++ b/packages/base/src/C/lapack-aux.h
@@ -0,0 +1,60 @@
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 FVEC(A) int A##n, float*A##p
40#define DVEC(A) int A##n, double*A##p
41#define QVEC(A) int A##n, complex*A##p
42#define CVEC(A) int A##n, doublecomplex*A##p
43#define PVEC(A) int A##n, void* A##p, int A##s
44#define FMAT(A) int A##r, int A##c, float* A##p
45#define DMAT(A) int A##r, int A##c, double* A##p
46#define QMAT(A) int A##r, int A##c, complex* A##p
47#define CMAT(A) int A##r, int A##c, doublecomplex* A##p
48#define PMAT(A) int A##r, int A##c, void* A##p, int A##s
49
50#define KFVEC(A) int A##n, const float*A##p
51#define KDVEC(A) int A##n, const double*A##p
52#define KQVEC(A) int A##n, const complex*A##p
53#define KCVEC(A) int A##n, const doublecomplex*A##p
54#define KPVEC(A) int A##n, const void* A##p, int A##s
55#define KFMAT(A) int A##r, int A##c, const float* A##p
56#define KDMAT(A) int A##r, int A##c, const double* A##p
57#define KQMAT(A) int A##r, int A##c, const complex* A##p
58#define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p
59#define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s
60
diff --git a/packages/hmatrix/src/Data/Packed.hs b/packages/base/src/Data/Packed.hs
index 957aab8..c66718a 100644
--- a/packages/hmatrix/src/Data/Packed.hs
+++ b/packages/base/src/Data/Packed.hs
@@ -1,29 +1,25 @@
1----------------------------------------------------------------------------- 1-----------------------------------------------------------------------------
2{- | 2{- |
3Module : Data.Packed 3Module : Data.Packed
4Copyright : (c) Alberto Ruiz 2006-2010 4Copyright : (c) Alberto Ruiz 2006-2014
5License : GPL-style 5License : BSD3
6 6Maintainer : Alberto Ruiz
7Maintainer : Alberto Ruiz (aruiz at um dot es)
8Stability : provisional 7Stability : provisional
9Portability : uses ffi
10 8
11Types for dense 'Vector' and 'Matrix' of 'Storable' elements. 9Types for dense 'Vector' and 'Matrix' of 'Storable' elements.
12 10
13-} 11-}
14----------------------------------------------------------------------------- 12-----------------------------------------------------------------------------
15{-# OPTIONS_HADDOCK hide #-}
16 13
17module Data.Packed ( 14module Data.Packed (
15 -- * Vector
16 --
17 -- | Vectors are @Data.Vector.Storable.Vector@ from the \"vector\" package.
18 module Data.Packed.Vector, 18 module Data.Packed.Vector,
19 -- * Matrix
19 module Data.Packed.Matrix, 20 module Data.Packed.Matrix,
20-- module Numeric.Conversion,
21-- module Data.Packed.Random,
22-- module Data.Complex
23) where 21) where
24 22
25import Data.Packed.Vector 23import Data.Packed.Vector
26import Data.Packed.Matrix 24import Data.Packed.Matrix
27--import Data.Packed.Random 25
28--import Data.Complex
29--import Numeric.Conversion
diff --git a/packages/hmatrix/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs
index 471e560..777b6c5 100644
--- a/packages/hmatrix/src/Data/Packed/Development.hs
+++ b/packages/base/src/Data/Packed/Development.hs
@@ -3,9 +3,8 @@
3-- | 3-- |
4-- Module : Data.Packed.Development 4-- Module : Data.Packed.Development
5-- Copyright : (c) Alberto Ruiz 2009 5-- Copyright : (c) Alberto Ruiz 2009
6-- License : GPL 6-- License : BSD3
7-- 7-- Maintainer : Alberto Ruiz
8-- Maintainer : Alberto Ruiz <aruiz@um.es>
9-- Stability : provisional 8-- Stability : provisional
10-- Portability : portable 9-- Portability : portable
11-- 10--
@@ -14,7 +13,6 @@
14-- in the @examples\/devel@ folder included in the package. 13-- in the @examples\/devel@ folder included in the package.
15-- 14--
16----------------------------------------------------------------------------- 15-----------------------------------------------------------------------------
17{-# OPTIONS_HADDOCK hide #-}
18 16
19module Data.Packed.Development ( 17module Data.Packed.Development (
20 createVector, createMatrix, 18 createVector, createMatrix,
@@ -30,3 +28,4 @@ module Data.Packed.Development (
30) where 28) where
31 29
32import Data.Packed.Internal 30import Data.Packed.Internal
31
diff --git a/packages/hmatrix/src/Data/Packed/Foreign.hs b/packages/base/src/Data/Packed/Foreign.hs
index 1ec3694..efa51ca 100644
--- a/packages/hmatrix/src/Data/Packed/Foreign.hs
+++ b/packages/base/src/Data/Packed/Foreign.hs
@@ -6,7 +6,6 @@
6-- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) 6-- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3)
7-- @ 7-- @
8-- 8--
9{-# OPTIONS_HADDOCK hide #-}
10module Data.Packed.Foreign 9module Data.Packed.Foreign
11 ( app 10 ( app
12 , appVector, appVectorLen 11 , appVector, appVectorLen
diff --git a/packages/hmatrix/src/Data/Packed/Internal.hs b/packages/base/src/Data/Packed/Internal.hs
index 537e51e..537e51e 100644
--- a/packages/hmatrix/src/Data/Packed/Internal.hs
+++ b/packages/base/src/Data/Packed/Internal.hs
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Common.hs b/packages/base/src/Data/Packed/Internal/Common.hs
index edef3c2..615bbdf 100644
--- a/packages/hmatrix/src/Data/Packed/Internal/Common.hs
+++ b/packages/base/src/Data/Packed/Internal/Common.hs
@@ -1,18 +1,15 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2-----------------------------------------------------------------------------
3-- | 2-- |
4-- Module : Data.Packed.Internal.Common 3-- Module : Data.Packed.Internal.Common
5-- Copyright : (c) Alberto Ruiz 2007 4-- Copyright : (c) Alberto Ruiz 2007
6-- License : GPL-style 5-- License : BSD3
7-- 6-- Maintainer : Alberto Ruiz
8-- Maintainer : Alberto Ruiz <aruiz@um.es>
9-- Stability : provisional 7-- Stability : provisional
10-- Portability : portable (uses FFI) 8--
11-- 9--
12-- Development utilities. 10-- Development utilities.
13-- 11--
14----------------------------------------------------------------------------- 12
15-- #hide
16 13
17module Data.Packed.Internal.Common( 14module Data.Packed.Internal.Common(
18 Adapt, 15 Adapt,
@@ -21,12 +18,11 @@ module Data.Packed.Internal.Common(
21 (//), check, mbCatch, 18 (//), check, mbCatch,
22 splitEvery, common, compatdim, 19 splitEvery, common, compatdim,
23 fi, 20 fi,
24 table 21 table,
22 finit
25) where 23) where
26 24
27import Foreign
28import Control.Monad(when) 25import Control.Monad(when)
29import Foreign.C.String(peekCString)
30import Foreign.C.Types 26import Foreign.C.Types
31import Foreign.Storable.Complex() 27import Foreign.Storable.Complex()
32import Data.List(transpose,intersperse) 28import Data.List(transpose,intersperse)
@@ -153,19 +149,12 @@ check msg f = do
153 finit 149 finit
154#endif 150#endif
155 err <- f 151 err <- f
156 when (err/=0) $ if err > 1024 152 when (err/=0) $ error (msg++": "++errorCode err)
157 then (error (msg++": "++errorCode err)) -- our errors
158 else do -- GSL errors
159 ps <- gsl_strerror err
160 s <- peekCString ps
161 error (msg++": "++s)
162 return () 153 return ()
163 154
164-- | description of GSL error codes
165foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
166
167-- | Error capture and conversion to Maybe 155-- | Error capture and conversion to Maybe
168mbCatch :: IO x -> IO (Maybe x) 156mbCatch :: IO x -> IO (Maybe x)
169mbCatch act = E.catch (Just `fmap` act) f 157mbCatch act = E.catch (Just `fmap` act) f
170 where f :: SomeException -> IO (Maybe x) 158 where f :: SomeException -> IO (Maybe x)
171 f _ = return Nothing 159 f _ = return Nothing
160
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs b/packages/base/src/Data/Packed/Internal/Matrix.hs
index 9719fc0..9b831cc 100644
--- a/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs
+++ b/packages/base/src/Data/Packed/Internal/Matrix.hs
@@ -2,20 +2,16 @@
2{-# LANGUAGE FlexibleContexts #-} 2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-} 3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE BangPatterns #-} 4{-# LANGUAGE BangPatterns #-}
5----------------------------------------------------------------------------- 5
6-- | 6-- |
7-- Module : Data.Packed.Internal.Matrix 7-- Module : Data.Packed.Internal.Matrix
8-- Copyright : (c) Alberto Ruiz 2007 8-- Copyright : (c) Alberto Ruiz 2007
9-- License : GPL-style 9-- License : BSD3
10-- 10-- Maintainer : Alberto Ruiz
11-- Maintainer : Alberto Ruiz <aruiz@um.es>
12-- Stability : provisional 11-- Stability : provisional
13-- Portability : portable (uses FFI)
14-- 12--
15-- Internal matrix representation 13-- Internal matrix representation
16-- 14--
17-----------------------------------------------------------------------------
18-- #hide
19 15
20module Data.Packed.Internal.Matrix( 16module Data.Packed.Internal.Matrix(
21 Matrix(..), rows, cols, cdat, fdat, 17 Matrix(..), rows, cols, cdat, fdat,
@@ -30,7 +26,6 @@ module Data.Packed.Internal.Matrix(
30 subMatrix, 26 subMatrix,
31 liftMatrix, liftMatrix2, 27 liftMatrix, liftMatrix2,
32 (@@>), atM', 28 (@@>), atM',
33 saveMatrix,
34 singleton, 29 singleton,
35 emptyM, 30 emptyM,
36 size, shSize, conformVs, conformMs, conformVTo, conformMTo 31 size, shSize, conformVs, conformMs, conformVTo, conformMTo
@@ -46,7 +41,6 @@ import Foreign.Ptr(Ptr, castPtr)
46import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf) 41import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf)
47import Data.Complex(Complex) 42import Data.Complex(Complex)
48import Foreign.C.Types 43import Foreign.C.Types
49import Foreign.C.String(newCString)
50import System.IO.Unsafe(unsafePerformIO) 44import System.IO.Unsafe(unsafePerformIO)
51import Control.DeepSeq 45import Control.DeepSeq
52 46
@@ -290,34 +284,6 @@ instance Element (Complex Double) where
290 284
291------------------------------------------------------------------- 285-------------------------------------------------------------------
292 286
293transdata' :: Storable a => Int -> Vector a -> Int -> Vector a
294transdata' c1 v c2 =
295 if noneed
296 then v
297 else unsafePerformIO $ do
298 w <- createVector (r2*c2)
299 unsafeWith v $ \p ->
300 unsafeWith w $ \q -> do
301 let go (-1) _ = return ()
302 go !i (-1) = go (i-1) (c1-1)
303 go !i !j = do x <- peekElemOff p (i*c1+j)
304 pokeElemOff q (j*c2+i) x
305 go i (j-1)
306 go (r1-1) (c1-1)
307 return w
308 where r1 = dim v `div` c1
309 r2 = dim v `div` c2
310 noneed = dim v == 0 || r1 == 1 || c1 == 1
311
312-- {-# SPECIALIZE transdata' :: Int -> Vector Double -> Int -> Vector Double #-}
313-- {-# SPECIALIZE transdata' :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double) #-}
314
315-- I don't know how to specialize...
316-- The above pragmas only seem to work on top level defs
317-- Fortunately everything seems to work using the above class
318
319-- C versions, still a little faster:
320
321transdataAux fun c1 d c2 = 287transdataAux fun c1 d c2 =
322 if noneed 288 if noneed
323 then d 289 then d
@@ -354,16 +320,6 @@ foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -
354 320
355---------------------------------------------------------------------- 321----------------------------------------------------------------------
356 322
357constant' v n = unsafePerformIO $ do
358 w <- createVector n
359 unsafeWith w $ \p -> do
360 let go (-1) = return ()
361 go !k = pokeElemOff p k v >> go (k-1)
362 go (n-1)
363 return w
364
365-- C versions
366
367constantAux fun x n = unsafePerformIO $ do 323constantAux fun x n = unsafePerformIO $ do
368 v <- createVector n 324 v <- createVector n
369 px <- newArray [x] 325 px <- newArray [x]
@@ -371,20 +327,12 @@ constantAux fun x n = unsafePerformIO $ do
371 free px 327 free px
372 return v 328 return v
373 329
374constantF :: Float -> Int -> Vector Float
375constantF = constantAux cconstantF
376foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF 330foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF
377 331
378constantR :: Double -> Int -> Vector Double
379constantR = constantAux cconstantR
380foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV 332foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV
381 333
382constantQ :: Complex Float -> Int -> Vector (Complex Float)
383constantQ = constantAux cconstantQ
384foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV 334foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
385 335
386constantC :: Complex Double -> Int -> Vector (Complex Double)
387constantC = constantAux cconstantC
388foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV 336foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV
389 337
390constantP :: Storable a => a -> Int -> Vector a 338constantP :: Storable a => a -> Int -> Vector a
@@ -429,23 +377,6 @@ subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m)
429 377
430-------------------------------------------------------------------------- 378--------------------------------------------------------------------------
431 379
432-- | Saves a matrix as 2D ASCII table.
433saveMatrix :: FilePath
434 -> String -- ^ format (%f, %g, %e)
435 -> Matrix Double
436 -> IO ()
437saveMatrix filename fmt m = do
438 charname <- newCString filename
439 charfmt <- newCString fmt
440 let o = if orderOf m == RowMajor then 1 else 0
441 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
442 free charname
443 free charfmt
444
445foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
446
447----------------------------------------------------------------------
448
449maxZ xs = if minimum xs == 0 then 0 else maximum xs 380maxZ xs = if minimum xs == 0 then 0 else maximum xs
450 381
451conformMs ms = map (conformMTo (r,c)) ms 382conformMs ms = map (conformMTo (r,c)) ms
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs b/packages/base/src/Data/Packed/Internal/Signatures.hs
index 2835720..acc3070 100644
--- a/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs
+++ b/packages/base/src/Data/Packed/Internal/Signatures.hs
@@ -1,16 +1,13 @@
1-----------------------------------------------------------------------------
2-- | 1-- |
3-- Module : Data.Packed.Internal.Signatures 2-- Module : Data.Packed.Internal.Signatures
4-- Copyright : (c) Alberto Ruiz 2009 3-- Copyright : (c) Alberto Ruiz 2009
5-- License : GPL-style 4-- License : BSD3
6-- 5-- Maintainer : Alberto Ruiz
7-- Maintainer : Alberto Ruiz <aruiz@um.es>
8-- Stability : provisional 6-- Stability : provisional
9-- Portability : portable (uses FFI)
10-- 7--
11-- Signatures of the C functions. 8-- Signatures of the C functions.
12-- 9--
13----------------------------------------------------------------------------- 10
14 11
15module Data.Packed.Internal.Signatures where 12module Data.Packed.Internal.Signatures where
16 13
@@ -70,3 +67,4 @@ type TVCV = CInt -> PD -> TCV --
70type TCVM = CInt -> PC -> TM -- 67type TCVM = CInt -> PC -> TM --
71type TMCVM = CInt -> CInt -> PD -> TCVM -- 68type TMCVM = CInt -> CInt -> PD -> TCVM --
72type TMMCVM = CInt -> CInt -> PD -> TMCVM -- 69type TMMCVM = CInt -> CInt -> PD -> TMCVM --
70
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Vector.hs b/packages/base/src/Data/Packed/Internal/Vector.hs
index 6d03438..d0bc143 100644
--- a/packages/hmatrix/src/Data/Packed/Internal/Vector.hs
+++ b/packages/base/src/Data/Packed/Internal/Vector.hs
@@ -1,17 +1,14 @@
1{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-} 1{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}
2-----------------------------------------------------------------------------
3-- | 2-- |
4-- Module : Data.Packed.Internal.Vector 3-- Module : Data.Packed.Internal.Vector
5-- Copyright : (c) Alberto Ruiz 2007 4-- Copyright : (c) Alberto Ruiz 2007
6-- License : GPL-style 5-- License : BSD3
7-- 6-- Maintainer : Alberto Ruiz
8-- Maintainer : Alberto Ruiz <aruiz@um.es>
9-- Stability : provisional 7-- Stability : provisional
10-- Portability : portable (uses FFI)
11-- 8--
12-- Vector implementation 9-- Vector implementation
13-- 10--
14----------------------------------------------------------------------------- 11--------------------------------------------------------------------------------
15 12
16module Data.Packed.Internal.Vector ( 13module Data.Packed.Internal.Vector (
17 Vector, dim, 14 Vector, dim,
@@ -24,7 +21,6 @@ module Data.Packed.Internal.Vector (
24 asComplex, asReal, float2DoubleV, double2FloatV, 21 asComplex, asReal, float2DoubleV, double2FloatV,
25 stepF, stepD, condF, condD, 22 stepF, stepD, condF, condD,
26 conjugateQ, conjugateC, 23 conjugateQ, conjugateC,
27 fwriteVector, freadVector, fprintfVector, fscanfVector,
28 cloneVector, 24 cloneVector,
29 unsafeToForeignPtr, 25 unsafeToForeignPtr,
30 unsafeFromForeignPtr, 26 unsafeFromForeignPtr,
@@ -33,12 +29,10 @@ module Data.Packed.Internal.Vector (
33 29
34import Data.Packed.Internal.Common 30import Data.Packed.Internal.Common
35import Data.Packed.Internal.Signatures 31import Data.Packed.Internal.Signatures
36import Foreign.Marshal.Alloc(free)
37import Foreign.Marshal.Array(peekArray, copyArray, advancePtr) 32import Foreign.Marshal.Array(peekArray, copyArray, advancePtr)
38import Foreign.ForeignPtr(ForeignPtr, castForeignPtr) 33import Foreign.ForeignPtr(ForeignPtr, castForeignPtr)
39import Foreign.Ptr(Ptr) 34import Foreign.Ptr(Ptr)
40import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf) 35import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf)
41import Foreign.C.String
42import Foreign.C.Types 36import Foreign.C.Types
43import Data.Complex 37import Data.Complex
44import Control.Monad(when) 38import Control.Monad(when)
@@ -474,48 +468,4 @@ mapVectorWithIndex f v = unsafePerformIO $ do
474 return w 468 return w
475{-# INLINE mapVectorWithIndex #-} 469{-# INLINE mapVectorWithIndex #-}
476 470
477-------------------------------------------------------------------
478
479
480-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
481fscanfVector :: FilePath -> Int -> IO (Vector Double)
482fscanfVector filename n = do
483 charname <- newCString filename
484 res <- createVector n
485 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
486 free charname
487 return res
488
489foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
490
491-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
492fprintfVector :: FilePath -> String -> Vector Double -> IO ()
493fprintfVector filename fmt v = do
494 charname <- newCString filename
495 charfmt <- newCString fmt
496 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
497 free charname
498 free charfmt
499
500foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
501
502-- | Loads a vector from a binary file (the number of elements must be known in advance).
503freadVector :: FilePath -> Int -> IO (Vector Double)
504freadVector filename n = do
505 charname <- newCString filename
506 res <- createVector n
507 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
508 free charname
509 return res
510
511foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
512
513-- | Saves the elements of a vector to a binary file.
514fwriteVector :: FilePath -> Vector Double -> IO ()
515fwriteVector filename v = do
516 charname <- newCString filename
517 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
518 free charname
519
520foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
521 471
diff --git a/packages/hmatrix/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs
index d94d167..d94d167 100644
--- a/packages/hmatrix/src/Data/Packed/Matrix.hs
+++ b/packages/base/src/Data/Packed/Matrix.hs
diff --git a/packages/hmatrix/src/Data/Packed/ST.hs b/packages/base/src/Data/Packed/ST.hs
index 1cef296..dae457c 100644
--- a/packages/hmatrix/src/Data/Packed/ST.hs
+++ b/packages/base/src/Data/Packed/ST.hs
@@ -6,9 +6,8 @@
6-- | 6-- |
7-- Module : Data.Packed.ST 7-- Module : Data.Packed.ST
8-- Copyright : (c) Alberto Ruiz 2008 8-- Copyright : (c) Alberto Ruiz 2008
9-- License : GPL-style 9-- License : BSD3
10-- 10-- Maintainer : Alberto Ruiz
11-- Maintainer : Alberto Ruiz <aruiz@um.es>
12-- Stability : provisional 11-- Stability : provisional
13-- Portability : portable 12-- Portability : portable
14-- 13--
@@ -16,7 +15,6 @@
16-- See examples/inplace.hs in the distribution. 15-- See examples/inplace.hs in the distribution.
17-- 16--
18----------------------------------------------------------------------------- 17-----------------------------------------------------------------------------
19{-# OPTIONS_HADDOCK hide #-}
20 18
21module Data.Packed.ST ( 19module Data.Packed.ST (
22 -- * Mutable Vectors 20 -- * Mutable Vectors
@@ -177,3 +175,4 @@ newUndefinedMatrix ord r c = unsafeIOToST $ fmap STMatrix $ createMatrix ord r c
177{-# NOINLINE newMatrix #-} 175{-# NOINLINE newMatrix #-}
178newMatrix :: Storable t => t -> Int -> Int -> ST s (STMatrix s t) 176newMatrix :: Storable t => t -> Int -> Int -> ST s (STMatrix s t)
179newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c) 177newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c)
178
diff --git a/packages/hmatrix/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs
index b5a4318..b5a4318 100644
--- a/packages/hmatrix/src/Data/Packed/Vector.hs
+++ b/packages/base/src/Data/Packed/Vector.hs
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 09:59:06 +0000