SCM

SCM Repository

[matrix] Annotation of /pkg/src/trMatrix.c
ViewVC logotype

Annotation of /pkg/src/trMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 70 - (view) (download) (as text)

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 :     CHAR(asChar(GET_SLOT(obj, install("uplo")))),
39 :     CHAR(asChar(GET_SLOT(obj, install("diag")))),
40 :     dims, dims+1,
41 :     REAL(GET_SLOT(obj, install("x"))),
42 :     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 :     CHAR(asChar(GET_SLOT(obj, install("uplo")))),
63 :     CHAR(asChar(GET_SLOT(obj, install("diag")))),
64 :     dims, REAL(GET_SLOT(obj, install("x"))),
65 :     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 :     F77_CALL(dtrtri)(CHAR(asChar(GET_SLOT(val, install("uplo")))),
84 :     CHAR(asChar(GET_SLOT(val, install("diag")))),
85 :     Dim, REAL(GET_SLOT(val, install("x"))), Dim, &info);
86 :     UNPROTECT(1);
87 :     return val;
88 :     }
89 :    
90 :     void make_array_triangular(double *to, SEXP from)
91 :     {
92 :     int i, j, *dims = INTEGER(GET_SLOT(from, Matrix_DimSym));
93 :     int n = dims[0], m = dims[1];
94 :    
95 :     if (toupper(*CHAR(asChar(GET_SLOT(from, install("uplo"))))) == 'U') {
96 :     for (j = 0; j < n; j++) {
97 :     for (i = j+1; i < m; i++) {
98 :     to[i + j*m] = 0.;
99 :     }
100 :     }
101 :     } else {
102 :     for (j = 1; j < n; j++) {
103 :     for (i = 0; i < j && i < m; i++) {
104 :     to[i + j*m] = 0.;
105 :     }
106 :     }
107 :     }
108 :     if (toupper(*CHAR(asChar(GET_SLOT(from, install("diag"))))) == 'U') {
109 :     j = (n < m) ? n : m;
110 :     for (i = 0; i < j; i++) {
111 :     to[i * (m + 1)] = 1.;
112 :     }
113 :     }
114 :     }
115 :    
116 :     SEXP trMatrix_as_geMatrix(SEXP from)
117 :     {
118 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));
119 :    
120 :     SET_SLOT(val, install("rcond"),
121 :     duplicate(GET_SLOT(from, install("rcond"))));
122 :     SET_SLOT(val, install("x"), duplicate(GET_SLOT(from, install("x"))));
123 :     SET_SLOT(val, Matrix_DimSym,
124 :     duplicate(GET_SLOT(from, Matrix_DimSym)));
125 :     make_array_triangular(REAL(GET_SLOT(val, install("x"))), from);
126 :     UNPROTECT(1);
127 :     return val;
128 :     }
129 :    
130 :     SEXP trMatrix_as_matrix(SEXP from)
131 :     {
132 :     int *Dim = INTEGER(GET_SLOT(from, Matrix_DimSym));
133 :     int m = Dim[0], n = Dim[1];
134 :     SEXP val = PROTECT(allocMatrix(REALSXP, m, n));
135 :    
136 :     make_array_triangular(Memcpy(REAL(val),
137 :     REAL(GET_SLOT(from, install("x"))), m * n),
138 :     from);
139 :     UNPROTECT(1);
140 :     return val;
141 :     }
142 :    
143 :     SEXP trMatrix_getDiag(SEXP x)
144 :     {
145 :     int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
146 :     SEXP ret = PROTECT(allocVector(REALSXP, n)),
147 :     xv = GET_SLOT(x, install("x"));
148 :    
149 :     if ('U' == toupper(CHAR(STRING_ELT(GET_SLOT(x, install("diag")), 0))[0])) {
150 :     for (i = 0; i < n; i++) REAL(ret)[i] = 1.;
151 :     } else {
152 :     for (i = 0; i < n; i++) {
153 :     REAL(ret)[i] = REAL(xv)[i * (n + 1)];
154 :     }
155 :     }
156 :     UNPROTECT(1);
157 :     return ret;
158 :     }
159 :    
160 :     SEXP trMatrix_geMatrix_mm(SEXP a, SEXP b)
161 :     {
162 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
163 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
164 :     m = adims[0], n = bdims[1], k = adims[1];
165 :     SEXP val = PROTECT(duplicate(b));
166 :     double one = 1.;
167 :    
168 :     if (bdims[0] != k)
169 :     error("Matrices are not conformable for multiplication");
170 :     if (m < 1 || n < 1 || k < 1)
171 :     error("Matrices with zero extents cannot be multiplied");
172 :     F77_CALL(dtrmm)("L", CHAR(asChar(GET_SLOT(a, Matrix_uploSym))), "N",
173 :     CHAR(asChar(GET_SLOT(a, Matrix_diagSym))),
174 :     adims, bdims+1, &one,
175 :     REAL(GET_SLOT(a, Matrix_xSym)), adims,
176 :     REAL(GET_SLOT(val, Matrix_xSym)), bdims);
177 :     UNPROTECT(1);
178 :     return val;
179 :     }
180 :    
181 :     SEXP trMatrix_geMatrix_mm_R(SEXP a, SEXP b)
182 :     {
183 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
184 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
185 :     m = adims[0], n = bdims[1], k = adims[1];
186 :     SEXP val = PROTECT(duplicate(b));
187 :     double one = 1.;
188 :    
189 :     if (bdims[0] != k)
190 :     error("Matrices are not conformable for multiplication");
191 :     if (m < 1 || n < 1 || k < 1)
192 :     error("Matrices with zero extents cannot be multiplied");
193 :     F77_CALL(dtrmm)("R", CHAR(asChar(GET_SLOT(a, Matrix_uploSym))), "N",
194 :     CHAR(asChar(GET_SLOT(a, Matrix_diagSym))),
195 :     adims, bdims+1, &one,
196 :     REAL(GET_SLOT(a, Matrix_xSym)), adims,
197 :     REAL(GET_SLOT(val, Matrix_xSym)), bdims);
198 :     UNPROTECT(1);
199 :     return val;
200 :     }

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