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

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