SCM

SCM Repository

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

Annotation of /pkg/src/syMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 10 #include "syMatrix.h"
2 :    
3 :     SEXP syMatrix_validate(SEXP obj)
4 :     {
5 :     SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
6 :     int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));
7 :     char *val;
8 :    
9 :     if (length(uplo) != 1)
10 :     return ScalarString(mkChar("uplo slot must have length 1"));
11 :     val = CHAR(STRING_ELT(uplo, 0));
12 :     if (strlen(val) != 1)
13 :     return ScalarString(mkChar("uplo[1] must have string length 1"));
14 :     if (toupper(*val) != 'U' && toupper(*val) != 'L')
15 :     return ScalarString(mkChar("uplo[1] must be \"U\" or \"L\""));
16 :     if (Dim[0] != Dim[1])
17 :     return ScalarString(mkChar("Symmetric matrix must be square"));
18 :     return ScalarLogical(1);
19 :     }
20 :    
21 : bates 296 double get_norm_sy(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' || *typnm == 'O') {
29 :     work = (double *) R_alloc(dims[0], sizeof(double));
30 :     }
31 :     return F77_CALL(dlansy)(typnm,
32 :     CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),
33 :     dims, REAL(GET_SLOT(obj, Matrix_xSym)),
34 :     dims, work);
35 :     }
36 :    
37 :     SEXP syMatrix_norm(SEXP obj, SEXP type)
38 :     {
39 :     return ScalarReal(get_norm_sy(obj, CHAR(asChar(type))));
40 :     }
41 :    
42 :     static
43 :     double set_rcond_sy(SEXP obj, char *typstr)
44 :     {
45 :     char typnm[] = {'\0', '\0'};
46 :     SEXP rcv = GET_SLOT(obj, install("rcond"));
47 :     double rcond;
48 :    
49 :     typnm[0] = rcond_type(typstr);
50 :     rcond = get_double_by_name(rcv, typnm);
51 :    
52 :     /* FIXME: Need a factorization here. */
53 :     if (R_IsNA(rcond)) {
54 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
55 :     double anorm = get_norm_sy(obj, "O");
56 :    
57 :     error("Code for set_rcond_sy not yet written");
58 :     F77_CALL(dsycon)(CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),
59 :     dims, REAL(GET_SLOT(obj, Matrix_xSym)),
60 :     dims, INTEGER(GET_SLOT(obj, install("pivot"))),
61 :     &anorm, &rcond,
62 :     (double *) R_alloc(2*dims[0], sizeof(double)),
63 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
64 :     SET_SLOT(obj, install("rcond"),
65 :     set_double_by_name(rcv, rcond, typnm));
66 :     }
67 :     return rcond;
68 :     }
69 :    
70 :     SEXP syMatrix_rcond(SEXP obj, SEXP type)
71 :     {
72 :     /* FIXME: This is a stub */
73 :     /* return ScalarReal(set_rcond_sy(obj, CHAR(asChar(type)))); */
74 :     return ScalarReal(NA_REAL);
75 :     }
76 :    
77 : bates 10 static
78 :     void make_symmetric(double *to, SEXP from, int n)
79 :     {
80 :     int i, j;
81 :     if (toupper(*CHAR(asChar(GET_SLOT(from, Matrix_uploSym)))) == 'U') {
82 :     for (j = 0; j < n; j++) {
83 :     for (i = j+1; i < n; i++) {
84 :     to[i + j*n] = to[j + i*n];
85 :     }
86 :     }
87 :     } else {
88 :     for (j = 1; j < n; j++) {
89 :     for (i = 0; i < j && i < n; i++) {
90 :     to[i + j*n] = to[j + i*n];
91 :     }
92 :     }
93 :     }
94 :     }
95 :    
96 : bates 296 SEXP syMatrix_solve(SEXP a)
97 :     {
98 :     /* FIXME: Write the code */
99 :     error("code for syMatrix_solve not yet written");
100 : bates 311 return R_NilValue;
101 : bates 296 }
102 :    
103 :     SEXP syMatrix_matrix_solve(SEXP a, SEXP b)
104 :     {
105 :     /* FIXME: Write the code */
106 :     error("code for syMatrix_matrix_solve not yet written");
107 : bates 311 return R_NilValue;
108 : bates 296 }
109 :    
110 : bates 10 SEXP syMatrix_as_geMatrix(SEXP from)
111 :     {
112 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix"))),
113 :     rcondSym = install("rcond");
114 :    
115 :     SET_SLOT(val, rcondSym, duplicate(GET_SLOT(from, rcondSym)));
116 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(from, Matrix_xSym)));
117 :     SET_SLOT(val, Matrix_DimSym,
118 :     duplicate(GET_SLOT(from, Matrix_DimSym)));
119 :     make_symmetric(REAL(GET_SLOT(val, Matrix_xSym)), from,
120 :     INTEGER(GET_SLOT(val, Matrix_DimSym))[0]);
121 :     UNPROTECT(1);
122 :     return val;
123 :     }
124 :    
125 :     SEXP syMatrix_as_matrix(SEXP from)
126 :     {
127 :     int n = INTEGER(GET_SLOT(from, Matrix_DimSym))[0];
128 :     SEXP val = PROTECT(allocMatrix(REALSXP, n, n));
129 :    
130 :     make_symmetric(Memcpy(REAL(val),
131 :     REAL(GET_SLOT(from, Matrix_xSym)), n * n),
132 :     from, n);
133 :     UNPROTECT(1);
134 :     return val;
135 :     }
136 :    
137 :     SEXP syMatrix_geMatrix_mm(SEXP a, SEXP b)
138 :     {
139 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
140 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
141 :     *cdims,
142 :     m = adims[0], n = bdims[1], k = adims[1];
143 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));
144 :     double one = 1., zero = 0.;
145 :    
146 :     if (bdims[0] != k)
147 :     error("Matrices are not conformable for multiplication");
148 :     if (m < 1 || n < 1 || k < 1)
149 :     error("Matrices with zero extents cannot be multiplied");
150 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
151 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
152 :     cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
153 :     cdims[0] = m; cdims[1] = n;
154 :     F77_CALL(dsymm)("L", CHAR(asChar(GET_SLOT(a, Matrix_uploSym))),
155 :     adims, bdims+1, &one,
156 :     REAL(GET_SLOT(a, Matrix_xSym)), adims,
157 :     REAL(GET_SLOT(b, Matrix_xSym)), bdims,
158 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);
159 :     UNPROTECT(1);
160 :     return val;
161 :     }
162 :    
163 :     SEXP syMatrix_geMatrix_mm_R(SEXP a, SEXP b)
164 :     {
165 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
166 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
167 :     *cdims,
168 :     m = adims[0], n = bdims[1], k = adims[1];
169 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix")));
170 :     double one = 1., zero = 0.;
171 :    
172 :     if (bdims[0] != k)
173 :     error("Matrices are not conformable for multiplication");
174 :     if (m < 1 || n < 1 || k < 1)
175 :     error("Matrices with zero extents cannot be multiplied");
176 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
177 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
178 :     cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
179 :     cdims[0] = m; cdims[1] = n;
180 :     F77_CALL(dsymm)("R", CHAR(asChar(GET_SLOT(a, Matrix_uploSym))),
181 :     adims, bdims+1, &one,
182 :     REAL(GET_SLOT(a, Matrix_xSym)), adims,
183 :     REAL(GET_SLOT(b, Matrix_xSym)), bdims,
184 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);
185 :     UNPROTECT(1);
186 :     return val;
187 :     }
188 :    

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