SCM

SCM Repository

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

Annotation of /pkg/src/dsyMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 478 #include "dsyMatrix.h"
2 : bates 10
3 : maechler 890 SEXP symmetricMatrix_validate(SEXP obj)
4 : bates 10 {
5 : maechler 890 SEXP val = GET_SLOT(obj, Matrix_DimSym);
6 :     if (LENGTH(val) < 2)
7 :     return mkString(_("'Dim' slot has length less than two"));
8 :     if (INTEGER(val)[0] != INTEGER(val)[1])
9 :     return mkString(_("Matrix is not square"));
10 : bates 592 if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
11 :     "LU", "uplo"))) return val;
12 : bates 10 return ScalarLogical(1);
13 :     }
14 :    
15 : maechler 890 SEXP dsyMatrix_validate(SEXP obj)
16 :     {
17 :     return symmetricMatrix_validate(obj);
18 :     /* see ./dspMatrix.c for how to do further checks */
19 :     }
20 :    
21 : bates 296 double get_norm_sy(SEXP obj, char *typstr)
22 : maechler 534 {
23 : bates 296 char typnm[] = {'\0', '\0'};
24 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
25 :     double *work = (double *) NULL;
26 :    
27 :     typnm[0] = norm_type(typstr);
28 :     if (*typnm == 'I' || *typnm == 'O') {
29 :     work = (double *) R_alloc(dims[0], sizeof(double));
30 :     }
31 : maechler 951 return F77_CALL(dlansy)(typnm, uplo_P(obj),
32 : bates 296 dims, REAL(GET_SLOT(obj, Matrix_xSym)),
33 :     dims, work);
34 :     }
35 :    
36 : bates 478 SEXP dsyMatrix_norm(SEXP obj, SEXP type)
37 : bates 296 {
38 :     return ScalarReal(get_norm_sy(obj, CHAR(asChar(type))));
39 :     }
40 :    
41 :     static
42 :     double set_rcond_sy(SEXP obj, char *typstr)
43 :     {
44 :     char typnm[] = {'\0', '\0'};
45 : bates 338 SEXP rcv = GET_SLOT(obj, Matrix_rcondSym);
46 : bates 296 double rcond;
47 :    
48 :     typnm[0] = rcond_type(typstr);
49 :     rcond = get_double_by_name(rcv, typnm);
50 :    
51 :     if (R_IsNA(rcond)) {
52 : bates 643 SEXP trf = dsyMatrix_trf(obj);
53 : bates 296 int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
54 :     double anorm = get_norm_sy(obj, "O");
55 :    
56 : maechler 951 F77_CALL(dsycon)(uplo_P(trf), dims,
57 :     REAL (GET_SLOT(trf, Matrix_xSym)), dims,
58 :     INTEGER(GET_SLOT(trf, Matrix_permSym)),
59 : maechler 534 &anorm, &rcond,
60 : bates 296 (double *) R_alloc(2*dims[0], sizeof(double)),
61 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
62 : bates 338 SET_SLOT(obj, Matrix_rcondSym,
63 : bates 296 set_double_by_name(rcv, rcond, typnm));
64 :     }
65 :     return rcond;
66 :     }
67 :    
68 : bates 478 SEXP dsyMatrix_rcond(SEXP obj, SEXP type)
69 : bates 296 {
70 : bates 643 return ScalarReal(set_rcond_sy(obj, CHAR(asChar(type))));
71 : bates 296 }
72 :    
73 : maechler 534 static
74 : bates 10 void make_symmetric(double *to, SEXP from, int n)
75 :     {
76 :     int i, j;
77 : maechler 951 if (*uplo_P(from) == 'U') {
78 : bates 10 for (j = 0; j < n; j++) {
79 :     for (i = j+1; i < n; i++) {
80 :     to[i + j*n] = to[j + i*n];
81 :     }
82 :     }
83 :     } else {
84 :     for (j = 1; j < n; j++) {
85 :     for (i = 0; i < j && i < n; i++) {
86 :     to[i + j*n] = to[j + i*n];
87 :     }
88 :     }
89 :     }
90 :     }
91 : maechler 534
92 : bates 478 SEXP dsyMatrix_solve(SEXP a)
93 : bates 296 {
94 : bates 643 SEXP trf = dsyMatrix_trf(a);
95 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix")));
96 :     int *dims = INTEGER(GET_SLOT(trf, Matrix_DimSym)), info;
97 :    
98 :     SET_SLOT(val, Matrix_uploSym, duplicate(GET_SLOT(trf, Matrix_uploSym)));
99 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(trf, Matrix_xSym)));
100 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(trf, Matrix_DimSym)));
101 :     SET_SLOT(val, Matrix_rcondSym, duplicate(GET_SLOT(a, Matrix_rcondSym)));
102 : maechler 951 F77_CALL(dsytri)(uplo_P(val), dims,
103 :     REAL(GET_SLOT(val, Matrix_xSym)), dims,
104 : bates 643 INTEGER(GET_SLOT(trf, Matrix_permSym)),
105 :     (double *) R_alloc((long) dims[0], sizeof(double)),
106 :     &info);
107 :     UNPROTECT(1);
108 :     return val;
109 : bates 296 }
110 :    
111 : bates 643 SEXP dsyMatrix_dgeMatrix_solve(SEXP a, SEXP b)
112 :     {
113 :     SEXP trf = dsyMatrix_trf(a),
114 :     val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
115 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
116 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
117 :     info;
118 :    
119 :     if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
120 :     error(_("Dimensions of system to be solved are inconsistent"));
121 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(b, Matrix_DimSym)));
122 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(b, Matrix_xSym)));
123 : maechler 951 F77_CALL(dsytrs)(uplo_P(trf), adims, bdims + 1,
124 : bates 643 REAL(GET_SLOT(trf, Matrix_xSym)), adims,
125 :     INTEGER(GET_SLOT(trf, Matrix_permSym)),
126 :     REAL(GET_SLOT(val, Matrix_xSym)),
127 :     bdims, &info);
128 :     UNPROTECT(1);
129 :     return val;
130 :     }
131 :    
132 : bates 478 SEXP dsyMatrix_matrix_solve(SEXP a, SEXP b)
133 : bates 296 {
134 : bates 643 SEXP trf = dsyMatrix_trf(a),
135 :     val = PROTECT(duplicate(b));
136 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
137 :     *bdims = INTEGER(getAttrib(b, R_DimSymbol)),
138 :     info;
139 :    
140 :     if (!(isReal(b) && isMatrix(b)))
141 :     error(_("Argument b must be a numeric matrix"));
142 :     if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
143 :     error(_("Dimensions of system to be solved are inconsistent"));
144 : maechler 951 F77_CALL(dsytrs)(uplo_P(trf), adims, bdims + 1,
145 : bates 643 REAL(GET_SLOT(trf, Matrix_xSym)), adims,
146 :     INTEGER(GET_SLOT(trf, Matrix_permSym)),
147 :     REAL(val), bdims, &info);
148 :     UNPROTECT(1);
149 :     return val;
150 : bates 296 }
151 :    
152 : bates 478 SEXP dsyMatrix_as_dgeMatrix(SEXP from)
153 : bates 10 {
154 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
155 : bates 338 rcondSym = Matrix_rcondSym;
156 : maechler 534
157 : bates 338 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
158 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
159 : bates 10 SET_SLOT(val, rcondSym, duplicate(GET_SLOT(from, rcondSym)));
160 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(from, Matrix_xSym)));
161 :     SET_SLOT(val, Matrix_DimSym,
162 :     duplicate(GET_SLOT(from, Matrix_DimSym)));
163 :     make_symmetric(REAL(GET_SLOT(val, Matrix_xSym)), from,
164 :     INTEGER(GET_SLOT(val, Matrix_DimSym))[0]);
165 :     UNPROTECT(1);
166 :     return val;
167 :     }
168 :    
169 : bates 478 SEXP dsyMatrix_as_matrix(SEXP from)
170 : bates 10 {
171 :     int n = INTEGER(GET_SLOT(from, Matrix_DimSym))[0];
172 :     SEXP val = PROTECT(allocMatrix(REALSXP, n, n));
173 : maechler 534
174 : bates 10 make_symmetric(Memcpy(REAL(val),
175 :     REAL(GET_SLOT(from, Matrix_xSym)), n * n),
176 :     from, n);
177 : maechler 628 setAttrib(val, R_DimNamesSymbol, GET_SLOT(from, Matrix_DimNamesSym));
178 : bates 10 UNPROTECT(1);
179 :     return val;
180 :     }
181 :    
182 : bates 478 SEXP dsyMatrix_dgeMatrix_mm(SEXP a, SEXP b)
183 : bates 10 {
184 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
185 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
186 :     *cdims,
187 :     m = adims[0], n = bdims[1], k = adims[1];
188 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
189 : bates 10 double one = 1., zero = 0.;
190 : maechler 534
191 : bates 10 if (bdims[0] != k)
192 : bates 582 error(_("Matrices are not conformable for multiplication"));
193 : bates 10 if (m < 1 || n < 1 || k < 1)
194 : bates 582 error(_("Matrices with zero extents cannot be multiplied"));
195 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
196 : bates 338 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
197 : bates 10 SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
198 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
199 :     cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
200 :     cdims[0] = m; cdims[1] = n;
201 : maechler 951 F77_CALL(dsymm)("L", uplo_P(a), adims, bdims+1, &one,
202 : bates 10 REAL(GET_SLOT(a, Matrix_xSym)), adims,
203 :     REAL(GET_SLOT(b, Matrix_xSym)), bdims,
204 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);
205 :     UNPROTECT(1);
206 :     return val;
207 :     }
208 :    
209 : bates 478 SEXP dsyMatrix_dgeMatrix_mm_R(SEXP a, SEXP b)
210 : bates 10 {
211 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
212 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
213 :     *cdims,
214 :     m = adims[0], n = bdims[1], k = adims[1];
215 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
216 : bates 10 double one = 1., zero = 0.;
217 : maechler 534
218 : bates 10 if (bdims[0] != k)
219 : bates 582 error(_("Matrices are not conformable for multiplication"));
220 : bates 10 if (m < 1 || n < 1 || k < 1)
221 : bates 582 error(_("Matrices with zero extents cannot be multiplied"));
222 : bates 338 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
223 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
224 : bates 10 SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
225 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
226 :     cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
227 :     cdims[0] = m; cdims[1] = n;
228 : maechler 951 F77_CALL(dsymm)("R", uplo_P(a), adims, bdims+1, &one,
229 : bates 10 REAL(GET_SLOT(a, Matrix_xSym)), adims,
230 :     REAL(GET_SLOT(b, Matrix_xSym)), bdims,
231 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);
232 :     UNPROTECT(1);
233 :     return val;
234 :     }
235 : maechler 534
236 : bates 631 SEXP dsyMatrix_trf(SEXP x)
237 :     {
238 :     SEXP val = get_factors(x, "BunchKaufman"),
239 :     dimP = GET_SLOT(x, Matrix_DimSym),
240 :     uploP = GET_SLOT(x, Matrix_uploSym);
241 :     int *dims = INTEGER(dimP), *perm, info;
242 :     int lwork = -1, n = dims[0];
243 :     char *uplo = CHAR(STRING_ELT(uploP, 0));
244 :     double tmp, *vx, *work;
245 :    
246 :     if (val != R_NilValue) return val;
247 :     dims = INTEGER(dimP);
248 :     val = PROTECT(NEW_OBJECT(MAKE_CLASS("BunchKaufman")));
249 :     SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
250 :     SET_SLOT(val, Matrix_diagSym, mkString("N"));
251 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
252 :     vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * n));
253 :     AZERO(vx, n * n);
254 :     F77_CALL(dlacpy)(uplo, &n, &n, REAL(GET_SLOT(x, Matrix_xSym)), &n, vx, &n);
255 :     perm = INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, n));
256 :     F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, &tmp, &lwork, &info);
257 :     lwork = (int) tmp;
258 :     work = Calloc(lwork, double);
259 :     F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, work, &lwork, &info);
260 :     if (info) error(_("Lapack routine dsytrf returned error code %d"), info);
261 :     UNPROTECT(1);
262 :     Free(work);
263 :     return set_factors(x, val, "BunchKaufman");
264 :     }
265 :    
266 : bates 643 SEXP dsyMatrix_as_dspMatrix(SEXP from)
267 :     {
268 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dspMatrix"))),
269 :     uplo = GET_SLOT(from, Matrix_uploSym),
270 :     dimP = GET_SLOT(from, Matrix_DimSym);
271 :     int n = *INTEGER(dimP);
272 :    
273 :     SET_SLOT(val, Matrix_rcondSym,
274 :     duplicate(GET_SLOT(from, Matrix_rcondSym)));
275 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
276 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
277 : maechler 952 full_to_packed_double(
278 :     REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
279 :     REAL(GET_SLOT(from, Matrix_xSym)), n,
280 :     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW, NUN);
281 : bates 643 UNPROTECT(1);
282 :     return val;
283 :     }

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