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

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