Changeset 28194 in project
 Timestamp:
 01/30/13 08:16:33 (8 years ago)
 Location:
 release/4
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

release/4/matrixutils/trunk/matrixutils.scm
r26824 r28194 4 4 ;; 5 5 ;; 6 ;; Copyright 2007201 2Ivan Raikov and the Okinawa Institute of6 ;; Copyright 20072013 Ivan Raikov and the Okinawa Institute of 7 7 ;; Science and Technology. 8 8 ;; … … 38 38 39 39 defineutilitymatrices 40 withutilitymatrices) 41 42 (import scheme chicken datastructures srfi4) 40 withutilitymatrices 41 42 f64vectorrepmat f64vectormeshgrid 43 ) 44 45 (import scheme chicken foreign datastructures srfi4) 43 46 (requireextension srfi4 srfi42 srfi4comprehensions blas) 44 47 … … 428 431 ) 429 432 433 434 #> 435 436 /* repeat a block of memory rep times */ 437 void cmemrep(double *dest, size_t chunk, int rep) 438 { 439 if(rep == 1) return; 440 memcpy(dest + chunk, dest, sizeof(double)*chunk); 441 if(rep & 1) { 442 dest += chunk; 443 memcpy(dest + chunk, dest, sizeof(double)*chunk); 444 } 445 /* now repeat using a block twice as big */ 446 cmemrep(dest, chunk<<1, rep>>1); 447 } 448 449 void crepmat(double *dest, const double *src, int ndim, 450 int *destdimsize, int *dimsize, const int *dims, int *rep) 451 { 452 int i, chunk; 453 int d = ndim1; 454 /* copy the first repetition into dest */ 455 if(d == 0) 456 { 457 chunk = dimsize[0]; 458 memcpy(dest,src,sizeof(double)*chunk); 459 } 460 else 461 { 462 /* recursively repeat each slice of src */ 463 for(i=0; i<dims[d]; i++) { 464 crepmat(dest + i*destdimsize[d1], src + i*dimsize[d1], 465 ndim1, destdimsize, dimsize, dims, rep); 466 } 467 chunk = destdimsize[d1]*dims[d]; 468 } 469 /* copy the result rep1 times */ 470 cmemrep(dest,chunk,rep[d]); 471 } 472 473 <# 474 475 (define crepmat (foreignsafelambda void "crepmat" f64vector f64vector int s32vector s32vector s32vector s32vector)) 476 477 (define cmemrep (foreignsafelambda void "cmemrep" f64vector int int)) 478 479 (define (f64vectorrepmat src dims reps) 480 (let* ((ndim (length dims)) 481 (dimsize (makes32vector ndim)) 482 (vdims (list>s32vector dims))) 483 484 (assert (every positive? reps)) 485 486 ;; compute dimension sizes 487 (s32vectorset! dimsize 0 (car dims)) 488 (let recur ((i 1) (dims (cdr dims))) 489 (if (< i ndim) 490 (begin 491 (s32vectorset! dimsize i (* (car dims) (s32vectorref dimsize ( i 1)))) 492 (recur (+ 1 i) (cdr dims))) 493 )) 494 495 ;; determine repetition vector 496 (let* ((nrep (length reps)) 497 (ndimdest (if (> nrep ndim) nrep ndim)) 498 (rep (makes32vector ndimdest))) 499 (let recur ((i 0) (reps reps)) 500 (if (< i nrep) 501 (let ((repv (car reps))) 502 (s32vectorset! rep i repv) 503 (recur (+ 1 i) (cdr reps))))) 504 505 ;; compute output size 506 (let ((destdims (makes32vector ndimdest 0))) 507 (let recur ((i 0)) 508 (if (< i ndim) 509 (begin (s32vectorset! destdims i 510 (* (s32vectorref vdims i) 511 (s32vectorref rep i ))) 512 (recur (+ 1 i))))) 513 514 (let ((extrarep 515 (let recur ((i ndim) (extrarep 1)) 516 (if (< i ndimdest) 517 (let ((ri (s32vectorref rep i))) 518 (s32vectorset! destdims i ri) 519 (recur (+ 1 i) (* extrarep ri))) 520 extrarep 521 )))) 522 523 (let ((destdimsize (makes32vector ndimdest))) 524 (s32vectorset! destdimsize 0 (s32vectorref destdims 0)) 525 (let recur ((i 1)) 526 (if (< i ndimdest) 527 (begin 528 (s32vectorset! destdimsize i 529 (* (s32vectorref destdimsize ( i 1)) 530 (s32vectorref destdims i))) 531 (recur (+ 1 i))))) 532 533 ;; return array 534 (let* ((destlen (s32vectorref destdimsize ( ndimdest 1))) 535 (dest (makef64vector destlen))) 536 537 (crepmat dest src ndim destdimsize dimsize vdims rep) 538 (if (> ndimdest ndim) 539 (let ((n (s32vectorref destdimsize ( ndim 1)))) 540 (cmemrep dest n extrarep))) 541 dest 542 )) 543 )) 544 ))) 545 546 547 548 549 (define (f64vectormeshgrid x y z) 550 (let ((lenx (f64vectorlength x)) 551 (leny (f64vectorlength y)) 552 (lenz (f64vectorlength z))) 553 (let ((dimx (list 1 lenx)) 554 (dimy (list 1 leny)) 555 (dimz (list 1 lenz))) 556 (let ((xx (f64vectorrepmat 557 (f64vectorrepmat x dimx (list leny 1)) 558 (list leny lenx) 559 (list 1 1 lenz))) 560 (yy (f64vectorrepmat 561 (f64vectorrepmat y dimy (list 1 lenx)) 562 (list lenx leny) 563 (list 1 1 lenz))) 564 (zz (f64vectorrepmat z dimz (list (* lenx leny) 1)))) 565 (list xx yy zz)) 566 ))) 567 568 430 569 ) 431 432 433 434 570 571 572 
release/4/srfi4utils/trunk/srfi4utils.meta
r23287 r28194 16 16 (category data) 17 17 18 ; A list of eggs matrixutils depends on.18 ; A list of eggs srfi4utils depends on. 19 19 20 20 (needs srfi42 srfi4comprehensions ) 21 22 (testdepends test) 21 23 22 24 (author "Ivan Raikov") 
release/4/srfi4utils/trunk/srfi4utils.scm
r28189 r28194 46 46 s16vectormergesort! u16vectormergesort! s8vectormergesort! u8vectormergesort! 47 47 48 f64vectorrepmat f64vectormeshgrid49 48 ) 50 49 51 (import scheme chicken foreign)50 (import scheme chicken) 52 51 (requirelibrary srfi1 extras) 53 52 (import (only srfi1 fold every) … … 157 156 158 157 (define (srfi4vectorquicksort! vectorref vectorset! vectorlength) 159 (lambda ( v elt<. rest)158 (lambda (elt< v . rest) 160 159 (letoptionals rest ((start 0) (end (vectorlength v))) 161 160 (let recur ((l start) (r end)) ; Sort the range [l,r). … … 331 330 (define u8vectormergesort! 332 331 (srfi4vectormergesort! u8vectorlength u8vectorref u8vectorset! u8vectormerge! u8vectorblit!)) 333 334 335 336 337 338 #>339 340 /* repeat a block of memory rep times */341 void cmemrep(double *dest, size_t chunk, int rep)342 {343 if(rep == 1) return;344 memcpy(dest + chunk, dest, sizeof(double)*chunk);345 if(rep & 1) {346 dest += chunk;347 memcpy(dest + chunk, dest, sizeof(double)*chunk);348 }349 /* now repeat using a block twice as big */350 cmemrep(dest, chunk<<1, rep>>1);351 }352 353 void crepmat(double *dest, const double *src, int ndim,354 int *destdimsize, int *dimsize, const int *dims, int *rep)355 {356 int i, chunk;357 int d = ndim1;358 /* copy the first repetition into dest */359 if(d == 0)360 {361 chunk = dimsize[0];362 memcpy(dest,src,sizeof(double)*chunk);363 }364 else365 {366 /* recursively repeat each slice of src */367 for(i=0; i<dims[d]; i++) {368 crepmat(dest + i*destdimsize[d1], src + i*dimsize[d1],369 ndim1, destdimsize, dimsize, dims, rep);370 }371 chunk = destdimsize[d1]*dims[d];372 }373 /* copy the result rep1 times */374 cmemrep(dest,chunk,rep[d]);375 }376 377 <#378 379 (define crepmat (foreignsafelambda void "crepmat" f64vector f64vector int s32vector s32vector s32vector s32vector))380 381 (define cmemrep (foreignsafelambda void "cmemrep" f64vector int int))382 383 (define (f64vectorrepmat src dims reps)384 (let* ((ndim (length dims))385 (dimsize (makes32vector ndim))386 (vdims (list>s32vector dims)))387 388 (assert (every positive? reps))389 390 ;; compute dimension sizes391 (s32vectorset! dimsize 0 (car dims))392 (let recur ((i 1) (dims (cdr dims)))393 (if (< i ndim)394 (begin395 (s32vectorset! dimsize i (* (car dims) (s32vectorref dimsize ( i 1))))396 (recur (+ 1 i) (cdr dims)))397 ))398 399 ;; determine repetition vector400 (let* ((nrep (length reps))401 (ndimdest (if (> nrep ndim) nrep ndim))402 (rep (makes32vector ndimdest)))403 (let recur ((i 0) (reps reps))404 (if (< i nrep)405 (let ((repv (car reps)))406 (s32vectorset! rep i repv)407 (recur (+ 1 i) (cdr reps)))))408 409 ;; compute output size410 (let ((destdims (makes32vector ndimdest 0)))411 (let recur ((i 0))412 (if (< i ndim)413 (begin (s32vectorset! destdims i414 (* (s32vectorref vdims i)415 (s32vectorref rep i )))416 (recur (+ 1 i)))))417 418 (let ((extrarep419 (let recur ((i ndim) (extrarep 1))420 (if (< i ndimdest)421 (let ((ri (s32vectorref rep i)))422 (s32vectorset! destdims i ri)423 (recur (+ 1 i) (* extrarep ri)))424 extrarep425 ))))426 427 (let ((destdimsize (makes32vector ndimdest)))428 (s32vectorset! destdimsize 0 (s32vectorref destdims 0))429 (let recur ((i 1))430 (if (< i ndimdest)431 (begin432 (s32vectorset! destdimsize i433 (* (s32vectorref destdimsize ( i 1))434 (s32vectorref destdims i)))435 (recur (+ 1 i)))))436 437 ;; return array438 (let* ((destlen (s32vectorref destdimsize ( ndimdest 1)))439 (dest (makef64vector destlen)))440 441 (crepmat dest src ndim destdimsize dimsize vdims rep)442 (if (> ndimdest ndim)443 (let ((n (s32vectorref destdimsize ( ndim 1))))444 (cmemrep dest n extrarep)))445 dest446 ))447 ))448 )))449 450 451 452 453 (define (f64vectormeshgrid x y z)454 (let ((lenx (f64vectorlength x))455 (leny (f64vectorlength y))456 (lenz (f64vectorlength z)))457 (let ((dimx (list 1 lenx))458 (dimy (list 1 leny))459 (dimz (list 1 lenz)))460 (let ((xx (f64vectorrepmat461 (f64vectorrepmat x dimx (list leny 1))462 (list leny lenx)463 (list 1 1 lenz)))464 (yy (f64vectorrepmat465 (f64vectorrepmat y dimy (list 1 lenx))466 (list lenx leny)467 (list 1 1 lenz)))468 (zz (f64vectorrepmat z dimz (list (* lenx leny) 1))))469 (list xx yy zz))470 )))471 472 473 332 474 333 ) 
release/4/srfi4utils/trunk/tests/run.scm
r28191 r28194 4 4 (define elt< (lambda (i ai j bj) (< ai bj))) 5 5 (define elt> (lambda (i ai j bj) (> ai bj))) 6 7 (print (f64vectormerge! elt< (f64vector 1 3 5 2 4) 0 3 5 (makef64vector 5 0.) 0))8 9 (print (f64vectormerge! elt< (f64vector 1 3 5 2 4 6) 0 3 6 (makef64vector 6 0.) 0))10 6 11 7 (testgroup
Note: See TracChangeset
for help on using the changeset viewer.