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 342 - (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 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("poMatrix")));
80 :     int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
81 : bates 342 *vDims;
82 :     int n = Dims[1];
83 : bates 10 double one = 1.0, zero = 0.0;
84 :    
85 : bates 332 SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
86 :     SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
87 :     SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar("U")));
88 : bates 342 SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
89 :     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
90 : bates 332 vDims[0] = vDims[1] = n;
91 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, n * n));
92 : bates 342 if (Dims[0] > 0) {
93 :     F77_CALL(dsyrk)("U", "T", vDims, Dims,
94 :     &one, REAL(GET_SLOT(x, Matrix_xSym)), Dims,
95 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)), vDims);
96 : bates 332 } else {
97 :     int i, nsqr = n * n;
98 :     double *xvals = REAL(GET_SLOT(val, Matrix_xSym));
99 :     for (i = 0; i < nsqr; i++) xvals[i] = 0.;
100 : bates 10 }
101 : bates 332
102 : bates 10 UNPROTECT(1);
103 :     return val;
104 :     }
105 :    
106 :     SEXP geMatrix_geMatrix_crossprod(SEXP x, SEXP y)
107 :     {
108 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));
109 :     int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
110 :     *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),
111 : bates 332 *vDims;
112 : bates 10 int m = xDims[1], n = yDims[1];
113 :     double one = 1.0, zero = 0.0;
114 :    
115 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
116 :     SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
117 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
118 :     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
119 : bates 10 if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
120 :     if (*xDims != *yDims)
121 :     error("Dimensions of x and y are not compatible for crossprod");
122 :     vDims[0] = m; vDims[1] = n;
123 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
124 :     F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
125 :     REAL(GET_SLOT(x, Matrix_xSym)), xDims,
126 :     REAL(GET_SLOT(y, Matrix_xSym)), yDims,
127 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)),
128 :     xDims + 1);
129 :     }
130 :     UNPROTECT(1);
131 :     return val;
132 :     }
133 :    
134 :     SEXP geMatrix_matrix_crossprod(SEXP x, SEXP y)
135 :     {
136 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));
137 :     int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
138 :     *yDims = INTEGER(getAttrib(y, R_DimSymbol)),
139 : bates 332 *vDims;
140 : bates 10 int m = xDims[1], n = yDims[1];
141 :     double one = 1.0, zero = 0.0;
142 :    
143 :     if (!(isMatrix(y) && isReal(y)))
144 :     error("Argument y must be a numeric (real) matrix");
145 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
146 :     SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
147 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
148 :     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
149 : bates 10 if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
150 :     if (*xDims != *yDims)
151 :     error("Dimensions of x and y are not compatible for crossprod");
152 :     vDims[0] = m; vDims[1] = n;
153 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
154 :     F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
155 :     REAL(GET_SLOT(x, Matrix_xSym)), xDims,
156 :     REAL(y), yDims,
157 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)),
158 :     xDims + 1);
159 :     }
160 :     UNPROTECT(1);
161 :     return val;
162 :     }
163 :    
164 :     SEXP geMatrix_getDiag(SEXP x)
165 :     {
166 :     int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
167 :     int i, m = dims[0], nret = (m < dims[1]) ? m : dims[1];
168 :     SEXP ret = PROTECT(allocVector(REALSXP, nret)),
169 :     xv = GET_SLOT(x, Matrix_xSym);
170 :    
171 :     for (i = 0; i < nret; i++) {
172 :     REAL(ret)[i] = REAL(xv)[i * (m + 1)];
173 :     }
174 :     UNPROTECT(1);
175 :     return ret;
176 :     }
177 :    
178 :     SEXP geMatrix_LU(SEXP x)
179 :     {
180 :     SEXP val = get_factorization(x, "LU");
181 :     int *dims, npiv, info;
182 :    
183 :     if (val != R_NilValue) return val;
184 :     dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
185 :     if (dims[0] < 1 || dims[1] < 1)
186 :     error("Cannot factor a matrix with zero extents");
187 :     npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];
188 :     val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));
189 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
190 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
191 :     SET_SLOT(val, install("pivot"), allocVector(INTSXP, npiv));
192 :     F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),
193 :     dims, INTEGER(GET_SLOT(val, install("pivot"))),
194 :     &info);
195 :     if (info) error("Lapack routine dgetrf returned error code %d", info);
196 :     UNPROTECT(1);
197 :     return set_factorization(x, val, "LU");
198 :     }
199 :    
200 :     SEXP geMatrix_determinant(SEXP x, SEXP logarithm)
201 :     {
202 :     int lg = asLogical(logarithm);
203 :     SEXP lu = geMatrix_LU(x);
204 :     int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
205 :     *jpvt = INTEGER(GET_SLOT(lu, install("pivot"))),
206 :     i, n = dims[0], sign = 1;
207 :     double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;
208 :    
209 :     if (n != dims[1])
210 :     error("Determinant requires a square matrix");
211 :     for (i = 0; i < n; i++) if (jpvt[i] != (i + 1)) sign = -sign;
212 :     if (lg) {
213 :     modulus = 0.0;
214 :     for (i = 0; i < n; i++) {
215 :     double dii = luvals[i*(n + 1)]; /* ith diagonal element */
216 :     modulus += log(dii < 0 ? -dii : dii);
217 :     if (dii < 0) sign = -sign;
218 :     }
219 :     } else {
220 :     modulus = 1.0;
221 :     for (i = 0; i < n; i++)
222 :     modulus *= luvals[i*(n + 1)];
223 :     if (modulus < 0) {
224 :     modulus = -modulus;
225 :     sign = -sign;
226 :     }
227 :     }
228 :     return as_det_obj(modulus, lg, sign);
229 :     }
230 :    
231 :     SEXP geMatrix_solve(SEXP a)
232 :     {
233 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix"))),
234 :     lu = geMatrix_LU(a);
235 :     int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
236 :     *pivot = INTEGER(GET_SLOT(lu, install("pivot")));
237 :     double *x, tmp;
238 :     int info, lwork = -1;
239 :    
240 :    
241 :     if (dims[0] != dims[1]) error("Solve requires a square matrix");
242 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
243 :     SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
244 : bates 10 SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));
245 :     x = REAL(GET_SLOT(val, Matrix_xSym));
246 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(lu, Matrix_DimSym)));
247 :     F77_CALL(dgetri)(dims, x, dims, pivot, &tmp, &lwork, &info);
248 :     lwork = (int) tmp;
249 :     F77_CALL(dgetri)(dims, x, dims, pivot,
250 :     (double *) R_alloc((size_t) lwork, sizeof(double)),
251 :     &lwork, &info);
252 :     UNPROTECT(1);
253 :     return val;
254 :     }
255 :    
256 :     SEXP geMatrix_geMatrix_mm(SEXP a, SEXP b)
257 :     {
258 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
259 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
260 :     *cdims,
261 :     m = adims[0], n = bdims[1], k = adims[1];
262 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));
263 :     char *trans = "N";
264 :     double one = 1., zero = 0.;
265 :    
266 :     if (bdims[0] != k)
267 :     error("Matrices are not conformable for multiplication");
268 :     if (m < 1 || n < 1 || k < 1)
269 :     error("Matrices with zero extents cannot be multiplied");
270 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
271 :     SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
272 : bates 10 SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
273 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
274 :     cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
275 :     cdims[0] = m; cdims[1] = n;
276 :     F77_CALL(dgemm)(trans, trans, adims, bdims+1, bdims, &one,
277 :     REAL(GET_SLOT(a, Matrix_xSym)), adims,
278 :     REAL(GET_SLOT(b, Matrix_xSym)), bdims,
279 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);
280 :     UNPROTECT(1);
281 :     return val;
282 :     }
283 : bates 301
284 :     SEXP geMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
285 :     {
286 :     int nu = asInteger(nnu),
287 :     nv = asInteger(nnv),
288 :     *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
289 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym));
290 :     SEXP val = PROTECT(allocVector(VECSXP, 3));
291 :    
292 :     if (dims[0] && dims[1]) {
293 :     int m = dims[0], n = dims[1], mm = (m < n)?m:n,
294 :     lwork = -1, info;
295 :     int *iwork = Calloc(8 * mm, int);
296 :     double tmp, *work;
297 : bates 306 /* int bdspac = 3*m*m + 4*m, */
298 :     /* wrkbl, maxwrk, minwrk, itmp, */
299 :     /* ione = 1, iminus1 = -1; */
300 :     /* int i1, i2, i3; */
301 : bates 301
302 :     SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));
303 :     SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));
304 :     SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));
305 :     F77_CALL(dgesdd)("S", &m, &n, xx, &m,
306 :     REAL(VECTOR_ELT(val, 0)),
307 :     REAL(VECTOR_ELT(val, 1)), &m,
308 :     REAL(VECTOR_ELT(val, 2)), &mm,
309 :     &tmp, &lwork, iwork, &info);
310 :     lwork = (int) tmp;
311 : bates 306 /* F77_CALL(foo)(&i1, &i2, &i3); */
312 :     /* wrkbl = 3*m+(m+n)*i1; */
313 :     /* if (wrkbl < (itmp = 3*m + m*i2)) wrkbl = itmp; */
314 :     /* if (wrkbl < (itmp = 3*m + m*i3)) wrkbl = itmp; */
315 :     /* itmp = bdspac+3*m; */
316 :     /* maxwrk = (wrkbl > itmp) ? wrkbl : itmp; */
317 :     /* minwrk = 3*m + ((bdspac > n) ? bdspac : n); */
318 : bates 301 work = Calloc(lwork, double);
319 :     F77_CALL(dgesdd)("S", &m, &n, xx, &m,
320 :     REAL(VECTOR_ELT(val, 0)),
321 :     REAL(VECTOR_ELT(val, 1)), &m,
322 :     REAL(VECTOR_ELT(val, 2)), &mm,
323 :     work, &lwork, iwork, &info);
324 :     Free(iwork); Free(work);
325 :     }
326 :     UNPROTECT(1);
327 :     return val;
328 :     }

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