The R Project SVN R

Rev

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

Rev 21875 Rev 23175
Line 320... Line 320...
320 
    char *transa = "N", *transb = "N";
 320 
    char *transa = "N", *transb = "N";
321 
    int i;
 321 
    int i;
322 
    double one = 1.0, zero = 0.0;
 322 
    double one = 1.0, zero = 0.0;
323 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 323 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
324 
        F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
 324 
        F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
325 
		    x, &nrx, y, &nry, &zero, z, &nrx);
 325 
			x, &nrx, y, &nry, &zero, z, &nrx);
326 
    }
 326 
    }
327 
    else { /* zero-extent operations should return zeroes */
 327 
    else { /* zero-extent operations should return zeroes */
328 
	for(i=0;i<nrx*ncy;i++)
 328 
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
329 
	    z[i]=0;
 - 
 
330 
    }
 329 
    }
331 
#else
 330 
#else
332 
 
 331 
 
333 
/* FIXME - What about non-IEEE overflow ??? */
 332 
/* FIXME - What about non-IEEE overflow ??? */
334 
/* Does it really matter? */
 333 
/* Does it really matter? */
Line 352... Line 351...
352 
	    ;
 351 
	    ;
353 
	}
 352 
	}
354 
#endif
 353 
#endif
355 
}
 354 
}
356 
 
 355 
 
- 
 
 356 
/* DGEMM - perform one of the matrix-matrix operations    */
- 
 
 357 
/* C := alpha*op( A )*op( B ) + beta*C */
- 
 
 358 
extern void
- 
 
 359 
F77_NAME(zgemm)(const char *transa, const char *transb, const int *m,
- 
 
 360 
		const int *n, const int *k, const Rcomplex *alpha,
- 
 
 361 
		const Rcomplex *a, const int *lda,
- 
 
 362 
		const Rcomplex *b, const int *ldb,
- 
 
 363 
		const Rcomplex *beta, Rcomplex *c, const int *ldc);
- 
 
 364 
 
357 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
 365 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
358 
		Rcomplex *y, int nry, int ncy, Rcomplex *z)
 366 
		     Rcomplex *y, int nry, int ncy, Rcomplex *z)
359 
{
 367 
{
- 
 
 368 
#ifdef IEEE_754
- 
 
 369 
    char *transa = "N", *transb = "N";
- 
 
 370 
    int i;
- 
 
 371 
    Rcomplex one, zero;
- 
 
 372 
 
- 
 
 373 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
- 
 
 374 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
- 
 
 375 
        F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
- 
 
 376 
			x, &nrx, y, &nry, &zero, z, &nrx);
- 
 
 377 
    } else { /* zero-extent operations should return zeroes */
- 
 
 378 
	for(i = 0; i < nrx*ncy; i++) z[i].r = z[i].i = 0;
- 
 
 379 
    }
- 
 
 380 
#else
360 
    int i, j, k;
 381 
    int i, j, k;
361 
    double xij_r, xij_i, yjk_r, yjk_i, sum_i, sum_r;
 382 
    double xij_r, xij_i, yjk_r, yjk_i, sum_i, sum_r;
362 
 
 383 
 
363 
    for (i = 0; i < nrx; i++)
 384 
    for (i = 0; i < nrx; i++)
364 
	for (k = 0; k < ncy; k++) {
 385 
	for (k = 0; k < ncy; k++) {
Line 369... Line 390...
369 
	    for (j = 0; j < ncx; j++) {
 390 
	    for (j = 0; j < ncx; j++) {
370 
		xij_r = x[i + j * nrx].r;
 391 
		xij_r = x[i + j * nrx].r;
371 
		xij_i = x[i + j * nrx].i;
 392 
		xij_i = x[i + j * nrx].i;
372 
		yjk_r = y[j + k * nry].r;
 393 
		yjk_r = y[j + k * nry].r;
373 
		yjk_i = y[j + k * nry].i;
 394 
		yjk_i = y[j + k * nry].i;
374 
#ifndef IEEE_754
 - 
 
375 
		if (ISNAN(xij_r) || ISNAN(xij_i)
 395 
		if (ISNAN(xij_r) || ISNAN(xij_i)
376 
		    || ISNAN(yjk_r) || ISNAN(yjk_i))
 396 
		    || ISNAN(yjk_r) || ISNAN(yjk_i))
377 
		    goto next_ik;
 397 
		    goto next_ik;
378 
#endif
 - 
 
379 
		sum_r += (xij_r * yjk_r - xij_i * yjk_i);
 398 
		sum_r += (xij_r * yjk_r - xij_i * yjk_i);
380 
		sum_i += (xij_r * yjk_i + xij_i * yjk_r);
 399 
		sum_i += (xij_r * yjk_i + xij_i * yjk_r);
381 
	    }
 400 
	    }
382 
	    z[i + k * nrx].r = sum_r;
 401 
	    z[i + k * nrx].r = sum_r;
383 
	    z[i + k * nrx].i = sum_i;
 402 
	    z[i + k * nrx].i = sum_i;
384 
#ifndef IEEE_754
 - 
 
385 
	next_ik:
 403 
	next_ik:
386 
	    ;
 404 
	    ;
387 
#endif
 - 
 
388 
	}
 405 
	}
- 
 
 406 
#endif
389 
}
 407 
}
390 
 
 408 
 
391 
static void symcrossprod(double *x, int nr, int nc, double *z)
 409 
static void symcrossprod(double *x, int nr, int nc, double *z)
392 
{
 410 
{
393 
    char *trans = "T", *uplo = "U";
 411 
    char *trans = "T", *uplo = "U";
394 
    double one = 1.0, zero = 0.0;
 412 
    double one = 1.0, zero = 0.0;
395 
    int i, j;
 413 
    int i, j;
396 
 
 - 
 
397 
    if (nr > 0 && nc > 0) {
 414 
    if (nr > 0 && nc > 0) {
398 
        F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
 415 
        F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
399 
	for (i = 1; i < nc; i++)
416 
	for (i = 1; i < nc; i++)
400 
	    for (j = 0; j < i; j++) z[i + nc *j] = z[j + nc * i];
 417 
	    for (j = 0; j < i; j++) z[i + nc *j] = z[j + nc * i];
401 
    }
 418 
    }
Line 407... Line 424...
407 
#ifdef IEEE_754
 424 
#ifdef IEEE_754
408 
    char *transa = "T", *transb = "N";
 425 
    char *transa = "T", *transb = "N";
409 
    double one = 1.0, zero = 0.0;
 426 
    double one = 1.0, zero = 0.0;
410 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 427 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
411 
        F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
 428 
        F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
412 
		    x, &nrx, y, &nry, &zero, z, &ncx);
 429 
			x, &nrx, y, &nry, &zero, z, &ncx);
413 
    }
 430 
    }
414 
#else
 431 
#else
415 
    int i, j, k;
 432 
    int i, j, k;
416 
    double xji, yjk, sum;
 433 
    double xji, yjk, sum;
417 
 
 434 
 
Line 434... Line 451...
434 
}
 451 
}
435 
 
 452 
 
436 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
 453 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
437 
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
 454 
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
438 
{
 455 
{
- 
 
 456 
#ifdef IEEE_754
- 
 
 457 
    char *transa = "T", *transb = "N";
- 
 
 458 
    Rcomplex one, zero;
- 
 
 459 
 
- 
 
 460 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
- 
 
 461 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
- 
 
 462 
        F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
- 
 
 463 
			x, &nrx, y, &nry, &zero, z, &ncx);
- 
 
 464 
    }
- 
 
 465 
#else
439 
    int i, j, k;
 466 
    int i, j, k;
440 
    double xji_r, xji_i, yjk_r, yjk_i, sum_r, sum_i;
 467 
    double xji_r, xji_i, yjk_r, yjk_i, sum_r, sum_i;
441 
 
 468 
 
442 
    for (i = 0; i < ncx; i++)
 469 
    for (i = 0; i < ncx; i++)
443 
	for (k = 0; k < ncy; k++) {
 470 
	for (k = 0; k < ncy; k++) {
Line 463... Line 490...
463 
#ifndef IEEE_754
 490 
#ifndef IEEE_754
464 
	next_ik:
 491 
	next_ik:
465 
	    ;
 492 
	    ;
466 
#endif
 493 
#endif
467 
	}
 494 
	}
- 
 
 495 
#endif
468 
}
 496 
}
469 
 
 - 
 
470 
SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
 497 
SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
471 
{
 498 
{
472 
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
 499 
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
473 
    SEXP x = CAR(args), y = CADR(args), xdims, ydims, ans;
 500 
    SEXP x = CAR(args), y = CADR(args), xdims, ydims, ans;
474 
    Rboolean sym;
 501 
    Rboolean sym;
Line 609... Line 636...
609 
	}
 636 
	}
610 
    }
 637 
    }
611 
    else {
 638 
    else {
612 
	PROTECT(ans = allocMatrix(mode, ncx, ncy));
 639 
	PROTECT(ans = allocMatrix(mode, ncx, ncy));
613 
	if (mode == CPLXSXP)
 640 
	if (mode == CPLXSXP)
- 
 
 641 
	    if(sym)
614 
	    ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
 642 
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
- 
 
 643 
			   COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans));
- 
 
 644 
	    else
- 
 
 645 
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
615 
		       COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
 646 
			   COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
616 
	else {
 647 
	else {
617 
#ifdef IEEE_754
 648 
#ifdef IEEE_754
618 
	    if(sym)
 649 
	    if(sym)
619 
		symcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));
 650 
		symcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));
620 
	    else
 651 
	    else