SCM

SCM Repository

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

Annotation of /pkg/src/dgeMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 343 - (view) (download) (as text)
Original Path: pkg/src/geMatrix.c

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

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