The R Project SVN R

Rev

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

Rev 13214 Rev 13942
Line 22... Line 22...
22 
#include <config.h>
 22 
#include <config.h>
23 
#endif
 23 
#endif
24 
 
 24 
 
25 
#include <Defn.h>
 25 
#include <Defn.h>
26 
#include <Rmath.h>
 26 
#include <Rmath.h>
27 
 
 - 
 
- 
 
 27 
#include <R_ext/RS.h>
28 
/* "GetRowNames" and "GetColNames" are utility routines which */
 28 
/* "GetRowNames" and "GetColNames" are utility routines which */
29 
/* locate and return the row names and column names from the */
 29 
/* locate and return the row names and column names from the */
30 
/* dimnames attribute of a matrix.  They are useful because */
 30 
/* dimnames attribute of a matrix.  They are useful because */
31 
/* old versions of R used pair-based lists for dimnames */
 31 
/* old versions of R used pair-based lists for dimnames */
32 
/* whereas recent versions use vector based lists */
 32 
/* whereas recent versions use vector based lists */
Line 304... Line 304...
304 
	break;
 304 
	break;
305 
    }
 305 
    }
306 
    return ans;
 306 
    return ans;
307 
}
 307 
}
308 
 
 308 
 
309 
/* FIXME - What about non-IEEE overflow ??? */
 - 
 
310 
/* Does it really matter? */
 - 
 
311 
 
 - 
 
312 
static void matprod(double *x, int nrx, int ncx,
 309 
static void matprod(double *x, int nrx, int ncx,
313 
		    double *y, int nry, int ncy, double *z)
 310 
		    double *y, int nry, int ncy, double *z)
314 
{
 311 
{
- 
 
 312 
#ifdef IEEE_754
- 
 
 313 
    char *transa = "N", *transb = "N";
- 
 
 314 
    double one = 1.0, zero = 0.0;
- 
 
 315 
    F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
- 
 
 316 
		    x, &nrx, y, &nry, &zero, z, &nrx);
- 
 
 317 
#else
- 
 
 318 
 
- 
 
 319 
/* FIXME - What about non-IEEE overflow ??? */
- 
 
 320 
/* Does it really matter? */
- 
 
 321 
 
315 
    int i, j, k;
 322 
    int i, j, k;
316 
    double xij, yjk, sum;
 323 
    double xij, yjk, sum;
317 
 
 324 
 
318 
    for (i = 0; i < nrx; i++)
 325 
    for (i = 0; i < nrx; i++)
319 
	for (k = 0; k < ncy; k++) {
 326 
	for (k = 0; k < ncy; k++) {
320 
	    z[i + k * nrx] = NA_REAL;
 327 
	    z[i + k * nrx] = NA_REAL;
321 
	    sum = 0.0;
 328 
	    sum = 0.0;
322 
	    for (j = 0; j < ncx; j++) {
 329 
	    for (j = 0; j < ncx; j++) {
323 
		xij = x[i + j * nrx];
 330 
		xij = x[i + j * nrx];
324 
		yjk = y[j + k * nry];
 331 
		yjk = y[j + k * nry];
325 
#ifndef IEEE_754
 - 
 
326 
		if (ISNAN(xij) || ISNAN(yjk))
 332 
		if (ISNAN(xij) || ISNAN(yjk))
327 
		    goto next_ik;
 333 
		    goto next_ik;
328 
#endif
 - 
 
329 
		sum += xij * yjk;
 334 
		sum += xij * yjk;
330 
	    }
 335 
	    }
331 
	    z[i + k * nrx] = sum;
 336 
	    z[i + k * nrx] = sum;
332 
#ifndef IEEE_754
 - 
 
333 
	next_ik:
 337 
	next_ik:
334 
	    ;
 338 
	    ;
335 
#endif
 - 
 
336 
	}
 339 
	}
- 
 
 340 
#endif
337 
}
 341 
}
338 
 
 342 
 
339 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
 343 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
340 
		Rcomplex *y, int nry, int ncy, Rcomplex *z)
 344 
		Rcomplex *y, int nry, int ncy, Rcomplex *z)
341 
{
 345 
{
Line 371... Line 375...
371 
}
 375 
}
372 
 
 376 
 
373 
static void crossprod(double *x, int nrx, int ncx,
 377 
static void crossprod(double *x, int nrx, int ncx,
374 
		      double *y, int nry, int ncy, double *z)
 378 
		      double *y, int nry, int ncy, double *z)
375 
{
 379 
{
- 
 
 380 
#ifdef IEEE_754
- 
 
 381 
    char *transa = "T", *transb = "N";
- 
 
 382 
    double one = 1.0, zero = 0.0;
- 
 
 383 
    F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
- 
 
 384 
		    x, &nrx, y, &nry, &zero, z, &ncx);
- 
 
 385 
#else
376 
    int i, j, k;
 386 
    int i, j, k;
377 
    double xji, yjk, sum;
 387 
    double xji, yjk, sum;
378 
 
 388 
 
379 
    for (i = 0; i < ncx; i++)
 389 
    for (i = 0; i < ncx; i++)
380 
	for (k = 0; k < ncy; k++) {
 390 
	for (k = 0; k < ncy; k++) {
381 
	    z[i + k * ncx] = NA_REAL;
 391 
	    z[i + k * ncx] = NA_REAL;
382 
	    sum = 0.0;
 392 
	    sum = 0.0;
383 
	    for (j = 0; j < nrx; j++) {
 393 
	    for (j = 0; j < nrx; j++) {
384 
		xji = x[j + i * nrx];
 394 
		xji = x[j + i * nrx];
385 
		yjk = y[j + k * nry];
 395 
		yjk = y[j + k * nry];
386 
#ifndef IEEE_754
 - 
 
387 
		if (ISNAN(xji) || ISNAN(yjk))
 396 
		if (ISNAN(xji) || ISNAN(yjk))
388 
		    goto next_ik;
 397 
		    goto next_ik;
389 
#endif
 - 
 
390 
		sum += xji * yjk;
 398 
		sum += xji * yjk;
391 
	    }
 399 
	    }
392 
	    z[i + k * ncx] = sum;
 400 
	    z[i + k * ncx] = sum;
393 
#ifndef IEEE_754
 - 
 
394 
	next_ik:
 401 
	next_ik:
395 
	    ;
 402 
	    ;
396 
#endif
 - 
 
397 
	}
 403 
	}
- 
 
 404 
#endif
398 
}
 405 
}
399 
 
 406 
 
400 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
 407 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
401 
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
 408 
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
402 
{
 409 
{