| 9 |
|
#include <stdlib.h> |
| 10 |
|
#include <string.h> |
| 11 |
|
#include <fcntl.h> |
| 12 |
+ |
#include <errno.h> |
| 13 |
|
#include "rtio.h" |
| 14 |
|
#include "platform.h" |
| 15 |
|
#include "resolu.h" |
| 146 |
|
static int |
| 147 |
|
rmx_load_double(RMATRIX *rm, FILE *fp) |
| 148 |
|
{ |
| 149 |
< |
int i, j, k; |
| 149 |
< |
double val[100]; |
| 149 |
> |
int i, j; |
| 150 |
|
|
| 151 |
– |
if (rm->ncomp > 100) { |
| 152 |
– |
fputs("Unsupported # components in rmx_load_double()\n", stderr); |
| 153 |
– |
exit(1); |
| 154 |
– |
} |
| 151 |
|
for (i = 0; i < rm->nrows; i++) |
| 152 |
< |
for (j = 0; j < rm->ncols; j++) { |
| 153 |
< |
if (getbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp) |
| 152 |
> |
for (j = 0; j < rm->ncols; j++) |
| 153 |
> |
if (getbinary(&rmx_lval(rm,i,j,0), sizeof(double), rm->ncomp, fp) != rm->ncomp) |
| 154 |
|
return(0); |
| 159 |
– |
for (k = rm->ncomp; k--; ) |
| 160 |
– |
rmx_lval(rm,i,j,k) = val[k]; |
| 161 |
– |
} |
| 155 |
|
return(1); |
| 156 |
|
} |
| 157 |
|
|
| 326 |
|
static int |
| 327 |
|
rmx_write_double(const RMATRIX *rm, FILE *fp) |
| 328 |
|
{ |
| 329 |
< |
int i, j, k; |
| 337 |
< |
double val[100]; |
| 329 |
> |
int i, j; |
| 330 |
|
|
| 339 |
– |
if (rm->ncomp > 100) { |
| 340 |
– |
fputs("Unsupported # components in rmx_write_double()\n", stderr); |
| 341 |
– |
exit(1); |
| 342 |
– |
} |
| 331 |
|
for (i = 0; i < rm->nrows; i++) |
| 332 |
< |
for (j = 0; j < rm->ncols; j++) { |
| 333 |
< |
for (k = rm->ncomp; k--; ) |
| 346 |
< |
val[k] = rmx_lval(rm,i,j,k); |
| 347 |
< |
if (putbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp) |
| 332 |
> |
for (j = 0; j < rm->ncols; j++) |
| 333 |
> |
if (putbinary(&rmx_lval(rm,i,j,0), sizeof(double), rm->ncomp, fp) != rm->ncomp) |
| 334 |
|
return(0); |
| 349 |
– |
} |
| 335 |
|
return(1); |
| 336 |
|
} |
| 337 |
|
|
| 505 |
|
rmx_lval(mres,i,j,k) = (double)d; |
| 506 |
|
} |
| 507 |
|
return(mres); |
| 508 |
+ |
} |
| 509 |
+ |
|
| 510 |
+ |
/* Element-wise multiplication (or division) of m2 into m1 */ |
| 511 |
+ |
int |
| 512 |
+ |
rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide) |
| 513 |
+ |
{ |
| 514 |
+ |
int zeroDivides = 0; |
| 515 |
+ |
int i, j, k; |
| 516 |
+ |
|
| 517 |
+ |
if ((m1 == NULL) | (m2 == NULL) || |
| 518 |
+ |
(m1->ncols != m2->ncols) | (m1->nrows != m2->nrows)) |
| 519 |
+ |
return(0); |
| 520 |
+ |
if ((m2->ncomp > 1) & (m2->ncomp != m1->ncomp)) |
| 521 |
+ |
return(0); |
| 522 |
+ |
i = rmx_newtype(m1->dtype, m2->dtype); |
| 523 |
+ |
if (i) |
| 524 |
+ |
m1->dtype = i; |
| 525 |
+ |
else |
| 526 |
+ |
rmx_addinfo(m1, rmx_mismatch_warn); |
| 527 |
+ |
for (i = m1->nrows; i--; ) |
| 528 |
+ |
for (j = m1->ncols; j--; ) |
| 529 |
+ |
if (divide) { |
| 530 |
+ |
double d; |
| 531 |
+ |
if (m2->ncomp == 1) { |
| 532 |
+ |
d = rmx_lval(m2,i,j,0); |
| 533 |
+ |
if (d == 0) { |
| 534 |
+ |
++zeroDivides; |
| 535 |
+ |
for (k = m1->ncomp; k--; ) |
| 536 |
+ |
rmx_lval(m1,i,j,k) = 0; |
| 537 |
+ |
} else { |
| 538 |
+ |
d = 1./d; |
| 539 |
+ |
for (k = m1->ncomp; k--; ) |
| 540 |
+ |
rmx_lval(m1,i,j,k) *= d; |
| 541 |
+ |
} |
| 542 |
+ |
} else |
| 543 |
+ |
for (k = m1->ncomp; k--; ) { |
| 544 |
+ |
d = rmx_lval(m2,i,j,k); |
| 545 |
+ |
if (d == 0) { |
| 546 |
+ |
++zeroDivides; |
| 547 |
+ |
rmx_lval(m1,i,j,k) = 0; |
| 548 |
+ |
} else |
| 549 |
+ |
rmx_lval(m1,i,j,k) /= d; |
| 550 |
+ |
} |
| 551 |
+ |
} else { |
| 552 |
+ |
if (m2->ncomp == 1) { |
| 553 |
+ |
const double d = rmx_lval(m2,i,j,0); |
| 554 |
+ |
for (k = m1->ncomp; k--; ) |
| 555 |
+ |
rmx_lval(m1,i,j,k) *= d; |
| 556 |
+ |
} else |
| 557 |
+ |
for (k = m1->ncomp; k--; ) |
| 558 |
+ |
rmx_lval(m1,i,j,k) *= rmx_lval(m2,i,j,k); |
| 559 |
+ |
} |
| 560 |
+ |
if (zeroDivides) { |
| 561 |
+ |
fputs("Divide by zero in rmx_elemult()\n", stderr); |
| 562 |
+ |
errno = ERANGE; |
| 563 |
+ |
} |
| 564 |
+ |
return(1); |
| 565 |
|
} |
| 566 |
|
|
| 567 |
|
/* Sum second matrix into first, applying scale factor beforehand */ |
