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

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