The R Project SVN R

Rev

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

Rev 300 Rev 571
Line 313... Line 313...
313 
		break;
 313 
		break;
314 
	}
 314 
	}
315 
	return ans;
 315 
	return ans;
316 
}
 316 
}
317 
 
 317 
 
- 
 
 318 
/* FIXME - What about non IEEE overflow ??? */
- 
 
 319 
 
318 
static void matprod(double *x, int nrx, int ncx, double *y, int nry, int ncy, double *z)
 320 
static void matprod(double *x, int nrx, int ncx, double *y, int nry, int ncy, double *z)
319 
{
 321 
{
320 
	int i, j, k;
 322 
	int i, j, k;
321 
	double xij, yjk, sum;
 323 
	double xij, yjk, sum;
322 
 
 324 
 
Line 325... Line 327...
325 
			z[i + k * nrx] = NA_REAL;
 327 
			z[i + k * nrx] = NA_REAL;
326 
			sum = 0.0;
 328 
			sum = 0.0;
327 
			for (j = 0; j < ncx; j++) {
 329 
			for (j = 0; j < ncx; j++) {
328 
				xij = x[i + j * nrx];
 330 
				xij = x[i + j * nrx];
329 
				yjk = y[j + k * nry];
 331 
				yjk = y[j + k * nry];
- 
 
 332 
#ifndef IEEE_754
330 
				if (!FINITE(xij) || !FINITE(yjk))
 333 
				if (NAN(xij) || NAN(yjk))
331 
					goto next_ik;
 334 
					goto next_ik;
- 
 
 335 
#endif
332 
				sum += xij * yjk;
 336 
				sum += xij * yjk;
333 
			}
 337 
			}
334 
			z[i + k * nrx] = sum;
 338 
			z[i + k * nrx] = sum;
335 
		next_ik:
 339 
		next_ik:
336 
			;
 340 
			;
Line 352... Line 356...
352 
			for (j=0; j<ncx; j++) {
 356 
			for (j=0; j<ncx; j++) {
353 
				xij_r = x[i+j*nrx].r;
 357 
				xij_r = x[i+j*nrx].r;
354 
				xij_i = x[i+j*nrx].i;
 358 
				xij_i = x[i+j*nrx].i;
355 
				yjk_r = y[j+k*nry].r;
 359 
				yjk_r = y[j+k*nry].r;
356 
				yjk_i = y[j+k*nry].i;
 360 
				yjk_i = y[j+k*nry].i;
- 
 
 361 
#ifndef IEEE_754
357 
				if (!FINITE(xij_r) || !FINITE(xij_i)
 362 
				if (NAN(xij_r) || NAN(xij_i)
358 
						|| !FINITE(yjk_r) || !FINITE(yjk_i))
 363 
					|| NAN(yjk_r) || NAN(yjk_i))
359 
					goto next_ik;
 364 
					goto next_ik;
- 
 
 365 
#endif
360 
				sum_r += (xij_r * yjk_r - xij_i * yjk_i);
 366 
				sum_r += (xij_r * yjk_r - xij_i * yjk_i);
361 
				sum_i += (xij_r * yjk_i + xij_i * yjk_r);
 367 
				sum_i += (xij_r * yjk_i + xij_i * yjk_r);
362 
			}
 368 
			}
363 
			z[i+k*nrx].r = sum_r;
 369 
			z[i+k*nrx].r = sum_r;
364 
			z[i+k*nrx].i = sum_i;
 370 
			z[i+k*nrx].i = sum_i;
Line 377... Line 383...
377 
			z[i + k * ncx] = NA_REAL;
 383 
			z[i + k * ncx] = NA_REAL;
378 
			sum = 0.0;
 384 
			sum = 0.0;
379 
			for (j = 0; j < nrx; j++) {
 385 
			for (j = 0; j < nrx; j++) {
380 
				xji = x[j + i * nrx];
 386 
				xji = x[j + i * nrx];
381 
				yjk = y[j + k * nry];
 387 
				yjk = y[j + k * nry];
- 
 
 388 
#ifndef IEEE_754
382 
				if (!FINITE(xji) || !FINITE(yjk))
 389 
				if (NAN(xji) || NAN(yjk))
383 
					goto next_ik;
 390 
					goto next_ik;
- 
 
 391 
#endif
384 
				sum += xji * yjk;
 392 
				sum += xji * yjk;
385 
			}
 393 
			}
386 
			z[i + k * ncx] = sum;
 394 
			z[i + k * ncx] = sum;
387 
		next_ik:
 395 
		next_ik:
388 
			;
 396 
			;
Line 404... Line 412...
404 
			for (j = 0; j < nrx; j++) {
 412 
			for (j = 0; j < nrx; j++) {
405 
				xji_r = x[j + i * nrx].r;
 413 
				xji_r = x[j + i * nrx].r;
406 
				xji_i = x[j + i * nrx].i;
 414 
				xji_i = x[j + i * nrx].i;
407 
				yjk_r = y[j + k * nry].r;
 415 
				yjk_r = y[j + k * nry].r;
408 
				yjk_i = y[j + k * nry].i;
 416 
				yjk_i = y[j + k * nry].i;
- 
 
 417 
#ifndef IEEE_754
409 
				if (!FINITE(xji_r) || !FINITE(xji_i)
 418 
				if (NAN(xji_r) || NAN(xji_i)
410 
						|| !FINITE(yjk_r) || !FINITE(yjk_i))
 419 
					|| NAN(yjk_r) || NAN(yjk_i))
411 
					goto next_ik;
 420 
					goto next_ik;
- 
 
 421 
#endif
412 
				sum_r += (xji_r * yjk_r - xji_i * yjk_i);
 422 
				sum_r += (xji_r * yjk_r - xji_i * yjk_i);
413 
				sum_i += (xji_r * yjk_i + xji_i * yjk_r);
 423 
				sum_i += (xji_r * yjk_i + xji_i * yjk_r);
414 
			}
 424 
			}
415 
			z[i + k * ncx].r = sum_r;
 425 
			z[i + k * ncx].r = sum_r;
416 
			z[i + k * ncx].i = sum_i;
 426 
			z[i + k * ncx].i = sum_i;