The R Project SVN – R – Diff – /trunk/src/modules/lapack/Lapack.c



			 &n, &p, xvals, &n, REAL(s),
			 REAL(u), INTEGER(getAttrib(u, R_DimSymbol)),
			 REAL(v), INTEGER(getAttrib(v, R_DimSymbol)),
			 &tmp, &lwork, &info);
	if (info != 0)
	    error(_("error code %d from Lapack routine '%s'"), info, "dgesvd");
	lwork = (int) tmp;
 
	work = (double *) R_alloc(lwork, sizeof(double));
	F77_CALL(dgesvd)(CHAR(STRING_ELT(jobu, 0)), CHAR(STRING_ELT(jobv, 0)),
			 &n, &p, xvals, &n, REAL(s),
			 REAL(u), INTEGER(getAttrib(u, R_DimSymbol)),
			 REAL(v), INTEGER(getAttrib(v, R_DimSymbol)),
			 work, &lwork, &info);
	if (info != 0)
	    error(_("error code %d from Lapack routine '%s'"), info, "dgesvd");
    } else {
	int ldu = INTEGER(getAttrib(u, R_DimSymbol))[0],
	    ldvt = INTEGER(getAttrib(v, R_DimSymbol))[0];
	int *iwork= (int *) R_alloc(8*(n<p ? n : p), sizeof(int));
 

			 &n, &p, xvals, &n, REAL(s),
			 REAL(u), &ldu,
			 REAL(v), &ldvt,
			 &tmp, &lwork, iwork, &info);
	if (info != 0)
	    error(_("error code %d from Lapack routine '%s'"), info, "dgesdd");
	lwork = (int) tmp;
	work = (double *) R_alloc(lwork, sizeof(double));
	F77_CALL(dgesdd)(CHAR(STRING_ELT(jobu, 0)),
			 &n, &p, xvals, &n, REAL(s),
			 REAL(u), &ldu,
			 REAL(v), &ldvt,
			 work, &lwork, iwork, &info);
	if (info != 0)
	    error(_("error code %d from Lapack routine '%s'"), info, "dgesdd");
    }
 
    val = PROTECT(allocVector(VECSXP, 3));
    nm = PROTECT(allocVector(STRSXP, 3));
    SET_STRING_ELT(nm, 0, mkChar("d"));

	F77_CALL(dsyev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, &info);
#else
        F77_CALL(rsyev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, &info);
#endif
	if (info != 0)
	    error(_("error code %d from Lapack routine '%s'"), info, "dsyev");
	lwork = (int) tmp;
	if (lwork < 3*n-1) lwork = 3*n-1;  /* Sanity check */
	work = (double *) R_alloc(lwork, sizeof(double));
#ifdef HAVE_LAPACK
	F77_CALL(dsyev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, &info);
#else
        F77_CALL(rsyev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, &info);
#endif
	if (info != 0)
	    error(_("error code %d from Lapack routine '%s'"), info, "dsyev");
    } else {
	int liwork, *iwork, itmp, m;
	double vl = 0.0, vu = 0.0, abstol = 0.0; 
	/* valgrind seems to think vu should be set, but it is documented 
	   not to be used if range='a' */

                         &vl, &vu, &il, &iu, &abstol, &m, rvalues,
                         REAL(z), &n, isuppz,
                         &tmp, &lwork, &itmp, &liwork, &info);
#endif
	if (info != 0)
	    error(_("error code %d from Lapack routine '%s'"), info, "dsyevr");
	lwork = (int) tmp;
	liwork = itmp;
 
	work = (double *) R_alloc(lwork, sizeof(double));
	iwork = (int *) R_alloc(liwork, sizeof(int));

                         &vl, &vu, &il, &iu, &abstol, &m, rvalues,
                         REAL(z), &n, isuppz,
                         work, &lwork, iwork, &liwork, &info);
#endif
	if (info != 0)
	    error(_("error code %d from Lapack routine '%s'"), info, "dsyevr");
    }
 
    if (!ov) {
	ret = PROTECT(allocVector(VECSXP, 2));
	nm = PROTECT(allocVector(STRSXP, 2));

#else
    F77_CALL(rgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI,
		    left, &n, right, &n, &tmp, &lwork, &info);
#endif
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dgeev");
    lwork = (int) tmp;
    work = (double *) R_alloc(lwork, sizeof(double));
#ifdef HAVE_LAPACK
    F77_CALL(dgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI,
		    left, &n, right, &n, work, &lwork, &info);
#else
    F77_CALL(rgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI,
		    left, &n, right, &n, work, &lwork, &info);
#endif
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dgeev");
 
    complexValues = FALSE;
    for (i = 0; i < n; i++)
	if (wI[i] != 0.0) { complexValues = TRUE; break; }
    ret = PROTECT(allocVector(VECSXP, 2));

    tau = PROTECT(allocVector(CPLXSXP, m < n ? m : n));
    lwork = -1;
    F77_CALL(zgeqp3)(&m, &n, COMPLEX(A), &m, INTEGER(jpvt), COMPLEX(tau),
		     &tmp, &lwork, rwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zgeqp3");
    lwork = (int) tmp.r;
    work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));
    F77_CALL(zgeqp3)(&m, &n, COMPLEX(A), &m, INTEGER(jpvt), COMPLEX(tau),
		     work, &lwork, rwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zgeqp3");
    val = PROTECT(allocVector(VECSXP, 4));
    nm = PROTECT(allocVector(STRSXP, 4));
    rank = PROTECT(allocVector(INTSXP, 1));
    INTEGER(rank)[0] = m < n ? m : n;
    SET_STRING_ELT(nm, 0, mkChar("qr"));

    lwork = -1;
    F77_CALL(zunmqr)("L", "C", &n, &nrhs, &k,
		     COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,
		     &tmp, &lwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zunmqr");
    lwork = (int) tmp.r;
    work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));
    F77_CALL(zunmqr)("L", "C", &n, &nrhs, &k,
		     COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,
		     work, &lwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zunmqr");
    F77_CALL(ztrtrs)("U", "N", "N", &k, &nrhs,
		     COMPLEX(qr), &n, COMPLEX(B), &n, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "ztrtrs");
    UNPROTECT(1);
    return B;
#else
    error(_("Fortran complex functions are not available on this platform"));
    return R_NilValue; /* -Wall */

    lwork = -1;
    F77_CALL(zunmqr)("L", tr ? "C" : "N", &n, &nrhs, &k,
		     COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,
		     &tmp, &lwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zunmqr");
    lwork = (int) tmp.r;
    work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));
    F77_CALL(zunmqr)("L", tr ? "C" : "N", &n, &nrhs, &k,
		     COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,
		     work, &lwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zunmqr");
    UNPROTECT(1);
    return B;
#else
    error(_("Fortran complex functions are not available on this platform"));
    return R_NilValue; /* -Wall */

		     &n, &p, COMPLEX(x), &n, REAL(s),
		     COMPLEX(u), INTEGER(getAttrib(u, R_DimSymbol)),
		     COMPLEX(v), INTEGER(getAttrib(v, R_DimSymbol)),
		     &tmp, &lwork, rwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zgesvd");
    lwork = (int) tmp.r;
    work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));
    F77_CALL(zgesvd)(CHAR(STRING_ELT(jobu, 0)), CHAR(STRING_ELT(jobv, 0)),
		     &n, &p, COMPLEX(x), &n, REAL(s),
		     COMPLEX(u), INTEGER(getAttrib(u, R_DimSymbol)),
		     COMPLEX(v), INTEGER(getAttrib(v, R_DimSymbol)),
		     work, &lwork, rwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zgesvd");
    val = PROTECT(allocVector(VECSXP, 3));
    nm = PROTECT(allocVector(STRSXP, 3));
    SET_STRING_ELT(nm, 0, mkChar("d"));
    SET_STRING_ELT(nm, 1, mkChar("u"));
    SET_STRING_ELT(nm, 2, mkChar("vt"));

    /* ask for optimal size of work array */
    lwork = -1;
    F77_CALL(zheev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, rwork,
		    &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zheev");
    lwork = (int) tmp.r;
    work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));
    F77_CALL(zheev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, rwork,
		    &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zheev");
    if (!ov) {
	ret = PROTECT(allocVector(VECSXP, 2));
	nm = PROTECT(allocVector(STRSXP, 2));
	SET_STRING_ELT(nm, 1, mkChar("vectors"));
	SET_VECTOR_ELT(ret, 1, x);

    /* ask for optimal size of work array */
    lwork = -1;
    F77_CALL(zgeev)(jobVL, jobVR, &n, xvals, &n, COMPLEX(values),
		    left, &n, right, &n, &tmp, &lwork, rwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zgeev");
    lwork = (int) tmp.r;
    work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));
    F77_CALL(zgeev)(jobVL, jobVR, &n, xvals, &n, COMPLEX(values),
		    left, &n, right, &n, work, &lwork, rwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "zgeev");
 
    if(!ov){
	ret = PROTECT(allocVector(VECSXP, 2));
	nm = PROTECT(allocVector(STRSXP, 2));
	SET_STRING_ELT(nm, 1, mkChar("vectors"));

    tau = PROTECT(allocVector(REALSXP, m < n ? m : n));
    lwork = -1;
    F77_CALL(dgeqp3)(&m, &n, REAL(A), &m, INTEGER(jpvt), REAL(tau),
		     &tmp, &lwork, &info);
    if (info < 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dgeqp3");
    lwork = (int) tmp;
    work = (double *) R_alloc(lwork, sizeof(double));
    F77_CALL(dgeqp3)(&m, &n, REAL(A), &m, INTEGER(jpvt), REAL(tau),
		     work, &lwork, &info);
    if (info < 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dgeqp3");
    val = PROTECT(allocVector(VECSXP, 4));
    nm = PROTECT(allocVector(STRSXP, 4));
    rank = PROTECT(allocVector(INTSXP, 1));
    INTEGER(rank)[0] = m < n ? m : n;
    SET_STRING_ELT(nm, 0, mkChar("qr"));

    lwork = -1;
    F77_CALL(dormqr)("L", "T", &n, &nrhs, &k,
		     REAL(qr), &n, REAL(tau), REAL(B), &n,
		     &tmp, &lwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dormqr");
    lwork = (int) tmp;
    work = (double *) R_alloc(lwork, sizeof(double));
    F77_CALL(dormqr)("L", "T", &n, &nrhs, &k,
		     REAL(qr), &n, REAL(tau), REAL(B), &n,
		     work, &lwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dormqr");
    F77_CALL(dtrtrs)("U", "N", "N", &k, &nrhs,
		     REAL(qr), &n, REAL(B), &n, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dtrtrs");
    UNPROTECT(1);
    return B;
}
 
static SEXP modqr_qy_real(SEXP Q, SEXP Bin, SEXP trans)

    lwork = -1;
    F77_CALL(dormqr)("L", tr ? "T" : "N", &n, &nrhs, &k,
		     REAL(qr), &n, REAL(tau), REAL(B), &n,
		     &tmp, &lwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dormqr");
    lwork = (int) tmp;
    work = (double *) R_alloc(lwork, sizeof(double));
    F77_CALL(dormqr)("L", tr ? "T" : "N", &n, &nrhs, &k,
		     REAL(qr), &n, REAL(tau), REAL(B), &n,
		     work, &lwork, &info);
    if (info != 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dormqr");
    UNPROTECT(1);
    return B;
}
 
static SEXP moddet_ge_real(SEXP Ain, SEXP logarithm)

	error(_("'A' must be a square matrix"));
    jpvt = (int *) R_alloc(n, sizeof(int));
    F77_CALL(dgetrf)(&n, &n, REAL(A), &n, jpvt, &info);
    sign = 1;
    if (info < 0)
	error(_("error code %d from Lapack routine '%s'"), info, "dgetrf");
    else if (info > 0) { /* Singular matrix:  U[i,i] (i := info) is 0 */
	/*warning("Lapack dgetrf(): singular matrix: U[%d,%d]=0", info,info);*/
	modulus = (useLog ? R_NegInf : 0.);
    }
    else {

Rev 33282	Rev 33387
Line 65...	Line 65...
65	`&n, &p, xvals, &n, REAL(s),`	65	`&n, &p, xvals, &n, REAL(s),`
66	`REAL(u), INTEGER(getAttrib(u, R_DimSymbol)),`	66	`REAL(u), INTEGER(getAttrib(u, R_DimSymbol)),`
67	`REAL(v), INTEGER(getAttrib(v, R_DimSymbol)),`	67	`REAL(v), INTEGER(getAttrib(v, R_DimSymbol)),`
68	`&tmp, &lwork, &info);`	68	`&tmp, &lwork, &info);`
69	`if (info != 0)`	69	`if (info != 0)`
70	`error(_("error code %d from Lapack routine %s"), info, "dgesvd");`	70	`error(_("error code %d from Lapack routine '%s'"), info, "dgesvd");`
71	`lwork = (int) tmp;`	71	`lwork = (int) tmp;`
72		72
73	`work = (double *) R_alloc(lwork, sizeof(double));`	73	`work = (double *) R_alloc(lwork, sizeof(double));`
74	`F77_CALL(dgesvd)(CHAR(STRING_ELT(jobu, 0)), CHAR(STRING_ELT(jobv, 0)),`	74	`F77_CALL(dgesvd)(CHAR(STRING_ELT(jobu, 0)), CHAR(STRING_ELT(jobv, 0)),`
75	`&n, &p, xvals, &n, REAL(s),`	75	`&n, &p, xvals, &n, REAL(s),`
76	`REAL(u), INTEGER(getAttrib(u, R_DimSymbol)),`	76	`REAL(u), INTEGER(getAttrib(u, R_DimSymbol)),`
77	`REAL(v), INTEGER(getAttrib(v, R_DimSymbol)),`	77	`REAL(v), INTEGER(getAttrib(v, R_DimSymbol)),`
78	`work, &lwork, &info);`	78	`work, &lwork, &info);`
79	`if (info != 0)`	79	`if (info != 0)`
80	`error(_("error code %d from Lapack routine %s"), info, "dgesvd");`	80	`error(_("error code %d from Lapack routine '%s'"), info, "dgesvd");`
81	`} else {`	81	`} else {`
82	`int ldu = INTEGER(getAttrib(u, R_DimSymbol))[0],`	82	`int ldu = INTEGER(getAttrib(u, R_DimSymbol))[0],`
83	`ldvt = INTEGER(getAttrib(v, R_DimSymbol))[0];`	83	`ldvt = INTEGER(getAttrib(v, R_DimSymbol))[0];`
84	`int iwork= (int ) R_alloc(8*(n<p ? n : p), sizeof(int));`	84	`int iwork= (int ) R_alloc(8*(n<p ? n : p), sizeof(int));`
85		85
Line 89...	Line 89...
89	`&n, &p, xvals, &n, REAL(s),`	89	`&n, &p, xvals, &n, REAL(s),`
90	`REAL(u), &ldu,`	90	`REAL(u), &ldu,`
91	`REAL(v), &ldvt,`	91	`REAL(v), &ldvt,`
92	`&tmp, &lwork, iwork, &info);`	92	`&tmp, &lwork, iwork, &info);`
93	`if (info != 0)`	93	`if (info != 0)`
94	`error(_("error code %d from Lapack routine %s"), info, "dgesdd");`	94	`error(_("error code %d from Lapack routine '%s'"), info, "dgesdd");`
95	`lwork = (int) tmp;`	95	`lwork = (int) tmp;`
96	`work = (double *) R_alloc(lwork, sizeof(double));`	96	`work = (double *) R_alloc(lwork, sizeof(double));`
97	`F77_CALL(dgesdd)(CHAR(STRING_ELT(jobu, 0)),`	97	`F77_CALL(dgesdd)(CHAR(STRING_ELT(jobu, 0)),`
98	`&n, &p, xvals, &n, REAL(s),`	98	`&n, &p, xvals, &n, REAL(s),`
99	`REAL(u), &ldu,`	99	`REAL(u), &ldu,`
100	`REAL(v), &ldvt,`	100	`REAL(v), &ldvt,`
101	`work, &lwork, iwork, &info);`	101	`work, &lwork, iwork, &info);`
102	`if (info != 0)`	102	`if (info != 0)`
103	`error(_("error code %d from Lapack routine %s"), info, "dgesdd");`	103	`error(_("error code %d from Lapack routine '%s'"), info, "dgesdd");`
104	`}`	104	`}`
105		105
106	`val = PROTECT(allocVector(VECSXP, 3));`	106	`val = PROTECT(allocVector(VECSXP, 3));`
107	`nm = PROTECT(allocVector(STRSXP, 3));`	107	`nm = PROTECT(allocVector(STRSXP, 3));`
108	`SET_STRING_ELT(nm, 0, mkChar("d"));`	108	`SET_STRING_ELT(nm, 0, mkChar("d"));`
Line 153...	Line 153...
153	`F77_CALL(dsyev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, &info);`	153	`F77_CALL(dsyev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, &info);`
154	`#else`	154	`#else`
155	`F77_CALL(rsyev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, &info);`	155	`F77_CALL(rsyev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, &info);`
156	`#endif`	156	`#endif`
157	`if (info != 0)`	157	`if (info != 0)`
158	`error(_("error code %d from Lapack routine %s"), info, "dsyev");`	158	`error(_("error code %d from Lapack routine '%s'"), info, "dsyev");`
159	`lwork = (int) tmp;`	159	`lwork = (int) tmp;`
160	`if (lwork < 3n-1) lwork = 3n-1; /* Sanity check */`	160	`if (lwork < 3n-1) lwork = 3n-1; /* Sanity check */`
161	`work = (double *) R_alloc(lwork, sizeof(double));`	161	`work = (double *) R_alloc(lwork, sizeof(double));`
162	`#ifdef HAVE_LAPACK`	162	`#ifdef HAVE_LAPACK`
163	`F77_CALL(dsyev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, &info);`	163	`F77_CALL(dsyev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, &info);`
164	`#else`	164	`#else`
165	`F77_CALL(rsyev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, &info);`	165	`F77_CALL(rsyev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, &info);`
166	`#endif`	166	`#endif`
167	`if (info != 0)`	167	`if (info != 0)`
168	`error(_("error code %d from Lapack routine %s"), info, "dsyev");`	168	`error(_("error code %d from Lapack routine '%s'"), info, "dsyev");`
169	`} else {`	169	`} else {`
170	`int liwork, *iwork, itmp, m;`	170	`int liwork, *iwork, itmp, m;`
171	`double vl = 0.0, vu = 0.0, abstol = 0.0;`	171	`double vl = 0.0, vu = 0.0, abstol = 0.0;`
172	`/* valgrind seems to think vu should be set, but it is documented`	172	`/* valgrind seems to think vu should be set, but it is documented`
173	`not to be used if range='a' */`	173	`not to be used if range='a' */`
Line 188...	Line 188...
188	`&vl, &vu, &il, &iu, &abstol, &m, rvalues,`	188	`&vl, &vu, &il, &iu, &abstol, &m, rvalues,`
189	`REAL(z), &n, isuppz,`	189	`REAL(z), &n, isuppz,`
190	`&tmp, &lwork, &itmp, &liwork, &info);`	190	`&tmp, &lwork, &itmp, &liwork, &info);`
191	`#endif`	191	`#endif`
192	`if (info != 0)`	192	`if (info != 0)`
193	`error(_("error code %d from Lapack routine %s"), info, "dsyevr");`	193	`error(_("error code %d from Lapack routine '%s'"), info, "dsyevr");`
194	`lwork = (int) tmp;`	194	`lwork = (int) tmp;`
195	`liwork = itmp;`	195	`liwork = itmp;`
196		196
197	`work = (double *) R_alloc(lwork, sizeof(double));`	197	`work = (double *) R_alloc(lwork, sizeof(double));`
198	`iwork = (int *) R_alloc(liwork, sizeof(int));`	198	`iwork = (int *) R_alloc(liwork, sizeof(int));`
Line 206...	Line 206...
206	`&vl, &vu, &il, &iu, &abstol, &m, rvalues,`	206	`&vl, &vu, &il, &iu, &abstol, &m, rvalues,`
207	`REAL(z), &n, isuppz,`	207	`REAL(z), &n, isuppz,`
208	`work, &lwork, iwork, &liwork, &info);`	208	`work, &lwork, iwork, &liwork, &info);`
209	`#endif`	209	`#endif`
210	`if (info != 0)`	210	`if (info != 0)`
211	`error(_("error code %d from Lapack routine %s"), info, "dsyevr");`	211	`error(_("error code %d from Lapack routine '%s'"), info, "dsyevr");`
212	`}`	212	`}`
213		213
214	`if (!ov) {`	214	`if (!ov) {`
215	`ret = PROTECT(allocVector(VECSXP, 2));`	215	`ret = PROTECT(allocVector(VECSXP, 2));`
216	`nm = PROTECT(allocVector(STRSXP, 2));`	216	`nm = PROTECT(allocVector(STRSXP, 2));`
Line 291...	Line 291...
291	`#else`	291	`#else`
292	`F77_CALL(rgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI,`	292	`F77_CALL(rgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI,`
293	`left, &n, right, &n, &tmp, &lwork, &info);`	293	`left, &n, right, &n, &tmp, &lwork, &info);`
294	`#endif`	294	`#endif`
295	`if (info != 0)`	295	`if (info != 0)`
296	`error(_("error code %d from Lapack routine %s"), info, "dgeev");`	296	`error(_("error code %d from Lapack routine '%s'"), info, "dgeev");`
297	`lwork = (int) tmp;`	297	`lwork = (int) tmp;`
298	`work = (double *) R_alloc(lwork, sizeof(double));`	298	`work = (double *) R_alloc(lwork, sizeof(double));`
299	`#ifdef HAVE_LAPACK`	299	`#ifdef HAVE_LAPACK`
300	`F77_CALL(dgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI,`	300	`F77_CALL(dgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI,`
301	`left, &n, right, &n, work, &lwork, &info);`	301	`left, &n, right, &n, work, &lwork, &info);`
302	`#else`	302	`#else`
303	`F77_CALL(rgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI,`	303	`F77_CALL(rgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI,`
304	`left, &n, right, &n, work, &lwork, &info);`	304	`left, &n, right, &n, work, &lwork, &info);`
305	`#endif`	305	`#endif`
306	`if (info != 0)`	306	`if (info != 0)`
307	`error(_("error code %d from Lapack routine %s"), info, "dgeev");`	307	`error(_("error code %d from Lapack routine '%s'"), info, "dgeev");`
308		308
309	`complexValues = FALSE;`	309	`complexValues = FALSE;`
310	`for (i = 0; i < n; i++)`	310	`for (i = 0; i < n; i++)`
311	`if (wI[i] != 0.0) { complexValues = TRUE; break; }`	311	`if (wI[i] != 0.0) { complexValues = TRUE; break; }`
312	`ret = PROTECT(allocVector(VECSXP, 2));`	312	`ret = PROTECT(allocVector(VECSXP, 2));`
Line 406...	Line 406...
406	`tau = PROTECT(allocVector(CPLXSXP, m < n ? m : n));`	406	`tau = PROTECT(allocVector(CPLXSXP, m < n ? m : n));`
407	`lwork = -1;`	407	`lwork = -1;`
408	`F77_CALL(zgeqp3)(&m, &n, COMPLEX(A), &m, INTEGER(jpvt), COMPLEX(tau),`	408	`F77_CALL(zgeqp3)(&m, &n, COMPLEX(A), &m, INTEGER(jpvt), COMPLEX(tau),`
409	`&tmp, &lwork, rwork, &info);`	409	`&tmp, &lwork, rwork, &info);`
410	`if (info != 0)`	410	`if (info != 0)`
411	`error(_("error code %d from Lapack routine %s"), info, "zgeqp3");`	411	`error(_("error code %d from Lapack routine '%s'"), info, "zgeqp3");`
412	`lwork = (int) tmp.r;`	412	`lwork = (int) tmp.r;`
413	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`	413	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`
414	`F77_CALL(zgeqp3)(&m, &n, COMPLEX(A), &m, INTEGER(jpvt), COMPLEX(tau),`	414	`F77_CALL(zgeqp3)(&m, &n, COMPLEX(A), &m, INTEGER(jpvt), COMPLEX(tau),`
415	`work, &lwork, rwork, &info);`	415	`work, &lwork, rwork, &info);`
416	`if (info != 0)`	416	`if (info != 0)`
417	`error(_("error code %d from Lapack routine %s"), info, "zgeqp3");`	417	`error(_("error code %d from Lapack routine '%s'"), info, "zgeqp3");`
418	`val = PROTECT(allocVector(VECSXP, 4));`	418	`val = PROTECT(allocVector(VECSXP, 4));`
419	`nm = PROTECT(allocVector(STRSXP, 4));`	419	`nm = PROTECT(allocVector(STRSXP, 4));`
420	`rank = PROTECT(allocVector(INTSXP, 1));`	420	`rank = PROTECT(allocVector(INTSXP, 1));`
421	`INTEGER(rank)[0] = m < n ? m : n;`	421	`INTEGER(rank)[0] = m < n ? m : n;`
422	`SET_STRING_ELT(nm, 0, mkChar("qr"));`	422	`SET_STRING_ELT(nm, 0, mkChar("qr"));`
Line 457...	Line 457...
457	`lwork = -1;`	457	`lwork = -1;`
458	`F77_CALL(zunmqr)("L", "C", &n, &nrhs, &k,`	458	`F77_CALL(zunmqr)("L", "C", &n, &nrhs, &k,`
459	`COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,`	459	`COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,`
460	`&tmp, &lwork, &info);`	460	`&tmp, &lwork, &info);`
461	`if (info != 0)`	461	`if (info != 0)`
462	`error(_("error code %d from Lapack routine %s"), info, "zunmqr");`	462	`error(_("error code %d from Lapack routine '%s'"), info, "zunmqr");`
463	`lwork = (int) tmp.r;`	463	`lwork = (int) tmp.r;`
464	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`	464	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`
465	`F77_CALL(zunmqr)("L", "C", &n, &nrhs, &k,`	465	`F77_CALL(zunmqr)("L", "C", &n, &nrhs, &k,`
466	`COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,`	466	`COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,`
467	`work, &lwork, &info);`	467	`work, &lwork, &info);`
468	`if (info != 0)`	468	`if (info != 0)`
469	`error(_("error code %d from Lapack routine %s"), info, "zunmqr");`	469	`error(_("error code %d from Lapack routine '%s'"), info, "zunmqr");`
470	`F77_CALL(ztrtrs)("U", "N", "N", &k, &nrhs,`	470	`F77_CALL(ztrtrs)("U", "N", "N", &k, &nrhs,`
471	`COMPLEX(qr), &n, COMPLEX(B), &n, &info);`	471	`COMPLEX(qr), &n, COMPLEX(B), &n, &info);`
472	`if (info != 0)`	472	`if (info != 0)`
473	`error(_("error code %d from Lapack routine %s"), info, "ztrtrs");`	473	`error(_("error code %d from Lapack routine '%s'"), info, "ztrtrs");`
474	`UNPROTECT(1);`	474	`UNPROTECT(1);`
475	`return B;`	475	`return B;`
476	`#else`	476	`#else`
477	`error(_("Fortran complex functions are not available on this platform"));`	477	`error(_("Fortran complex functions are not available on this platform"));`
478	`return R_NilValue; /* -Wall */`	478	`return R_NilValue; /* -Wall */`
Line 502...	Line 502...
502	`lwork = -1;`	502	`lwork = -1;`
503	`F77_CALL(zunmqr)("L", tr ? "C" : "N", &n, &nrhs, &k,`	503	`F77_CALL(zunmqr)("L", tr ? "C" : "N", &n, &nrhs, &k,`
504	`COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,`	504	`COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,`
505	`&tmp, &lwork, &info);`	505	`&tmp, &lwork, &info);`
506	`if (info != 0)`	506	`if (info != 0)`
507	`error(_("error code %d from Lapack routine %s"), info, "zunmqr");`	507	`error(_("error code %d from Lapack routine '%s'"), info, "zunmqr");`
508	`lwork = (int) tmp.r;`	508	`lwork = (int) tmp.r;`
509	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`	509	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`
510	`F77_CALL(zunmqr)("L", tr ? "C" : "N", &n, &nrhs, &k,`	510	`F77_CALL(zunmqr)("L", tr ? "C" : "N", &n, &nrhs, &k,`
511	`COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,`	511	`COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n,`
512	`work, &lwork, &info);`	512	`work, &lwork, &info);`
513	`if (info != 0)`	513	`if (info != 0)`
514	`error(_("error code %d from Lapack routine %s"), info, "zunmqr");`	514	`error(_("error code %d from Lapack routine '%s'"), info, "zunmqr");`
515	`UNPROTECT(1);`	515	`UNPROTECT(1);`
516	`return B;`	516	`return B;`
517	`#else`	517	`#else`
518	`error(_("Fortran complex functions are not available on this platform"));`	518	`error(_("Fortran complex functions are not available on this platform"));`
519	`return R_NilValue; /* -Wall */`	519	`return R_NilValue; /* -Wall */`
Line 540...	Line 540...
540	`&n, &p, COMPLEX(x), &n, REAL(s),`	540	`&n, &p, COMPLEX(x), &n, REAL(s),`
541	`COMPLEX(u), INTEGER(getAttrib(u, R_DimSymbol)),`	541	`COMPLEX(u), INTEGER(getAttrib(u, R_DimSymbol)),`
542	`COMPLEX(v), INTEGER(getAttrib(v, R_DimSymbol)),`	542	`COMPLEX(v), INTEGER(getAttrib(v, R_DimSymbol)),`
543	`&tmp, &lwork, rwork, &info);`	543	`&tmp, &lwork, rwork, &info);`
544	`if (info != 0)`	544	`if (info != 0)`
545	`error(_("error code %d from Lapack routine %s"), info, "zgesvd");`	545	`error(_("error code %d from Lapack routine '%s'"), info, "zgesvd");`
546	`lwork = (int) tmp.r;`	546	`lwork = (int) tmp.r;`
547	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`	547	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`
548	`F77_CALL(zgesvd)(CHAR(STRING_ELT(jobu, 0)), CHAR(STRING_ELT(jobv, 0)),`	548	`F77_CALL(zgesvd)(CHAR(STRING_ELT(jobu, 0)), CHAR(STRING_ELT(jobv, 0)),`
549	`&n, &p, COMPLEX(x), &n, REAL(s),`	549	`&n, &p, COMPLEX(x), &n, REAL(s),`
550	`COMPLEX(u), INTEGER(getAttrib(u, R_DimSymbol)),`	550	`COMPLEX(u), INTEGER(getAttrib(u, R_DimSymbol)),`
551	`COMPLEX(v), INTEGER(getAttrib(v, R_DimSymbol)),`	551	`COMPLEX(v), INTEGER(getAttrib(v, R_DimSymbol)),`
552	`work, &lwork, rwork, &info);`	552	`work, &lwork, rwork, &info);`
553	`if (info != 0)`	553	`if (info != 0)`
554	`error(_("error code %d from Lapack routine %s"), info, "zgesvd");`	554	`error(_("error code %d from Lapack routine '%s'"), info, "zgesvd");`
555	`val = PROTECT(allocVector(VECSXP, 3));`	555	`val = PROTECT(allocVector(VECSXP, 3));`
556	`nm = PROTECT(allocVector(STRSXP, 3));`	556	`nm = PROTECT(allocVector(STRSXP, 3));`
557	`SET_STRING_ELT(nm, 0, mkChar("d"));`	557	`SET_STRING_ELT(nm, 0, mkChar("d"));`
558	`SET_STRING_ELT(nm, 1, mkChar("u"));`	558	`SET_STRING_ELT(nm, 1, mkChar("u"));`
559	`SET_STRING_ELT(nm, 2, mkChar("vt"));`	559	`SET_STRING_ELT(nm, 2, mkChar("vt"));`
Line 595...	Line 595...
595	`/* ask for optimal size of work array */`	595	`/* ask for optimal size of work array */`
596	`lwork = -1;`	596	`lwork = -1;`
597	`F77_CALL(zheev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, rwork,`	597	`F77_CALL(zheev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, rwork,`
598	`&info);`	598	`&info);`
599	`if (info != 0)`	599	`if (info != 0)`
600	`error(_("error code %d from Lapack routine %s"), info, "zheev");`	600	`error(_("error code %d from Lapack routine '%s'"), info, "zheev");`
601	`lwork = (int) tmp.r;`	601	`lwork = (int) tmp.r;`
602	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`	602	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`
603	`F77_CALL(zheev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, rwork,`	603	`F77_CALL(zheev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, rwork,`
604	`&info);`	604	`&info);`
605	`if (info != 0)`	605	`if (info != 0)`
606	`error(_("error code %d from Lapack routine %s"), info, "zheev");`	606	`error(_("error code %d from Lapack routine '%s'"), info, "zheev");`
607	`if (!ov) {`	607	`if (!ov) {`
608	`ret = PROTECT(allocVector(VECSXP, 2));`	608	`ret = PROTECT(allocVector(VECSXP, 2));`
609	`nm = PROTECT(allocVector(STRSXP, 2));`	609	`nm = PROTECT(allocVector(STRSXP, 2));`
610	`SET_STRING_ELT(nm, 1, mkChar("vectors"));`	610	`SET_STRING_ELT(nm, 1, mkChar("vectors"));`
611	`SET_VECTOR_ELT(ret, 1, x);`	611	`SET_VECTOR_ELT(ret, 1, x);`
Line 656...	Line 656...
656	`/* ask for optimal size of work array */`	656	`/* ask for optimal size of work array */`
657	`lwork = -1;`	657	`lwork = -1;`
658	`F77_CALL(zgeev)(jobVL, jobVR, &n, xvals, &n, COMPLEX(values),`	658	`F77_CALL(zgeev)(jobVL, jobVR, &n, xvals, &n, COMPLEX(values),`
659	`left, &n, right, &n, &tmp, &lwork, rwork, &info);`	659	`left, &n, right, &n, &tmp, &lwork, rwork, &info);`
660	`if (info != 0)`	660	`if (info != 0)`
661	`error(_("error code %d from Lapack routine %s"), info, "zgeev");`	661	`error(_("error code %d from Lapack routine '%s'"), info, "zgeev");`
662	`lwork = (int) tmp.r;`	662	`lwork = (int) tmp.r;`
663	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`	663	`work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex));`
664	`F77_CALL(zgeev)(jobVL, jobVR, &n, xvals, &n, COMPLEX(values),`	664	`F77_CALL(zgeev)(jobVL, jobVR, &n, xvals, &n, COMPLEX(values),`
665	`left, &n, right, &n, work, &lwork, rwork, &info);`	665	`left, &n, right, &n, work, &lwork, rwork, &info);`
666	`if (info != 0)`	666	`if (info != 0)`
667	`error(_("error code %d from Lapack routine %s"), info, "zgeev");`	667	`error(_("error code %d from Lapack routine '%s'"), info, "zgeev");`
668		668
669	`if(!ov){`	669	`if(!ov){`
670	`ret = PROTECT(allocVector(VECSXP, 2));`	670	`ret = PROTECT(allocVector(VECSXP, 2));`
671	`nm = PROTECT(allocVector(STRSXP, 2));`	671	`nm = PROTECT(allocVector(STRSXP, 2));`
672	`SET_STRING_ELT(nm, 1, mkChar("vectors"));`	672	`SET_STRING_ELT(nm, 1, mkChar("vectors"));`
Line 822...	Line 822...
822	`tau = PROTECT(allocVector(REALSXP, m < n ? m : n));`	822	`tau = PROTECT(allocVector(REALSXP, m < n ? m : n));`
823	`lwork = -1;`	823	`lwork = -1;`
824	`F77_CALL(dgeqp3)(&m, &n, REAL(A), &m, INTEGER(jpvt), REAL(tau),`	824	`F77_CALL(dgeqp3)(&m, &n, REAL(A), &m, INTEGER(jpvt), REAL(tau),`
825	`&tmp, &lwork, &info);`	825	`&tmp, &lwork, &info);`
826	`if (info < 0)`	826	`if (info < 0)`
827	`error(_("error code %d from Lapack routine %s"), info, "dgeqp3");`	827	`error(_("error code %d from Lapack routine '%s'"), info, "dgeqp3");`
828	`lwork = (int) tmp;`	828	`lwork = (int) tmp;`
829	`work = (double *) R_alloc(lwork, sizeof(double));`	829	`work = (double *) R_alloc(lwork, sizeof(double));`
830	`F77_CALL(dgeqp3)(&m, &n, REAL(A), &m, INTEGER(jpvt), REAL(tau),`	830	`F77_CALL(dgeqp3)(&m, &n, REAL(A), &m, INTEGER(jpvt), REAL(tau),`
831	`work, &lwork, &info);`	831	`work, &lwork, &info);`
832	`if (info < 0)`	832	`if (info < 0)`
833	`error(_("error code %d from Lapack routine %s"), info, "dgeqp3");`	833	`error(_("error code %d from Lapack routine '%s'"), info, "dgeqp3");`
834	`val = PROTECT(allocVector(VECSXP, 4));`	834	`val = PROTECT(allocVector(VECSXP, 4));`
835	`nm = PROTECT(allocVector(STRSXP, 4));`	835	`nm = PROTECT(allocVector(STRSXP, 4));`
836	`rank = PROTECT(allocVector(INTSXP, 1));`	836	`rank = PROTECT(allocVector(INTSXP, 1));`
837	`INTEGER(rank)[0] = m < n ? m : n;`	837	`INTEGER(rank)[0] = m < n ? m : n;`
838	`SET_STRING_ELT(nm, 0, mkChar("qr"));`	838	`SET_STRING_ELT(nm, 0, mkChar("qr"));`
Line 868...	Line 868...
868	`lwork = -1;`	868	`lwork = -1;`
869	`F77_CALL(dormqr)("L", "T", &n, &nrhs, &k,`	869	`F77_CALL(dormqr)("L", "T", &n, &nrhs, &k,`
870	`REAL(qr), &n, REAL(tau), REAL(B), &n,`	870	`REAL(qr), &n, REAL(tau), REAL(B), &n,`
871	`&tmp, &lwork, &info);`	871	`&tmp, &lwork, &info);`
872	`if (info != 0)`	872	`if (info != 0)`
873	`error(_("error code %d from Lapack routine %s"), info, "dormqr");`	873	`error(_("error code %d from Lapack routine '%s'"), info, "dormqr");`
874	`lwork = (int) tmp;`	874	`lwork = (int) tmp;`
875	`work = (double *) R_alloc(lwork, sizeof(double));`	875	`work = (double *) R_alloc(lwork, sizeof(double));`
876	`F77_CALL(dormqr)("L", "T", &n, &nrhs, &k,`	876	`F77_CALL(dormqr)("L", "T", &n, &nrhs, &k,`
877	`REAL(qr), &n, REAL(tau), REAL(B), &n,`	877	`REAL(qr), &n, REAL(tau), REAL(B), &n,`
878	`work, &lwork, &info);`	878	`work, &lwork, &info);`
879	`if (info != 0)`	879	`if (info != 0)`
880	`error(_("error code %d from Lapack routine %s"), info, "dormqr");`	880	`error(_("error code %d from Lapack routine '%s'"), info, "dormqr");`
881	`F77_CALL(dtrtrs)("U", "N", "N", &k, &nrhs,`	881	`F77_CALL(dtrtrs)("U", "N", "N", &k, &nrhs,`
882	`REAL(qr), &n, REAL(B), &n, &info);`	882	`REAL(qr), &n, REAL(B), &n, &info);`
883	`if (info != 0)`	883	`if (info != 0)`
884	`error(_("error code %d from Lapack routine %s"), info, "dtrtrs");`	884	`error(_("error code %d from Lapack routine '%s'"), info, "dtrtrs");`
885	`UNPROTECT(1);`	885	`UNPROTECT(1);`
886	`return B;`	886	`return B;`
887	`}`	887	`}`
888		888
889	`static SEXP modqr_qy_real(SEXP Q, SEXP Bin, SEXP trans)`	889	`static SEXP modqr_qy_real(SEXP Q, SEXP Bin, SEXP trans)`
Line 908...	Line 908...
908	`lwork = -1;`	908	`lwork = -1;`
909	`F77_CALL(dormqr)("L", tr ? "T" : "N", &n, &nrhs, &k,`	909	`F77_CALL(dormqr)("L", tr ? "T" : "N", &n, &nrhs, &k,`
910	`REAL(qr), &n, REAL(tau), REAL(B), &n,`	910	`REAL(qr), &n, REAL(tau), REAL(B), &n,`
911	`&tmp, &lwork, &info);`	911	`&tmp, &lwork, &info);`
912	`if (info != 0)`	912	`if (info != 0)`
913	`error(_("error code %d from Lapack routine %s"), info, "dormqr");`	913	`error(_("error code %d from Lapack routine '%s'"), info, "dormqr");`
914	`lwork = (int) tmp;`	914	`lwork = (int) tmp;`
915	`work = (double *) R_alloc(lwork, sizeof(double));`	915	`work = (double *) R_alloc(lwork, sizeof(double));`
916	`F77_CALL(dormqr)("L", tr ? "T" : "N", &n, &nrhs, &k,`	916	`F77_CALL(dormqr)("L", tr ? "T" : "N", &n, &nrhs, &k,`
917	`REAL(qr), &n, REAL(tau), REAL(B), &n,`	917	`REAL(qr), &n, REAL(tau), REAL(B), &n,`
918	`work, &lwork, &info);`	918	`work, &lwork, &info);`
919	`if (info != 0)`	919	`if (info != 0)`
920	`error(_("error code %d from Lapack routine %s"), info, "dormqr");`	920	`error(_("error code %d from Lapack routine '%s'"), info, "dormqr");`
921	`UNPROTECT(1);`	921	`UNPROTECT(1);`
922	`return B;`	922	`return B;`
923	`}`	923	`}`
924		924
925	`static SEXP moddet_ge_real(SEXP Ain, SEXP logarithm)`	925	`static SEXP moddet_ge_real(SEXP Ain, SEXP logarithm)`
Line 939...	Line 939...
939	`error(_("'A' must be a square matrix"));`	939	`error(_("'A' must be a square matrix"));`
940	`jpvt = (int *) R_alloc(n, sizeof(int));`	940	`jpvt = (int *) R_alloc(n, sizeof(int));`
941	`F77_CALL(dgetrf)(&n, &n, REAL(A), &n, jpvt, &info);`	941	`F77_CALL(dgetrf)(&n, &n, REAL(A), &n, jpvt, &info);`
942	`sign = 1;`	942	`sign = 1;`
943	`if (info < 0)`	943	`if (info < 0)`
944	`error(_("error code %d from Lapack routine %s"), info, "dgetrf");`	944	`error(_("error code %d from Lapack routine '%s'"), info, "dgetrf");`
945	`else if (info > 0) { /* Singular matrix: U[i,i] (i := info) is 0 */`	945	`else if (info > 0) { /* Singular matrix: U[i,i] (i := info) is 0 */`
946	`/warning("Lapack dgetrf(): singular matrix: U[%d,%d]=0", info,info);/`	946	`/warning("Lapack dgetrf(): singular matrix: U[%d,%d]=0", info,info);/`
947	`modulus = (useLog ? R_NegInf : 0.);`	947	`modulus = (useLog ? R_NegInf : 0.);`
948	`}`	948	`}`
949	`else {`	949	`else {`

The R Project SVN R

(root)/trunk/src/modules/lapack/Lapack.c – Rev 33282 → 33387