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 296 - (view) (download) (as text)
Original Path: pkg/src/trMatrix.c

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

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