SCM

SCM Repository

[matrix] Annotation of /pkg/src/geMatrix.c
ViewVC logotype

Annotation of /pkg/src/geMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 345 - (view) (download) (as text)

1 : bates 10 #include "geMatrix.h"
2 :    
3 :     SEXP geMatrix_validate(SEXP obj)
4 :     {
5 :     SEXP x = GET_SLOT(obj, Matrix_xSym),
6 : bates 332 Dim = GET_SLOT(obj, Matrix_DimSym),
7 :     fact = GET_SLOT(obj, Matrix_factorization),
8 :     rc = GET_SLOT(obj, Matrix_rcondSym);
9 : bates 10 int m, n;
10 :    
11 :     if (length(Dim) != 2)
12 :     return ScalarString(mkChar("Dim slot must have length 2"));
13 :     m = INTEGER(Dim)[0]; n = INTEGER(Dim)[1];
14 :     if (m < 0 || n < 0)
15 : bates 332 return ScalarString(mkChar("Negative value(s) in Dim"));
16 : bates 10 if (length(x) != m * n)
17 : bates 332 return ScalarString(mkChar("length of x slot != prod(Dim)"));
18 :     if (length(fact) > 0 && getAttrib(fact, R_NamesSymbol) == R_NilValue)
19 :     return ScalarString(mkChar("factorization slot must be named list"));
20 :     if (length(rc) > 0 && getAttrib(rc, R_NamesSymbol) == R_NilValue)
21 :     return ScalarString(mkChar("rcond slot must be named numeric vector"));
22 : bates 10 return ScalarLogical(1);
23 :     }
24 :    
25 :     static
26 :     double get_norm(SEXP obj, char *typstr)
27 :     {
28 :     char typnm[] = {'\0', '\0'};
29 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
30 :     double *work = (double *) NULL;
31 :    
32 :     typnm[0] = norm_type(typstr);
33 :     if (*typnm == 'I') {
34 :     work = (double *) R_alloc(dims[0], sizeof(double));
35 :     }
36 :     return F77_CALL(dlange)(typstr, dims, dims+1,
37 :     REAL(GET_SLOT(obj, Matrix_xSym)),
38 :     dims, work);
39 :     }
40 :    
41 :     SEXP geMatrix_norm(SEXP obj, SEXP type)
42 :     {
43 :     return ScalarReal(get_norm(obj, CHAR(asChar(type))));
44 :     }
45 :    
46 :     static
47 :     double set_rcond(SEXP obj, char *typstr)
48 :     {
49 :     char typnm[] = {'\0', '\0'};
50 : bates 332 SEXP rcv = GET_SLOT(obj, Matrix_rcondSym);
51 : bates 10 double rcond = get_double_by_name(rcv, typnm);
52 :    
53 :     typnm[0] = rcond_type(typstr);
54 :     if (R_IsNA(rcond)) {
55 :     SEXP LU = geMatrix_LU(obj);
56 :     int *dims = INTEGER(GET_SLOT(LU, Matrix_DimSym)), info;
57 :     double anorm = get_norm(obj, typstr);
58 :    
59 :     if (dims[0] != dims[1] || dims[0] < 1)
60 :     error("rcond requires a square, non-empty matrix");
61 :     F77_CALL(dgecon)(typnm,
62 :     dims, REAL(GET_SLOT(LU, Matrix_xSym)),
63 :     dims, &anorm, &rcond,
64 :     (double *) R_alloc(4*dims[0], sizeof(double)),
65 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
66 : bates 332 SET_SLOT(obj, Matrix_rcondSym,
67 : bates 10 set_double_by_name(rcv, rcond, typnm));
68 :     }
69 :     return rcond;
70 :     }
71 :    
72 :     SEXP geMatrix_rcond(SEXP obj, SEXP type)
73 :     {
74 :     return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
75 :     }
76 :    
77 :     SEXP geMatrix_crossprod(SEXP x)
78 :     {
79 : bates 345 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("poMatrix"))),
80 :     xslot;
81 : bates 10 int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
82 : bates 342 *vDims;
83 : bates 343 int i, n = Dims[1];
84 :     int nsqr = n * n;
85 :     double one = 1.0, *xvals, zero = 0.0;
86 : bates 10
87 : bates 332 SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
88 :     SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
89 :     SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar("U")));
90 : bates 342 SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
91 :     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
92 : bates 332 vDims[0] = vDims[1] = n;
93 : bates 345 SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nsqr));
94 : bates 343 xvals = REAL(GET_SLOT(val, Matrix_xSym));
95 :     for (i = 0; i < nsqr; i++) xvals[i] = 0.; /* keep valgrind happy */
96 : bates 345 F77_CALL(dsyrk)("U", "T", vDims, Dims,
97 :     &one, REAL(GET_SLOT(x, Matrix_xSym)), Dims,
98 :     &zero, xvals, vDims);
99 : bates 10 UNPROTECT(1);
100 :     return val;
101 :     }
102 :    
103 :     SEXP geMatrix_geMatrix_crossprod(SEXP x, SEXP y)
104 :     {
105 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));
106 :     int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
107 :     *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),
108 : bates 332 *vDims;
109 : bates 10 int m = xDims[1], n = yDims[1];
110 :     double one = 1.0, zero = 0.0;
111 :    
112 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
113 :     SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
114 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
115 :     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
116 : bates 10 if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
117 :     if (*xDims != *yDims)
118 :     error("Dimensions of x and y are not compatible for crossprod");
119 :     vDims[0] = m; vDims[1] = n;
120 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
121 :     F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
122 :     REAL(GET_SLOT(x, Matrix_xSym)), xDims,
123 :     REAL(GET_SLOT(y, Matrix_xSym)), yDims,
124 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)),
125 :     xDims + 1);
126 :     }
127 :     UNPROTECT(1);
128 :     return val;
129 :     }
130 :    
131 :     SEXP geMatrix_matrix_crossprod(SEXP x, SEXP y)
132 :     {
133 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));
134 :     int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
135 :     *yDims = INTEGER(getAttrib(y, R_DimSymbol)),
136 : bates 332 *vDims;
137 : bates 10 int m = xDims[1], n = yDims[1];
138 :     double one = 1.0, zero = 0.0;
139 :    
140 :     if (!(isMatrix(y) && isReal(y)))
141 :     error("Argument y must be a numeric (real) matrix");
142 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
143 :     SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
144 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
145 :     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
146 : bates 10 if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
147 :     if (*xDims != *yDims)
148 :     error("Dimensions of x and y are not compatible for crossprod");
149 :     vDims[0] = m; vDims[1] = n;
150 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
151 :     F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
152 :     REAL(GET_SLOT(x, Matrix_xSym)), xDims,
153 :     REAL(y), yDims,
154 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)),
155 :     xDims + 1);
156 :     }
157 :     UNPROTECT(1);
158 :     return val;
159 :     }
160 :    
161 :     SEXP geMatrix_getDiag(SEXP x)
162 :     {
163 :     int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
164 :     int i, m = dims[0], nret = (m < dims[1]) ? m : dims[1];
165 :     SEXP ret = PROTECT(allocVector(REALSXP, nret)),
166 :     xv = GET_SLOT(x, Matrix_xSym);
167 :    
168 :     for (i = 0; i < nret; i++) {
169 :     REAL(ret)[i] = REAL(xv)[i * (m + 1)];
170 :     }
171 :     UNPROTECT(1);
172 :     return ret;
173 :     }
174 :    
175 :     SEXP geMatrix_LU(SEXP x)
176 :     {
177 :     SEXP val = get_factorization(x, "LU");
178 :     int *dims, npiv, info;
179 :    
180 :     if (val != R_NilValue) return val;
181 :     dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
182 :     if (dims[0] < 1 || dims[1] < 1)
183 :     error("Cannot factor a matrix with zero extents");
184 :     npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];
185 :     val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));
186 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
187 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
188 :     SET_SLOT(val, install("pivot"), allocVector(INTSXP, npiv));
189 :     F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),
190 :     dims, INTEGER(GET_SLOT(val, install("pivot"))),
191 :     &info);
192 :     if (info) error("Lapack routine dgetrf returned error code %d", info);
193 :     UNPROTECT(1);
194 :     return set_factorization(x, val, "LU");
195 :     }
196 :    
197 :     SEXP geMatrix_determinant(SEXP x, SEXP logarithm)
198 :     {
199 :     int lg = asLogical(logarithm);
200 :     SEXP lu = geMatrix_LU(x);
201 :     int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
202 :     *jpvt = INTEGER(GET_SLOT(lu, install("pivot"))),
203 :     i, n = dims[0], sign = 1;
204 :     double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;
205 :    
206 :     if (n != dims[1])
207 :     error("Determinant requires a square matrix");
208 :     for (i = 0; i < n; i++) if (jpvt[i] != (i + 1)) sign = -sign;
209 :     if (lg) {
210 :     modulus = 0.0;
211 :     for (i = 0; i < n; i++) {
212 :     double dii = luvals[i*(n + 1)]; /* ith diagonal element */
213 :     modulus += log(dii < 0 ? -dii : dii);
214 :     if (dii < 0) sign = -sign;
215 :     }
216 :     } else {
217 :     modulus = 1.0;
218 :     for (i = 0; i < n; i++)
219 :     modulus *= luvals[i*(n + 1)];
220 :     if (modulus < 0) {
221 :     modulus = -modulus;
222 :     sign = -sign;
223 :     }
224 :     }
225 :     return as_det_obj(modulus, lg, sign);
226 :     }
227 :    
228 :     SEXP geMatrix_solve(SEXP a)
229 :     {
230 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix"))),
231 :     lu = geMatrix_LU(a);
232 :     int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
233 :     *pivot = INTEGER(GET_SLOT(lu, install("pivot")));
234 :     double *x, tmp;
235 :     int info, lwork = -1;
236 :    
237 :    
238 :     if (dims[0] != dims[1]) error("Solve requires a square matrix");
239 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
240 :     SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
241 : bates 10 SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));
242 :     x = REAL(GET_SLOT(val, Matrix_xSym));
243 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(lu, Matrix_DimSym)));
244 :     F77_CALL(dgetri)(dims, x, dims, pivot, &tmp, &lwork, &info);
245 :     lwork = (int) tmp;
246 :     F77_CALL(dgetri)(dims, x, dims, pivot,
247 :     (double *) R_alloc((size_t) lwork, sizeof(double)),
248 :     &lwork, &info);
249 :     UNPROTECT(1);
250 :     return val;
251 :     }
252 :    
253 :     SEXP geMatrix_geMatrix_mm(SEXP a, SEXP b)
254 :     {
255 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
256 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
257 :     *cdims,
258 :     m = adims[0], n = bdims[1], k = adims[1];
259 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));
260 :     char *trans = "N";
261 :     double one = 1., zero = 0.;
262 :    
263 :     if (bdims[0] != k)
264 :     error("Matrices are not conformable for multiplication");
265 :     if (m < 1 || n < 1 || k < 1)
266 :     error("Matrices with zero extents cannot be multiplied");
267 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
268 :     SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
269 : bates 10 SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
270 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
271 :     cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
272 :     cdims[0] = m; cdims[1] = n;
273 :     F77_CALL(dgemm)(trans, trans, adims, bdims+1, bdims, &one,
274 :     REAL(GET_SLOT(a, Matrix_xSym)), adims,
275 :     REAL(GET_SLOT(b, Matrix_xSym)), bdims,
276 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);
277 :     UNPROTECT(1);
278 :     return val;
279 :     }
280 : bates 301
281 :     SEXP geMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
282 :     {
283 :     int nu = asInteger(nnu),
284 :     nv = asInteger(nnv),
285 :     *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
286 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym));
287 :     SEXP val = PROTECT(allocVector(VECSXP, 3));
288 :    
289 :     if (dims[0] && dims[1]) {
290 :     int m = dims[0], n = dims[1], mm = (m < n)?m:n,
291 :     lwork = -1, info;
292 :     int *iwork = Calloc(8 * mm, int);
293 :     double tmp, *work;
294 : bates 306 /* int bdspac = 3*m*m + 4*m, */
295 :     /* wrkbl, maxwrk, minwrk, itmp, */
296 :     /* ione = 1, iminus1 = -1; */
297 :     /* int i1, i2, i3; */
298 : bates 301
299 :     SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));
300 :     SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));
301 :     SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));
302 :     F77_CALL(dgesdd)("S", &m, &n, xx, &m,
303 :     REAL(VECTOR_ELT(val, 0)),
304 :     REAL(VECTOR_ELT(val, 1)), &m,
305 :     REAL(VECTOR_ELT(val, 2)), &mm,
306 :     &tmp, &lwork, iwork, &info);
307 :     lwork = (int) tmp;
308 : bates 306 /* F77_CALL(foo)(&i1, &i2, &i3); */
309 :     /* wrkbl = 3*m+(m+n)*i1; */
310 :     /* if (wrkbl < (itmp = 3*m + m*i2)) wrkbl = itmp; */
311 :     /* if (wrkbl < (itmp = 3*m + m*i3)) wrkbl = itmp; */
312 :     /* itmp = bdspac+3*m; */
313 :     /* maxwrk = (wrkbl > itmp) ? wrkbl : itmp; */
314 :     /* minwrk = 3*m + ((bdspac > n) ? bdspac : n); */
315 : bates 301 work = Calloc(lwork, double);
316 :     F77_CALL(dgesdd)("S", &m, &n, xx, &m,
317 :     REAL(VECTOR_ELT(val, 0)),
318 :     REAL(VECTOR_ELT(val, 1)), &m,
319 :     REAL(VECTOR_ELT(val, 2)), &mm,
320 :     work, &lwork, iwork, &info);
321 :     Free(iwork); Free(work);
322 :     }
323 :     UNPROTECT(1);
324 :     return val;
325 :     }

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge