The R Project SVN R

Rev

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

Rev 17521 Rev 18753
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  R : A Computer Language for Statistical Data Analysis
 2 
 *  R : A Computer Language for Statistical Data Analysis
3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
4 
 *  Copyright (C) 1998-2001   The R Development Core Team.
 4 
 *  Copyright (C) 1998-2002   The R Development Core Team.
5 
 *
 5 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 6 
 *  This program is free software; you can redistribute it and/or modify
7 
 *  it under the terms of the GNU General Public License as published by
 7 
 *  it under the terms of the GNU General Public License as published by
8 
 *  the Free Software Foundation; either version 2 of the License, or
 8 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 9 
 *  (at your option) any later version.
Line 386... Line 386...
386 
	    ;
 386 
	    ;
387 
#endif
 387 
#endif
388 
	}
 388 
	}
389 
}
 389 
}
390 
 
 390 
 
- 
 
 391 
static void symcrossprod(double *x, int nr, int nc, double *z)
- 
 
 392 
{
- 
 
 393 
    char *trans = "T", *uplo = "U";
- 
 
 394 
    double one = 1.0, zero = 0.0;
- 
 
 395 
    int i, j;
- 
 
 396 
 
- 
 
 397 
    if (nr > 0 && nc > 0) {
- 
 
 398 
        F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
- 
 
 399 
	for (i = 1; i < nc; i++)
- 
 
 400 
	    for (j = 0; j < i; j++) z[i + nc *j] = z[j + nc * i];
- 
 
 401 
    }
- 
 
 402 
}
- 
 
 403 
 
391 
static void crossprod(double *x, int nrx, int ncx,
 404 
static void crossprod(double *x, int nrx, int ncx,
392 
		      double *y, int nry, int ncy, double *z)
 405 
		      double *y, int nry, int ncy, double *z)
393 
{
 406 
{
394 
#ifdef IEEE_754
 407 
#ifdef IEEE_754
395 
    char *transa = "T", *transb = "N";
 408 
    char *transa = "T", *transb = "N";
Line 455... Line 468...
455 
}
 468 
}
456 
 
 469 
 
457 
SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
 470 
SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
458 
{
 471 
{
459 
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
 472 
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
460 
    SEXP x, y, xdims, ydims, ans;
 473 
    SEXP x = CAR(args), y = CADR(args), xdims, ydims, ans;
- 
 
 474 
    Rboolean sym;
461 
 
 475 
 
- 
 
 476 
    sym = isNull(y);
462 
    if (!(isNumeric(CAR(args)) || isComplex(CAR(args))) ||
 477 
    if (sym && (PRIMVAL(op) == 1)) y = x;
463 
	!(isNumeric(CADR(args)) || isComplex(CADR(args))))
 478 
    if ( !(isNumeric(x) || isComplex(x)) || !(isNumeric(y) || isComplex(y)) )
464 
	errorcall(call, "requires numeric matrix/vector arguments");
 479 
	errorcall(call, "requires numeric matrix/vector arguments");
465 
 
 480 
 
466 
    x = CAR(args);
 - 
 
467 
    y = CADR(args);
 - 
 
468 
    xdims = getAttrib(x, R_DimSymbol);
 481 
    xdims = getAttrib(x, R_DimSymbol);
469 
    ydims = getAttrib(y, R_DimSymbol);
 482 
    ydims = getAttrib(y, R_DimSymbol);
470 
    ldx = length(xdims);
 483 
    ldx = length(xdims);
471 
    ldy = length(ydims);
 484 
    ldy = length(ydims);
472 
 
 485 
 
Line 592... Line 605...
592 
    else {
 605 
    else {
593 
	PROTECT(ans = allocMatrix(mode, ncx, ncy));
 606 
	PROTECT(ans = allocMatrix(mode, ncx, ncy));
594 
	if (mode == CPLXSXP)
 607 
	if (mode == CPLXSXP)
595 
	    ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
 608 
	    ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
596 
		       COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
 609 
		       COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
597 
	else
 610 
	else {
- 
 
 611 
#ifdef IEEE_754
- 
 
 612 
	    if(sym)
- 
 
 613 
		symcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));
- 
 
 614 
	    else
- 
 
 615 
#endif
598 
	    crossprod(REAL(CAR(args)), nrx, ncx,
 616 
		crossprod(REAL(CAR(args)), nrx, ncx,
599 
		      REAL(CADR(args)), nry, ncy, REAL(ans));
 617 
			  REAL(CADR(args)), nry, ncy, REAL(ans));
- 
 
 618 
	}
600 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
 619 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
601 
	PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
 620 
	PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
602 
	if (xdims != R_NilValue || ydims != R_NilValue) {
 621 
	if (xdims != R_NilValue || ydims != R_NilValue) {
603 
	    SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;
 622 
	    SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;
604 
	    PROTECT(dimnames = allocVector(VECSXP, 2));
 623 
	    PROTECT(dimnames = allocVector(VECSXP, 2));