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