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

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

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