The R Project SVN R-packages

Rev

Rev 4529 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 4529 Rev 4560
Line 343... Line 343...
343 
    SEXP val = PROTECT(allocVector(VECSXP, 3));
 343 
    SEXP val = PROTECT(allocVector(VECSXP, 3));
344 
 
 344 
 
345 
    if (dims[0] && dims[1]) {
 345 
    if (dims[0] && dims[1]) {
346 
	int m = dims[0], n = dims[1], mm = (m < n)?m:n,
 346 
	int m = dims[0], n = dims[1], mm = (m < n)?m:n,
347 
	    lwork = -1, info;
 347 
	    lwork = -1, info;
348 
	int *iwork = Calloc(8 * mm, int);
 - 
 
349 
	double tmp, *work;
 348 
	double tmp, *work;
350 
/* 	int bdspac = 3*m*m + 4*m, */
 349 
	int *iwork = Alloca(8 * mm, int);
351 
/* 	    wrkbl, maxwrk, minwrk, itmp, */
 - 
 
352 
/* 	    ione = 1, iminus1 = -1; */
 - 
 
353 
/* 	int i1, i2, i3; */
 350 
	R_CheckStack();
354 
 
 351 
 
355 
	SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));
 352 
	SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));
356 
	SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));
 353 
	SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));
357 
	SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));
 354 
	SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));
358 
	F77_CALL(dgesdd)("S", &m, &n, xx, &m,
 355 
	F77_CALL(dgesdd)("S", &m, &n, xx, &m,
359 
			 REAL(VECTOR_ELT(val, 0)),
 356 
			 REAL(VECTOR_ELT(val, 0)),
360 
			 REAL(VECTOR_ELT(val, 1)), &m,
 357 
			 REAL(VECTOR_ELT(val, 1)), &m,
361 
			 REAL(VECTOR_ELT(val, 2)), &mm,
 358 
			 REAL(VECTOR_ELT(val, 2)), &mm,
362 
			 &tmp, &lwork, iwork, &info);
 359 
			 &tmp, &lwork, iwork, &info);
363 
	lwork = (int) tmp;
 360 
	lwork = (int) tmp;
364 
/* 	F77_CALL(foo)(&i1, &i2, &i3); */
 - 
 
365 
/* 	wrkbl = 3*m+(m+n)*i1; */
 - 
 
366 
/* 	if (wrkbl < (itmp = 3*m + m*i2)) wrkbl = itmp; */
 - 
 
367 
/* 	if (wrkbl < (itmp = 3*m + m*i3)) wrkbl = itmp; */
 - 
 
368 
/* 	itmp = bdspac+3*m; */
 - 
 
369 
/* 	maxwrk = (wrkbl > itmp) ? wrkbl : itmp; */
 - 
 
370 
/* 	minwrk = 3*m + ((bdspac > n) ?  bdspac : n); */
 - 
 
371 
	work = Calloc(lwork, double);
 361 
	work = Alloca(lwork, double);
- 
 
 362 
	R_CheckStack();
372 
	F77_CALL(dgesdd)("S", &m, &n, xx, &m,
 363 
	F77_CALL(dgesdd)("S", &m, &n, xx, &m,
373 
			 REAL(VECTOR_ELT(val, 0)),
 364 
			 REAL(VECTOR_ELT(val, 0)),
374 
			 REAL(VECTOR_ELT(val, 1)), &m,
 365 
			 REAL(VECTOR_ELT(val, 1)), &m,
375 
			 REAL(VECTOR_ELT(val, 2)), &mm,
 366 
			 REAL(VECTOR_ELT(val, 2)), &mm,
376 
			 work, &lwork, iwork, &info);
 367 
			 work, &lwork, iwork, &info);
377 
	Free(iwork); Free(work);
 - 
 
- 
 
 368 
 
378 
    }
 369 
    }
379 
    UNPROTECT(1);
 370 
    UNPROTECT(1);
380 
    return val;
 371 
    return val;
381 
}
 372 
}
382 
 
 373 
 
Line 403... Line 394...
403 
{
 394 
{
404 
    SEXP val = PROTECT(duplicate(x));
 395 
    SEXP val = PROTECT(duplicate(x));
405 
    int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
 396 
    int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
406 
    int i, ilo, ilos, ihi, ihis, j, nc = Dims[1], sqpow;
 397 
    int i, ilo, ilos, ihi, ihis, j, nc = Dims[1], sqpow;
407 
    int ncp1 = Dims[1] + 1, ncsqr = nc * nc;
 398 
    int ncp1 = Dims[1] + 1, ncsqr = nc * nc;
408 
    int *pivot = Calloc(nc, int);
 399 
    int *pivot = Alloca(nc, int);
409 
    int *iperm = Calloc(nc, int);
 400 
    int *iperm = Alloca(nc, int);
410 
    double *dpp = Calloc(ncsqr, double), /* denominator power Pade' */
 401 
    double *dpp = Alloca(ncsqr, double), /* denominator power Pade' */
411 
	*npp = Calloc(ncsqr, double), /* numerator power Pade' */
 402 
	*npp = Alloca(ncsqr, double), /* numerator power Pade' */
412 
	*perm = Calloc(nc, double),
 403 
	*perm = Alloca(nc, double),
413 
	*scale = Calloc(nc, double),
 404 
	*scale = Alloca(nc, double),
414 
	*v = REAL(GET_SLOT(val, Matrix_xSym)),
 405 
	*v = REAL(GET_SLOT(val, Matrix_xSym)),
415 
	*work = Calloc(ncsqr, double), inf_norm, m1_j, /* (-1)^j */
 406 
	*work = Alloca(ncsqr, double), inf_norm, m1_j, /* (-1)^j */
416 
	one = 1., trshift, zero = 0.;
 407 
	one = 1., trshift, zero = 0.;
- 
 
 408 
    R_CheckStack();
417 
 
 409 
 
418 
    if (nc < 1 || Dims[0] != nc)
 410 
    if (nc < 1 || Dims[0] != nc)
419 
	error(_("Matrix exponential requires square, non-null matrix"));
 411 
	error(_("Matrix exponential requires square, non-null matrix"));
420 
 
 412 
 
421 
    /* FIXME: Add special treatment for nc == 1 */
 413 
    /* FIXME: Add special treatment for nc == 1 */
Line 516... Line 508...
516 
	double mult = exp(trshift);
 508 
	double mult = exp(trshift);
517 
	for (i = 0; i < ncsqr; i++) v[i] *= mult;
 509 
	for (i = 0; i < ncsqr; i++) v[i] *= mult;
518 
    }
 510 
    }
519 
 
 511 
 
520 
    /* Clean up */
 512 
    /* Clean up */
521 
    Free(dpp); Free(npp); Free(perm); Free(iperm); Free(pivot); Free(scale); Free(work);
 - 
 
522 
    UNPROTECT(1);
 513 
    UNPROTECT(1);
523 
    return val;
 514 
    return val;
524 
}
 515 
}
525 
 
 516 
 
526 
SEXP dgeMatrix_Schur(SEXP x, SEXP vectors)
 517 
SEXP dgeMatrix_Schur(SEXP x, SEXP vectors)
Line 541... Line 532...
541 
    F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,
 532 
    F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,
542 
		    (double *) NULL, (double *) NULL, (double *) NULL, dims,
 533 
		    (double *) NULL, (double *) NULL, (double *) NULL, dims,
543 
		    &tmp, &lwork, (int *) NULL, &info);
 534 
		    &tmp, &lwork, (int *) NULL, &info);
544 
    if (info) error(_("dgeMatrix_Schur: first call to dgees failed"));
 535 
    if (info) error(_("dgeMatrix_Schur: first call to dgees failed"));
545 
    lwork = (int) tmp;
 536 
    lwork = (int) tmp;
546 
    work = Calloc(lwork, double);
 537 
    work = Alloca(lwork, double);
- 
 
 538 
    R_CheckStack();
547 
    F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,
 539 
    F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,
548 
		    &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),
 540 
		    &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),
549 
		    REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,
 541 
		    REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,
550 
		    (int *) NULL, &info);
 542 
		    (int *) NULL, &info);
551 
    if (info) error(_("dgeMatrix_Schur: dgees returned code %d"), info);
 543 
    if (info) error(_("dgeMatrix_Schur: dgees returned code %d"), info);
552 
    Free(work);
 - 
 
553 
    UNPROTECT(1);
 544 
    UNPROTECT(1);
554 
    return val;
 545 
    return val;
555 
}
 546 
}
556 
 
 547 
 
557 
SEXP dgeMatrix_colsums(SEXP x, SEXP naRmP, SEXP cols, SEXP mean)
 548 
SEXP dgeMatrix_colsums(SEXP x, SEXP naRmP, SEXP cols, SEXP mean)
Line 580... Line 571...
580 
	    REAL(ans)[j] = sum;
 571 
	    REAL(ans)[j] = sum;
581 
	}
 572 
	}
582 
    } else {
 573 
    } else {
583 
	double *rans = REAL(ans), *ra = rans, *rx = xx, *Cnt = NULL, *c;
 574 
	double *rans = REAL(ans), *ra = rans, *rx = xx, *Cnt = NULL, *c;
584 
	cnt = p;
 575 
	cnt = p;
585 
	if (!keepNA && doMean) Cnt = Calloc(n, double);
 576 
	if (!keepNA && doMean) Cnt = Alloca(n, double);
- 
 
 577 
	R_CheckStack();
586 
	for (ra = rans, i = 0; i < n; i++) *ra++ = 0.0;
 578 
	for (ra = rans, i = 0; i < n; i++) *ra++ = 0.0;
587 
	for (j = 0; j < p; j++) {
 579 
	for (j = 0; j < p; j++) {
588 
	    ra = rans;
 580 
	    ra = rans;
589 
	    if (keepNA)
 581 
	    if (keepNA)
590 
		for (i = 0; i < n; i++) *ra++ += *rx++;
 582 
		for (i = 0; i < n; i++) *ra++ += *rx++;
Line 600... Line 592...
600 
		for (ra = rans, i = 0; i < n; i++)
 592 
		for (ra = rans, i = 0; i < n; i++)
601 
		    *ra++ /= p;
 593 
		    *ra++ /= p;
602 
	    else {
 594 
	    else {
603 
		for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)
 595 
		for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)
604 
		    if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
 596 
		    if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
605 
		Free(Cnt);
 - 
 
606 
	    }
 597 
	    }
607 
	}
 598 
	}
608 
    }
 599 
    }
609 
 
 600 
 
610 
    UNPROTECT(1);
 601 
    UNPROTECT(1);