The R Project SVN R

Rev

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

Rev 26666 Rev 27297
Line 313... Line 313...
313 
    return ans;
 313 
    return ans;
314 
}
 314 
}
315 
 
 315 
 
316 
static void matprod(double *x, int nrx, int ncx,
 316 
static void matprod(double *x, int nrx, int ncx,
317 
		    double *y, int nry, int ncy, double *z)
 317 
		    double *y, int nry, int ncy, double *z)
318 
{
 - 
 
319 
#ifdef IEEE_754
 318 
#ifdef IEEE_754
- 
 
 319 
{
320 
    char *transa = "N", *transb = "N";
 320 
    char *transa = "N", *transb = "N";
321 
    int i;
 321 
    int i,  j, k;
322 
    double one = 1.0, zero = 0.0;
 322 
    double one = 1.0, zero = 0.0, sum;
- 
 
 323 
    Rboolean have_na = FALSE;
- 
 
 324 
 
323 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 325 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
- 
 
 326 
	/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582
- 
 
 327 
	 * The test is only O(n) here
- 
 
 328 
	 */
- 
 
 329 
	for (i = 0; i < nrx*ncx; i++)
- 
 
 330 
	    if (ISNAN(x[i])) {have_na = TRUE; break;}
- 
 
 331 
	if (!have_na)
- 
 
 332 
	    for (i = 0; i < nry*ncy; i++)
- 
 
 333 
		if (ISNAN(y[i])) {have_na = TRUE; break;}
- 
 
 334 
	if (have_na) {
- 
 
 335 
	    for (i = 0; i < nrx; i++)
- 
 
 336 
		for (k = 0; k < ncy; k++) {
- 
 
 337 
		    sum = 0.0;
- 
 
 338 
		    for (j = 0; j < ncx; j++)
- 
 
 339 
			sum += x[i + j * nrx] * y[j + k * nry];
- 
 
 340 
		    z[i + k * nrx] = sum;
- 
 
 341 
		}
- 
 
 342 
	} else
324 
        F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
 343 
	    F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
325 
			x, &nrx, y, &nry, &zero, z, &nrx);
 344 
			    x, &nrx, y, &nry, &zero, z, &nrx);
326 
    }
 - 
 
327 
    else { /* zero-extent operations should return zeroes */
 345 
    } else /* zero-extent operations should return zeroes */
328 
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
 346 
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
329 
    }
 347 
}
330 
#else
 348 
#else
331 
 
 349 
{
332 
/* FIXME - What about non-IEEE overflow ??? */
 350 
/* FIXME - What about non-IEEE overflow ??? */
333 
/* Does it really matter? */
 351 
/* Does it really matter? */
334 
 
 352 
 
335 
    int i, j, k;
 353 
    int i, j, k;
336 
    double xij, yjk, sum;
 354 
    double xij, yjk, sum;
337 
 
 355 
 
- 
 
 356 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
338 
    for (i = 0; i < nrx; i++)
 357 
	for (i = 0; i < nrx; i++)
339 
	for (k = 0; k < ncy; k++) {
 358 
	    for (k = 0; k < ncy; k++) {
340 
	    z[i + k * nrx] = NA_REAL;
 359 
		z[i + k * nrx] = NA_REAL;
341 
	    sum = 0.0;
 360 
		sum = 0.0;
342 
	    for (j = 0; j < ncx; j++) {
 361 
		for (j = 0; j < ncx; j++) {
343 
		xij = x[i + j * nrx];
 362 
		    xij = x[i + j * nrx];
344 
		yjk = y[j + k * nry];
 363 
		    yjk = y[j + k * nry];
345 
		if (ISNAN(xij) || ISNAN(yjk))
 364 
		    if (ISNAN(xij) || ISNAN(yjk)) goto next_ik;
346 
		    goto next_ik;
 365 
		    sum += xij * yjk;
- 
 
 366 
		}
347 
		sum += xij * yjk;
 367 
		z[i + k * nrx] = sum;
- 
 
 368 
	    next_ik:
- 
 
 369 
		;
348 
	    }
 370 
	    }
- 
 
 371 
    } else /* zero-extent operations should return zeroes */
349 
	    z[i + k * nrx] = sum;
 372 
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
350 
	next_ik:
 - 
 
351 
	    ;
 - 
 
352 
	}
 - 
 
353 
#endif
 - 
 
354 
}
 373 
}
- 
 
 374 
#endif
355 
 
 375 
 
356 
#ifdef HAVE_DOUBLE_COMPLEX
 376 
#ifdef HAVE_DOUBLE_COMPLEX
357 
/* ZGEMM - perform one of the matrix-matrix operations    */
 377 
/* ZGEMM - perform one of the matrix-matrix operations    */
358 
/* C := alpha*op( A )*op( B ) + beta*C */
 378 
/* C := alpha*op( A )*op( B ) + beta*C */
359 
extern void
 379 
extern void