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 633 - (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 478 SEXP dtrMatrix_matrix_solve(SEXP a, SEXP b)
85 : bates 296 {
86 :     SEXP val = PROTECT(duplicate(b));
87 :     int *Dim = INTEGER(GET_SLOT(a, Matrix_DimSym)),
88 :     *bDim = INTEGER(getAttrib(val, R_DimSymbol));
89 :     double one = 1.0;
90 :    
91 :     if (bDim[0] != Dim[1])
92 : bates 582 error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
93 : bates 296 Dim[0], Dim[1], bDim[0], bDim[1]);
94 :     F77_CALL(dtrsm)("L", CHAR(asChar(GET_SLOT(val, Matrix_uploSym))),
95 :     "N", CHAR(asChar(GET_SLOT(val, Matrix_diagSym))),
96 :     bDim, bDim+1, &one,
97 :     REAL(GET_SLOT(a, Matrix_xSym)), Dim,
98 :     REAL(val), bDim);
99 :     UNPROTECT(1);
100 :     return val;
101 :     }
102 :    
103 : bates 478 SEXP dtrMatrix_as_dgeMatrix(SEXP from)
104 : bates 10 {
105 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
106 : maechler 534
107 : maechler 628 SET_SLOT(val, Matrix_rcondSym, duplicate(GET_SLOT(from, Matrix_rcondSym)));
108 : bates 296 SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(from, Matrix_xSym)));
109 : maechler 633 /* Dim < 2 can give a seg.fault problem in make_array_triangular(),
110 :     * by new("dtrMatrix", Dim = 2:2, x=as.double(1:4)) )# length(Dim) !=2 */
111 : maechler 576 if (LENGTH(GET_SLOT(from, Matrix_DimSym)) < 2)
112 : maechler 628 error(_("'Dim' slot has length less than two"));
113 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(from, Matrix_DimSym)));
114 :     SET_SLOT(val, Matrix_DimNamesSym, duplicate(GET_SLOT(from, Matrix_DimNamesSym)));
115 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
116 : bates 296 make_array_triangular(REAL(GET_SLOT(val, Matrix_xSym)), from);
117 : bates 10 UNPROTECT(1);
118 :     return val;
119 :     }
120 :    
121 : bates 478 SEXP dtrMatrix_as_matrix(SEXP from)
122 : bates 10 {
123 :     int *Dim = INTEGER(GET_SLOT(from, Matrix_DimSym));
124 :     int m = Dim[0], n = Dim[1];
125 :     SEXP val = PROTECT(allocMatrix(REALSXP, m, n));
126 :     make_array_triangular(Memcpy(REAL(val),
127 : bates 296 REAL(GET_SLOT(from, Matrix_xSym)), m * n),
128 : bates 10 from);
129 : maechler 628 setAttrib(val, R_DimNamesSymbol, GET_SLOT(from, Matrix_DimNamesSym));
130 : bates 10 UNPROTECT(1);
131 :     return val;
132 :     }
133 :    
134 : bates 478 SEXP dtrMatrix_getDiag(SEXP x)
135 : bates 10 {
136 :     int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
137 :     SEXP ret = PROTECT(allocVector(REALSXP, n)),
138 : bates 296 xv = GET_SLOT(x, Matrix_xSym);
139 : bates 10
140 : maechler 534 if ('U' == CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0]) {
141 : bates 10 for (i = 0; i < n; i++) REAL(ret)[i] = 1.;
142 :     } else {
143 :     for (i = 0; i < n; i++) {
144 :     REAL(ret)[i] = REAL(xv)[i * (n + 1)];
145 :     }
146 :     }
147 :     UNPROTECT(1);
148 :     return ret;
149 :     }
150 : maechler 534
151 : bates 478 SEXP dtrMatrix_dgeMatrix_mm(SEXP a, SEXP b)
152 : bates 10 {
153 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
154 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
155 :     m = adims[0], n = bdims[1], k = adims[1];
156 :     SEXP val = PROTECT(duplicate(b));
157 :     double one = 1.;
158 : maechler 534
159 : bates 10 if (bdims[0] != k)
160 : bates 582 error(_("Matrices are not conformable for multiplication"));
161 : bates 10 if (m < 1 || n < 1 || k < 1)
162 : bates 582 error(_("Matrices with zero extents cannot be multiplied"));
163 : bates 10 F77_CALL(dtrmm)("L", CHAR(asChar(GET_SLOT(a, Matrix_uploSym))), "N",
164 :     CHAR(asChar(GET_SLOT(a, Matrix_diagSym))),
165 :     adims, bdims+1, &one,
166 :     REAL(GET_SLOT(a, Matrix_xSym)), adims,
167 :     REAL(GET_SLOT(val, Matrix_xSym)), bdims);
168 :     UNPROTECT(1);
169 :     return val;
170 :     }
171 :    
172 : bates 478 SEXP dtrMatrix_dgeMatrix_mm_R(SEXP a, SEXP b)
173 : bates 10 {
174 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
175 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
176 :     m = adims[0], n = bdims[1], k = adims[1];
177 :     SEXP val = PROTECT(duplicate(b));
178 :     double one = 1.;
179 : maechler 534
180 : bates 10 if (bdims[0] != k)
181 : bates 582 error(_("Matrices are not conformable for multiplication"));
182 : bates 10 if (m < 1 || n < 1 || k < 1)
183 : bates 582 error(_("Matrices with zero extents cannot be multiplied"));
184 : bates 10 F77_CALL(dtrmm)("R", CHAR(asChar(GET_SLOT(a, Matrix_uploSym))), "N",
185 :     CHAR(asChar(GET_SLOT(a, Matrix_diagSym))),
186 :     adims, bdims+1, &one,
187 :     REAL(GET_SLOT(a, Matrix_xSym)), adims,
188 :     REAL(GET_SLOT(val, Matrix_xSym)), bdims);
189 :     UNPROTECT(1);
190 :     return val;
191 :     }
192 : bates 597
193 :     SEXP dtrMatrix_as_dtpMatrix(SEXP from)
194 :     {
195 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtpMatrix"))),
196 :     uplo = GET_SLOT(from, Matrix_uploSym),
197 :     diag = GET_SLOT(from, Matrix_diagSym),
198 :     dimP = GET_SLOT(from, Matrix_DimSym);
199 :     int n = *INTEGER(dimP);
200 :    
201 :     SET_SLOT(val, Matrix_rcondSym,
202 :     duplicate(GET_SLOT(from, Matrix_rcondSym)));
203 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
204 :     SET_SLOT(val, Matrix_diagSym, duplicate(diag));
205 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
206 :     full_to_packed(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
207 :     REAL(GET_SLOT(from, Matrix_xSym)), n,
208 :     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW,
209 :     *CHAR(STRING_ELT(diag, 0)) == 'U' ? UNT : NUN);
210 :     UNPROTECT(1);
211 :     return val;
212 :     }

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