SCM

SCM Repository

[matrix] Annotation of /branches/Matrix-mer2/src/dtpMatrix.c
ViewVC logotype

Annotation of /branches/Matrix-mer2/src/dtpMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 628 - (view) (download) (as text)
Original Path: pkg/src/dtpMatrix.c

1 : maechler 628 /* double (precision) Triangular Packed Matrices
2 :     * Note: this means *square* {n x n} matrices
3 :     */
4 : bates 591
5 :     #include "dtpMatrix.h"
6 :    
7 :     SEXP dtpMatrix_validate(SEXP obj)
8 :     {
9 :     SEXP val;
10 : bates 597 int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
11 : bates 591
12 :     if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
13 :     "LU", "uplo"))) return val;
14 :     if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),
15 :     "NU", "diag"))) return val;
16 : maechler 628 if (dims[0] != dims[1]) return mkString(_("Matrix is not square"));
17 : bates 597 if (dims[0] != packed_ncol(length(GET_SLOT(obj, Matrix_xSym))))
18 :     return(mkString(_("Incorrect length of 'x' slot")));
19 : bates 591 return ScalarLogical(1);
20 :     }
21 :    
22 :     static
23 :     double get_norm(SEXP obj, char *typstr)
24 :     {
25 :     char typnm[] = {'\0', '\0'};
26 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
27 :     double *work = (double *) NULL;
28 :    
29 :     typnm[0] = norm_type(typstr);
30 :     if (*typnm == 'I') {
31 :     work = (double *) R_alloc(dims[0], sizeof(double));
32 :     }
33 :     return F77_CALL(dlantp)(typnm,
34 :     CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),
35 :     CHAR(asChar(GET_SLOT(obj, Matrix_diagSym))),
36 :     dims, REAL(GET_SLOT(obj, Matrix_xSym)), work);
37 :     }
38 :    
39 :     SEXP dtpMatrix_norm(SEXP obj, SEXP type)
40 :     {
41 :     return ScalarReal(get_norm(obj, CHAR(asChar(type))));
42 :     }
43 :    
44 :     static
45 :     double set_rcond(SEXP obj, char *typstr)
46 :     {
47 :     char typnm[] = {'\0', '\0'};
48 :     SEXP rcv = GET_SLOT(obj, Matrix_rcondSym);
49 :     double rcond = get_double_by_name(rcv, typnm);
50 :    
51 :     typnm[0] = rcond_type(typstr);
52 :     if (R_IsNA(rcond)) {
53 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
54 :     F77_CALL(dtpcon)(typnm,
55 :     CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),
56 :     CHAR(asChar(GET_SLOT(obj, Matrix_diagSym))),
57 :     dims, REAL(GET_SLOT(obj, Matrix_xSym)),
58 :     &rcond,
59 :     (double *) R_alloc(3*dims[0], sizeof(double)),
60 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
61 :     SET_SLOT(obj, Matrix_rcondSym,
62 :     set_double_by_name(rcv, rcond, typnm));
63 :     }
64 :     return rcond;
65 :     }
66 :    
67 :     SEXP dtpMatrix_rcond(SEXP obj, SEXP type)
68 :     {
69 :     return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
70 :     }
71 :    
72 :     SEXP dtpMatrix_solve(SEXP a)
73 :     {
74 :     SEXP val = PROTECT(duplicate(a));
75 :     int info, *Dim = INTEGER(GET_SLOT(val, Matrix_DimSym));
76 :     F77_CALL(dtptri)(CHAR(asChar(GET_SLOT(val, Matrix_uploSym))),
77 :     CHAR(asChar(GET_SLOT(val, Matrix_diagSym))),
78 :     Dim, REAL(GET_SLOT(val, Matrix_xSym)), &info);
79 :     UNPROTECT(1);
80 :     return val;
81 :     }
82 :    
83 : bates 597 SEXP dtpMatrix_getDiag(SEXP x)
84 :     {
85 :     int n = *INTEGER(GET_SLOT(x, Matrix_DimSym));
86 :     SEXP val = PROTECT(allocVector(REALSXP, n));
87 :    
88 :     if (*CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0)) == 'U') {
89 :     int j;
90 :     for (j = 0; j < n; j++) REAL(val)[j] = 1.;
91 :     } else {
92 :     packed_getDiag(REAL(val), x);
93 :     }
94 :     UNPROTECT(1);
95 :     return val;
96 :     }
97 :    
98 : bates 591 SEXP dtpMatrix_matrix_solve(SEXP a, SEXP b)
99 :     {
100 :     SEXP val = PROTECT(duplicate(b));
101 :     int *Dim = INTEGER(GET_SLOT(a, Matrix_DimSym)),
102 :     *bDim = INTEGER(getAttrib(val, R_DimSymbol));
103 : bates 597 char *uplo = CHAR(STRING_ELT(GET_SLOT(a, Matrix_uploSym), 0)),
104 :     *diag = CHAR(STRING_ELT(GET_SLOT(a, Matrix_diagSym), 0));
105 :     double *ax = REAL(GET_SLOT(a, Matrix_xSym));
106 : bates 591 int ione = 1, j;
107 :    
108 :     if (bDim[0] != Dim[1])
109 :     error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
110 :     Dim[0], Dim[1], bDim[0], bDim[1]);
111 :     for (j = 0; j < bDim[1]; j++)
112 : bates 597 F77_CALL(dtpsv)(uplo, "N", diag, bDim, ax,
113 : bates 591 REAL(val) + j * bDim[0], &ione);
114 :     UNPROTECT(1);
115 :     return val;
116 :     }
117 :    
118 : bates 597 SEXP dtpMatrix_dgeMatrix_mm(SEXP x, SEXP y)
119 : bates 591 {
120 : bates 597 SEXP val = PROTECT(duplicate(y));
121 : maechler 628 /* Since 'x' is square (n x n ), dim(x %*% y) = dim(y) */
122 : bates 597 int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
123 :     *yDim = INTEGER(GET_SLOT(y, Matrix_DimSym));
124 :     int ione = 1, j;
125 :     char *uplo = CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0)),
126 :     *diag = CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0));
127 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym)),
128 :     *vx = REAL(GET_SLOT(val, Matrix_xSym));
129 : bates 591
130 : bates 597 if (yDim[0] != xDim[1])
131 :     error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
132 :     xDim[0], xDim[1], yDim[0], yDim[1]);
133 : maechler 628 for (j = 0; j < yDim[1]; j++) /* A %*% x via BLAS 2 DTPMV(.) */
134 : bates 603 F77_CALL(dtpmv)(uplo, "N", diag, yDim, xx,
135 : bates 597 vx + j * yDim[0], &ione);
136 : bates 591 UNPROTECT(1);
137 :     return val;
138 :     }
139 :    
140 : maechler 628 SEXP dgeMatrix_dtpMatrix_mm(SEXP x, SEXP y)
141 :     {
142 :     SEXP val = PROTECT(duplicate(x));
143 :     /* Since 'y' is square (n x n ), dim(x %*% y) = dim(x) */
144 :     int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
145 :     *yDim = INTEGER(GET_SLOT(y, Matrix_DimSym));
146 :     int i;
147 :     char *uplo = CHAR(STRING_ELT(GET_SLOT(y, Matrix_uploSym), 0)),
148 :     *diag = CHAR(STRING_ELT(GET_SLOT(y, Matrix_diagSym), 0));
149 :     double *yx = REAL(GET_SLOT(y, Matrix_xSym)),
150 :     *vx = REAL(GET_SLOT(val, Matrix_xSym));
151 :    
152 :     if (yDim[0] != xDim[1])
153 :     error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
154 :     xDim[0], xDim[1], yDim[0], yDim[1]);
155 :     for (i = 0; i < xDim[0]; i++)/* val[i,] := Y' %*% x[i,] */
156 :     F77_CALL(dtpmv)(uplo, "T", diag, yDim, yx,
157 :     vx + i * xDim[1], /* incr = */ xDim);
158 :     UNPROTECT(1);
159 :     return val;
160 :     }
161 :    
162 : bates 603 SEXP dtpMatrix_matrix_mm(SEXP x, SEXP y)
163 :     {
164 :     SEXP val = PROTECT(duplicate(y));
165 : maechler 628 /* Since 'x' is square (n x n ), dim(x %*% y) = dim(y) */
166 : bates 603 int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
167 :     *yDim = INTEGER(getAttrib(y, R_DimSymbol));
168 :     int ione = 1, j;
169 :     char *uplo = CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0)),
170 :     *diag = CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0));
171 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym));
172 :    
173 :     if (yDim[0] != xDim[1])
174 :     error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
175 :     xDim[0], xDim[1], yDim[0], yDim[1]);
176 :     for (j = 0; j < yDim[1]; j++)
177 :     F77_CALL(dtpmv)(uplo, "N", diag, yDim, xx,
178 :     REAL(val) + j * yDim[0], &ione);
179 :     UNPROTECT(1);
180 :     return val;
181 :     }
182 :    
183 : bates 597 SEXP dtpMatrix_as_dtrMatrix(SEXP from)
184 : bates 591 {
185 : bates 597 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))),
186 : bates 591 uplo = GET_SLOT(from, Matrix_uploSym),
187 :     diag = GET_SLOT(from, Matrix_diagSym),
188 : maechler 628 dimP = GET_SLOT(from, Matrix_DimSym),
189 :     dmnP = GET_SLOT(from, Matrix_DimNamesSym);
190 : bates 591 int n = *INTEGER(dimP);
191 :    
192 :     SET_SLOT(val, Matrix_rcondSym,
193 :     duplicate(GET_SLOT(from, Matrix_rcondSym)));
194 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
195 : maechler 628 SET_SLOT(val, Matrix_DimNamesSym, duplicate(dmnP));
196 : bates 591 SET_SLOT(val, Matrix_diagSym, duplicate(diag));
197 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
198 : bates 597 packed_to_full(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),
199 : bates 591 REAL(GET_SLOT(from, Matrix_xSym)), n,
200 : bates 597 *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW);
201 : bates 591 UNPROTECT(1);
202 :     return val;
203 :     }
204 : bates 597

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