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 576 - (view) (download) (as text)

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

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