The R Project SVN R

Rev

Rev 29387 | Rev 29902 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 29387 Rev 29410
Line 371... Line 371...
371 
    return ans;
 371 
    return ans;
372 
}
 372 
}
373 
 
 373 
 
374 
static void matprod(double *x, int nrx, int ncx,
 374 
static void matprod(double *x, int nrx, int ncx,
375 
		    double *y, int nry, int ncy, double *z)
 375 
		    double *y, int nry, int ncy, double *z)
376 
#ifdef IEEE_754
 - 
 
377 
{
 376 
{
378 
    char *transa = "N", *transb = "N";
 377 
    char *transa = "N", *transb = "N";
379 
    int i,  j, k;
 378 
    int i,  j, k;
380 
    double one = 1.0, zero = 0.0, sum;
 379 
    double one = 1.0, zero = 0.0, sum;
381 
    Rboolean have_na = FALSE;
 380 
    Rboolean have_na = FALSE;
Line 401... Line 400...
401 
	    F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
 400 
	    F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
402 
			    x, &nrx, y, &nry, &zero, z, &nrx);
 401 
			    x, &nrx, y, &nry, &zero, z, &nrx);
403 
    } else /* zero-extent operations should return zeroes */
 402 
    } else /* zero-extent operations should return zeroes */
404 
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
 403 
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
405 
}
 404 
}
406 
#else
 - 
 
407 
{
 - 
 
408 
/* FIXME - What about non-IEEE overflow ??? */
 - 
 
409 
/* Does it really matter? */
 - 
 
410 
 
 - 
 
411 
    int i, j, k;
 - 
 
412 
    double xij, yjk, sum;
 - 
 
413 
 
 - 
 
414 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 - 
 
415 
	for (i = 0; i < nrx; i++)
 - 
 
416 
	    for (k = 0; k < ncy; k++) {
 - 
 
417 
		z[i + k * nrx] = NA_REAL;
 - 
 
418 
		sum = 0.0;
 - 
 
419 
		for (j = 0; j < ncx; j++) {
 - 
 
420 
		    xij = x[i + j * nrx];
 - 
 
421 
		    yjk = y[j + k * nry];
 - 
 
422 
		    if (ISNAN(xij) || ISNAN(yjk)) goto next_ik;
 - 
 
423 
		    sum += xij * yjk;
 - 
 
424 
		}
 - 
 
425 
		z[i + k * nrx] = sum;
 - 
 
426 
	    next_ik:
 - 
 
427 
		;
 - 
 
428 
	    }
 - 
 
429 
    } else /* zero-extent operations should return zeroes */
 - 
 
430 
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
 - 
 
431 
}
 - 
 
432 
#endif
 - 
 
433 
 
 405 
 
434 
#ifdef HAVE_DOUBLE_COMPLEX
 406 
#ifdef HAVE_DOUBLE_COMPLEX
435 
/* ZGEMM - perform one of the matrix-matrix operations    */
 407 
/* ZGEMM - perform one of the matrix-matrix operations    */
436 
/* C := alpha*op( A )*op( B ) + beta*C */
 408 
/* C := alpha*op( A )*op( B ) + beta*C */
437 
extern void
 409 
extern void
Line 443... Line 415...
443 
#endif
 415 
#endif
444 
 
 416 
 
445 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
 417 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
446 
		     Rcomplex *y, int nry, int ncy, Rcomplex *z)
 418 
		     Rcomplex *y, int nry, int ncy, Rcomplex *z)
447 
{
 419 
{
448 
#if defined(HAVE_DOUBLE_COMPLEX) && defined(IEEE_754)
 420 
#ifdef HAVE_DOUBLE_COMPLEX
449 
    char *transa = "N", *transb = "N";
 421 
    char *transa = "N", *transb = "N";
450 
    int i;
 422 
    int i;
451 
    Rcomplex one, zero;
 423 
    Rcomplex one, zero;
452 
 
 424 
 
453 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
 425 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
Line 502... Line 474...
502 
}
 474 
}
503 
 
 475 
 
504 
static void crossprod(double *x, int nrx, int ncx,
 476 
static void crossprod(double *x, int nrx, int ncx,
505 
		      double *y, int nry, int ncy, double *z)
 477 
		      double *y, int nry, int ncy, double *z)
506 
{
 478 
{
507 
#ifdef IEEE_754
 - 
 
508 
    char *transa = "T", *transb = "N";
 479 
    char *transa = "T", *transb = "N";
509 
    double one = 1.0, zero = 0.0;
 480 
    double one = 1.0, zero = 0.0;
510 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 481 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
511 
        F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
 482 
        F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
512 
			x, &nrx, y, &nry, &zero, z, &ncx);
 483 
			x, &nrx, y, &nry, &zero, z, &ncx);
513 
    }
 484 
    }
514 
    else { /* zero-extent operations should return zeroes */
 485 
    else { /* zero-extent operations should return zeroes */
515 
	int i;
 486 
	int i;
516 
	for(i = 0; i < ncx*ncy; i++) z[i] = 0;
 487 
	for(i = 0; i < ncx*ncy; i++) z[i] = 0;
517 
    }
 488 
    }
518 
#else
 - 
 
519 
    int i, j, k;
 - 
 
520 
    double xji, yjk, sum;
 - 
 
521 
 
 - 
 
522 
    for (i = 0; i < ncx; i++)
 - 
 
523 
	for (k = 0; k < ncy; k++) {
 - 
 
524 
	    z[i + k * ncx] = NA_REAL;
 - 
 
525 
	    sum = 0.0;
 - 
 
526 
	    for (j = 0; j < nrx; j++) {
 - 
 
527 
		xji = x[j + i * nrx];
 - 
 
528 
		yjk = y[j + k * nry];
 - 
 
529 
		if (ISNAN(xji) || ISNAN(yjk))
 - 
 
530 
		    goto next_ik;
 - 
 
531 
		sum += xji * yjk;
 - 
 
532 
	    }
 - 
 
533 
	    z[i + k * ncx] = sum;
 - 
 
534 
	next_ik:
 - 
 
535 
	    ;
 - 
 
536 
	}
 - 
 
537 
#endif
 - 
 
538 
}
 489 
}
539 
 
 490 
 
540 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
 491 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
541 
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
 492 
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
542 
{
 493 
{
543 
#ifdef IEEE_754
 - 
 
544 
    char *transa = "T", *transb = "N";
 494 
    char *transa = "T", *transb = "N";
545 
    Rcomplex one, zero;
 495 
    Rcomplex one, zero;
546 
 
 496 
 
547 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
 497 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
548 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 498 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
Line 551... Line 501...
551 
    }
 501 
    }
552 
    else { /* zero-extent operations should return zeroes */
 502 
    else { /* zero-extent operations should return zeroes */
553 
	int i;
 503 
	int i;
554 
	for(i = 0; i < ncx*ncy; i++) z[i].r = z[i].i = 0;
 504 
	for(i = 0; i < ncx*ncy; i++) z[i].r = z[i].i = 0;
555 
    }
 505 
    }
556 
#else
 - 
 
557 
    int i, j, k;
 - 
 
558 
    double xji_r, xji_i, yjk_r, yjk_i, sum_r, sum_i;
 - 
 
559 
 
 - 
 
560 
    for (i = 0; i < ncx; i++)
 - 
 
561 
	for (k = 0; k < ncy; k++) {
 - 
 
562 
	    z[i + k * ncx].r = NA_REAL;
 - 
 
563 
	    z[i + k * ncx].i = NA_REAL;
 - 
 
564 
	    sum_r = 0.0;
 - 
 
565 
	    sum_i = 0.0;
 - 
 
566 
	    for (j = 0; j < nrx; j++) {
 - 
 
567 
		xji_r = x[j + i * nrx].r;
 - 
 
568 
		xji_i = x[j + i * nrx].i;
 - 
 
569 
		yjk_r = y[j + k * nry].r;
 - 
 
570 
		yjk_i = y[j + k * nry].i;
 - 
 
571 
#ifndef IEEE_754
 - 
 
572 
		if (ISNAN(xji_r) || ISNAN(xji_i)
 - 
 
573 
		    || ISNAN(yjk_r) || ISNAN(yjk_i))
 - 
 
574 
		    goto next_ik;
 - 
 
575 
#endif
 - 
 
576 
		sum_r += (xji_r * yjk_r - xji_i * yjk_i);
 - 
 
577 
		sum_i += (xji_r * yjk_i + xji_i * yjk_r);
 - 
 
578 
	    }
 - 
 
579 
	    z[i + k * ncx].r = sum_r;
 - 
 
580 
	    z[i + k * ncx].i = sum_i;
 - 
 
581 
#ifndef IEEE_754
 - 
 
582 
	next_ik:
 - 
 
583 
	    ;
 - 
 
584 
#endif
 - 
 
585 
	}
 - 
 
586 
#endif
 - 
 
587 
}
 506 
}
588 
/* "%*%" (op = 0)  or  crossprod (op = 1) : */
 507 
/* "%*%" (op = 0)  or  crossprod (op = 1) : */
589 
SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
 508 
SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
590 
{
 509 
{
591 
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
 510 
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
Line 740... Line 659...
740 
			   COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans));
 659 
			   COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans));
741 
	    else
 660 
	    else
742 
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
 661 
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
743 
			   COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
 662 
			   COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
744 
	else {
 663 
	else {
745 
#ifdef IEEE_754
 - 
 
746 
	    if(sym)
 664 
	    if(sym)
747 
		symcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));
 665 
		symcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));
748 
	    else
 666 
	    else
749 
		crossprod(REAL(CAR(args)), nrx, ncx,
 667 
		crossprod(REAL(CAR(args)), nrx, ncx,
750 
			  REAL(CADR(args)), nry, ncy, REAL(ans));
 668 
			  REAL(CADR(args)), nry, ncy, REAL(ans));
751 
#endif
 - 
 
752 
	    if(sym)
 - 
 
753 
		crossprod(REAL(CAR(args)), nrx, ncx,
 - 
 
754 
			  REAL(CAR(args)), nry, ncy, REAL(ans));
 - 
 
755 
	    else
 - 
 
756 
		crossprod(REAL(CAR(args)), nrx, ncx,
 - 
 
757 
			  REAL(CADR(args)), nry, ncy, REAL(ans));
 - 
 
758 
	}
 669 
	}
759 
 
 670 
 
760 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
 671 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
761 
	if (sym)
 672 
	if (sym)
762 
	    PROTECT(ydims = xdims);
 673 
	    PROTECT(ydims = xdims);
Line 1109... Line 1020...
1109 
	PROTECT(ans = allocVector(REALSXP, p));
 1020 
	PROTECT(ans = allocVector(REALSXP, p));
1110 
	for (j = 0; j < p; j++) {
 1021 
	for (j = 0; j < p; j++) {
1111 
	    switch (type) {
 1022 
	    switch (type) {
1112 
	    case REALSXP:
 1023 
	    case REALSXP:
1113 
		rx = REAL(x) + n*j;
 1024 
		rx = REAL(x) + n*j;
1114 
#ifdef IEEE_754
 - 
 
1115 
		if (keepNA)
 1025 
		if (keepNA)
1116 
		    for (sum = 0., i = 0; i < n; i++) sum += *rx++;
 1026 
		    for (sum = 0., i = 0; i < n; i++) sum += *rx++;
1117 
		else {
 1027 
		else {
1118 
		    for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
 1028 
		    for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
1119 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
 1029 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
1120 
			else if (keepNA) {sum = NA_REAL; break;}
 1030 
			else if (keepNA) {sum = NA_REAL; break;}
1121 
		}
 1031 
		}
1122 
#else
 - 
 
1123 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
 - 
 
1124 
		    if (!ISNAN(*rx)) {cnt++; sum += *rx;}
 - 
 
1125 
		    else if (keepNA) {sum = NA_REAL; break;}
 - 
 
1126 
#endif
 - 
 
1127 
		break;
 1032 
		break;
1128 
	    case INTSXP:
 1033 
	    case INTSXP:
1129 
		ix = INTEGER(x) + n*j;
 1034 
		ix = INTEGER(x) + n*j;
1130 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
 1035 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
1131 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
 1036 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
Line 1147... Line 1052...
1147 
 
 1052 
 
1148 
    if (OP == 2 || OP == 3) { /* rows */
 1053 
    if (OP == 2 || OP == 3) { /* rows */
1149 
	cnt = p;
 1054 
	cnt = p;
1150 
	PROTECT(ans = allocVector(REALSXP, n));
 1055 
	PROTECT(ans = allocVector(REALSXP, n));
1151 
 
 1056 
 
1152 
#ifdef IEEE_754
 - 
 
1153 
	/* reverse summation order to improve cache hits */
 1057 
	/* reverse summation order to improve cache hits */
1154 
	if (type == REALSXP) {
 1058 
	if (type == REALSXP) {
1155 
	    double *rans = REAL(ans), *ra = rans, *cnt = NULL, *c;
 1059 
	    double *rans = REAL(ans), *ra = rans, *cnt = NULL, *c;
1156 
	    rx = REAL(x);
 1060 
	    rx = REAL(x);
1157 
	    if (!keepNA && OP == 3) cnt = Calloc(n, double);
 1061 
	    if (!keepNA && OP == 3) cnt = Calloc(n, double);
Line 1178... Line 1082...
1178 
		}
 1082 
		}
1179 
	    }
 1083 
	    }
1180 
	    UNPROTECT(1);
 1084 
	    UNPROTECT(1);
1181 
	    return ans;
 1085 
	    return ans;
1182 
	}
 1086 
	}
1183 
#endif
 - 
 
1184 
 
 1087 
 
1185 
	for (i = 0; i < n; i++) {
 1088 
	for (i = 0; i < n; i++) {
1186 
	    switch (type) {
 1089 
	    switch (type) {
1187 
	    case REALSXP:
 1090 
	    case REALSXP: /* this cannot be reached */
1188 
		rx = REAL(x) + i;
 1091 
		rx = REAL(x) + i;
1189 
#ifdef IEEE_754
 - 
 
1190 
		if (keepNA)
 1092 
		if (keepNA)
1191 
		    for (sum = 0., j = 0; j < p; j++, rx += n) sum += *rx;
 1093 
		    for (sum = 0., j = 0; j < p; j++, rx += n) sum += *rx;
1192 
		else {
 1094 
		else {
1193 
		    for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n)
 1095 
		    for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n)
1194 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
 1096 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
1195 
			else if (keepNA) {sum = NA_REAL; break;}
 1097 
			else if (keepNA) {sum = NA_REAL; break;}
1196 
		}
 1098 
		}
1197 
#else
 - 
 
1198 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n)
 - 
 
1199 
		    if (!ISNAN(*rx)) {cnt++; sum += *rx;}
 - 
 
1200 
		    else if (keepNA) {sum = NA_REAL; break;}
 - 
 
1201 
#endif
 - 
 
1202 
		break;
 1099 
		break;
1203 
	    case INTSXP:
 1100 
	    case INTSXP:
1204 
		ix = INTEGER(x) + i;
 1101 
		ix = INTEGER(x) + i;
1205 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
 1102 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
1206 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
 1103 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}