SCM

SCM Repository

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

Annotation of /pkg/src/dtrMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : maechler 534 /* double (precision) TRiangular Matrices */
2 :    
3 : bates 478 #include "dtrMatrix.h"
4 : bates 10
5 : maechler 633 /* FIXME: Validation works "funny": dtrMatrix_as_dgeMatrix() {below}
6 : maechler 576 * ----- is called *before* the following - presumably in order to
7 :     * apply the higher level validation first
8 :     */
9 : bates 478 SEXP dtrMatrix_validate(SEXP obj)
10 : bates 10 {
11 : bates 592 SEXP val;
12 : maechler 534
13 : bates 592 if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
14 :     "LU", "uplo"))) return val;
15 :     if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),
16 :     "NU", "diag"))) return val;
17 : bates 10 return ScalarLogical(1);
18 :     }
19 :    
20 :     static
21 :     double get_norm(SEXP obj, char *typstr)
22 :     {
23 :     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') {
29 :     work = (double *) R_alloc(dims[0], sizeof(double));
30 :     }
31 :     return F77_CALL(dlantr)(typnm,
32 : bates 296 CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),
33 :     CHAR(asChar(GET_SLOT(obj, Matrix_diagSym))),
34 : bates 10 dims, dims+1,
35 : bates 296 REAL(GET_SLOT(obj, Matrix_xSym)),
36 : bates 10 dims, work);
37 :     }
38 :    
39 :    
40 : bates 478 SEXP dtrMatrix_norm(SEXP obj, SEXP type)
41 : bates 10 {
42 :     return ScalarReal(get_norm(obj, CHAR(asChar(type))));
43 :     }
44 :    
45 :     static
46 :     double set_rcond(SEXP obj, char *typstr)
47 :     {
48 :     char typnm[] = {'\0', '\0'};
49 : bates 342 SEXP rcv = GET_SLOT(obj, Matrix_rcondSym);
50 : bates 10 double rcond = get_double_by_name(rcv, typnm);
51 :    
52 :     typnm[0] = rcond_type(typstr);
53 :     if (R_IsNA(rcond)) {
54 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
55 :     F77_CALL(dtrcon)(typnm,
56 : bates 296 CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),
57 :     CHAR(asChar(GET_SLOT(obj, Matrix_diagSym))),
58 :     dims, REAL(GET_SLOT(obj, Matrix_xSym)),
59 : bates 10 dims, &rcond,
60 :     (double *) R_alloc(3*dims[0], sizeof(double)),
61 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
62 : bates 342 SET_SLOT(obj, Matrix_rcondSym,
63 : bates 10 set_double_by_name(rcv, rcond, typnm));
64 :     }
65 :     return rcond;
66 :     }
67 :    
68 : bates 478 SEXP dtrMatrix_rcond(SEXP obj, SEXP type)
69 : bates 10 {
70 :     return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
71 :     }
72 :    
73 : bates 478 SEXP dtrMatrix_solve(SEXP a)
74 : bates 10 {
75 :     SEXP val = PROTECT(duplicate(a));
76 :     int info, *Dim = INTEGER(GET_SLOT(val, Matrix_DimSym));
77 : bates 296 F77_CALL(dtrtri)(CHAR(asChar(GET_SLOT(val, Matrix_uploSym))),
78 :     CHAR(asChar(GET_SLOT(val, Matrix_diagSym))),
79 :     Dim, REAL(GET_SLOT(val, Matrix_xSym)), Dim, &info);
80 : bates 10 UNPROTECT(1);
81 :     return val;
82 :     }
83 :    
84 : bates 654 SEXP dtrMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
85 : maechler 655 {
86 : bates 654 int cl = asLogical(classed);
87 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
88 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
89 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
90 :     getAttrib(b, R_DimSymbol));
91 : maechler 655 int n = bdims[0], nrhs = bdims[1];
92 : bates 654 int sz = n * nrhs;
93 : bates 296 double one = 1.0;
94 :    
95 : bates 654 if (*adims != *bdims || bdims[1] < 1 || *adims < 1 || *adims != adims[1])
96 :     error(_("Dimensions of system to be solved are inconsistent"));
97 : bates 296 F77_CALL(dtrsm)("L", CHAR(asChar(GET_SLOT(val, Matrix_uploSym))),
98 :     "N", CHAR(asChar(GET_SLOT(val, Matrix_diagSym))),
99 : bates 654 &n, &nrhs, &one, REAL(GET_SLOT(a, Matrix_xSym)), &n,
100 :     Memcpy(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz)),
101 :     REAL(cl ? GET_SLOT(b, Matrix_xSym):b), sz), &n);
102 : bates 296 UNPROTECT(1);
103 :     return val;
104 :     }
105 :    
106 : bates 654 SEXP dtrMatrix_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
107 :     {
108 :     int cl = asLogical(classed), rt = asLogical(right);
109 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
110 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
111 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
112 :     getAttrib(b, R_DimSymbol)),
113 :     *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
114 :     int m, n, sz;
115 :     double one = 1.;
116 :    
117 :     if (!cl && !(isReal(b) && isMatrix(b)))
118 :     error(_("Argument b must be a numeric matrix"));
119 : maechler 688 if (adims[0] != adims[1]) error(_("dtrMatrix in %*% must be square"));
120 : bates 654 m = rt ? bdims[0] : adims[0];
121 :     n = rt ? adims[1] : bdims[1];
122 :     if ((rt && (adims[0] != m)) || (!rt && (bdims[0] != m)))
123 :     error(_("Matrices are not conformable for multiplication"));
124 :     if (m < 1 || n < 1)
125 :     error(_("Matrices with zero extents cannot be multiplied"));
126 :     cdims[0] = m; cdims[1] = n; sz = m * n;
127 :     F77_CALL(dtrmm)(rt ? "R" : "L", CHAR(asChar(GET_SLOT(a, Matrix_uploSym))),
128 :     "N", CHAR(asChar(GET_SLOT(a, Matrix_diagSym))), &m, &n,
129 :     &one, REAL(GET_SLOT(a, Matrix_xSym)), rt ? &n : &m,
130 :     Memcpy(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz)),
131 :     REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), sz),
132 :     rt ? &m : &n);
133 :     UNPROTECT(1);
134 :     return val;
135 :     }
136 :    
137 : bates 478 SEXP dtrMatrix_as_dgeMatrix(SEXP from)
138 : bates 10 {
139 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
140 : maechler 534
141 : maechler 628 SET_SLOT(val, Matrix_rcondSym, duplicate(GET_SLOT(from, Matrix_rcondSym)));
142 : bates 296 SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(from, Matrix_xSym)));
143 : maechler 633 /* Dim < 2 can give a seg.fault problem in make_array_triangular(),
144 :     * by new("dtrMatrix", Dim = 2:2, x=as.double(1:4)) )# length(Dim) !=2 */
145 : maechler 576 if (LENGTH(GET_SLOT(from, Matrix_DimSym)) < 2)
146 : maechler 628 error(_("'Dim' slot has length less than two"));
147 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(from, Matrix_DimSym)));
148 :     SET_SLOT(val, Matrix_DimNamesSym, duplicate(GET_SLOT(from, Matrix_DimNamesSym)));
149 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
150 : bates 296 make_array_triangular(REAL(GET_SLOT(val, Matrix_xSym)), from);
151 : bates 10 UNPROTECT(1);
152 :     return val;
153 :     }
154 :    
155 : bates 478 SEXP dtrMatrix_as_matrix(SEXP from)
156 : bates 10 {
157 :     int *Dim = INTEGER(GET_SLOT(from, Matrix_DimSym));
158 :     int m = Dim[0], n = Dim[1];
159 :     SEXP val = PROTECT(allocMatrix(REALSXP, m, n));
160 :     make_array_triangular(Memcpy(REAL(val),
161 : bates 296 REAL(GET_SLOT(from, Matrix_xSym)), m * n),
162 : bates 10 from);
163 : maechler 628 setAttrib(val, R_DimNamesSymbol, GET_SLOT(from, Matrix_DimNamesSym));
164 : bates 10 UNPROTECT(1);
165 :     return val;
166 :     }
167 :    
168 : bates 478 SEXP dtrMatrix_getDiag(SEXP x)
169 : bates 10 {
170 :     int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
171 :     SEXP ret = PROTECT(allocVector(REALSXP, n)),
172 : bates 296 xv = GET_SLOT(x, Matrix_xSym);
173 : bates 10
174 : maechler 534 if ('U' == CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0]) {
175 : bates 10 for (i = 0; i < n; i++) REAL(ret)[i] = 1.;
176 :     } else {
177 :     for (i = 0; i < n; i++) {
178 :     REAL(ret)[i] = REAL(xv)[i * (n + 1)];
179 :     }
180 :     }
181 :     UNPROTECT(1);
182 :     return ret;
183 :     }
184 : maechler 534
185 : bates 478 SEXP dtrMatrix_dgeMatrix_mm_R(SEXP a, SEXP b)
186 : bates 10 {
187 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
188 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
189 :     m = adims[0], n = bdims[1], k = adims[1];
190 :     SEXP val = PROTECT(duplicate(b));
191 :     double one = 1.;
192 : maechler 534
193 : bates 10 if (bdims[0] != k)
194 : bates 582 error(_("Matrices are not conformable for multiplication"));
195 : bates 10 if (m < 1 || n < 1 || k < 1)
196 : bates 582 error(_("Matrices with zero extents cannot be multiplied"));
197 : bates 10 F77_CALL(dtrmm)("R", CHAR(asChar(GET_SLOT(a, Matrix_uploSym))), "N",
198 :     CHAR(asChar(GET_SLOT(a, Matrix_diagSym))),
199 :     adims, bdims+1, &one,
200 :     REAL(GET_SLOT(a, Matrix_xSym)), adims,
201 :     REAL(GET_SLOT(val, Matrix_xSym)), bdims);
202 :     UNPROTECT(1);
203 :     return val;
204 :     }
205 : bates 597
206 :     SEXP dtrMatrix_as_dtpMatrix(SEXP from)
207 :     {
208 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtpMatrix"))),
209 :     uplo = GET_SLOT(from, Matrix_uploSym),
210 :     diag = GET_SLOT(from, Matrix_diagSym),
211 :     dimP = GET_SLOT(from, Matrix_DimSym);
212 :     int n = *INTEGER(dimP);
213 :    
214 :     SET_SLOT(val, Matrix_rcondSym,
215 :     duplicate(GET_SLOT(from, Matrix_rcondSym)));
216 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
217 :     SET_SLOT(val, Matrix_diagSym, duplicate(diag));
218 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
219 :     full_to_packed(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
220 :     REAL(GET_SLOT(from, Matrix_xSym)), n,
221 :     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW,
222 :     *CHAR(STRING_ELT(diag, 0)) == 'U' ? UNT : NUN);
223 :     UNPROTECT(1);
224 :     return val;
225 :     }

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