SCM

SCM Repository

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

Annotation of /pkg/src/dtpMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 591 /* double (precision) Triangular Packed Matrices */
2 :    
3 :     #include "dtpMatrix.h"
4 :    
5 :     SEXP dtpMatrix_validate(SEXP obj)
6 :     {
7 :     SEXP val;
8 : bates 597 int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
9 : bates 591
10 :     if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
11 :     "LU", "uplo"))) return val;
12 :     if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),
13 :     "NU", "diag"))) return val;
14 : bates 597 if (dims[0] != dims[1]) return mkString(_("Matrix in not square"));
15 :     if (dims[0] != packed_ncol(length(GET_SLOT(obj, Matrix_xSym))))
16 :     return(mkString(_("Incorrect length of 'x' slot")));
17 : bates 591 return ScalarLogical(1);
18 :     }
19 :    
20 :     static
21 :     double get_norm(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') {
29 :     work = (double *) R_alloc(dims[0], sizeof(double));
30 :     }
31 :     return F77_CALL(dlantp)(typnm,
32 :     CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),
33 :     CHAR(asChar(GET_SLOT(obj, Matrix_diagSym))),
34 :     dims, REAL(GET_SLOT(obj, Matrix_xSym)), work);
35 :     }
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 :     F77_CALL(dtpcon)(typnm,
53 :     CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),
54 :     CHAR(asChar(GET_SLOT(obj, Matrix_diagSym))),
55 :     dims, REAL(GET_SLOT(obj, Matrix_xSym)),
56 :     &rcond,
57 :     (double *) R_alloc(3*dims[0], sizeof(double)),
58 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
59 :     SET_SLOT(obj, Matrix_rcondSym,
60 :     set_double_by_name(rcv, rcond, typnm));
61 :     }
62 :     return rcond;
63 :     }
64 :    
65 :     SEXP dtpMatrix_rcond(SEXP obj, SEXP type)
66 :     {
67 :     return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
68 :     }
69 :    
70 :     SEXP dtpMatrix_solve(SEXP a)
71 :     {
72 :     SEXP val = PROTECT(duplicate(a));
73 :     int info, *Dim = INTEGER(GET_SLOT(val, Matrix_DimSym));
74 :     F77_CALL(dtptri)(CHAR(asChar(GET_SLOT(val, Matrix_uploSym))),
75 :     CHAR(asChar(GET_SLOT(val, Matrix_diagSym))),
76 :     Dim, REAL(GET_SLOT(val, Matrix_xSym)), &info);
77 :     UNPROTECT(1);
78 :     return val;
79 :     }
80 :    
81 : bates 597 SEXP dtpMatrix_getDiag(SEXP x)
82 :     {
83 :     int n = *INTEGER(GET_SLOT(x, Matrix_DimSym));
84 :     SEXP val = PROTECT(allocVector(REALSXP, n));
85 :    
86 :     if (*CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0)) == 'U') {
87 :     int j;
88 :     for (j = 0; j < n; j++) REAL(val)[j] = 1.;
89 :     } else {
90 :     packed_getDiag(REAL(val), x);
91 :     }
92 :     UNPROTECT(1);
93 :     return val;
94 :     }
95 :    
96 : bates 591 SEXP dtpMatrix_matrix_solve(SEXP a, SEXP b)
97 :     {
98 :     SEXP val = PROTECT(duplicate(b));
99 :     int *Dim = INTEGER(GET_SLOT(a, Matrix_DimSym)),
100 :     *bDim = INTEGER(getAttrib(val, R_DimSymbol));
101 : bates 597 char *uplo = CHAR(STRING_ELT(GET_SLOT(a, Matrix_uploSym), 0)),
102 :     *diag = CHAR(STRING_ELT(GET_SLOT(a, Matrix_diagSym), 0));
103 :     double *ax = REAL(GET_SLOT(a, Matrix_xSym));
104 : bates 591 int ione = 1, j;
105 :    
106 :     if (bDim[0] != Dim[1])
107 :     error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
108 :     Dim[0], Dim[1], bDim[0], bDim[1]);
109 :     for (j = 0; j < bDim[1]; j++)
110 : bates 597 F77_CALL(dtpsv)(uplo, "N", diag, bDim, ax,
111 : bates 591 REAL(val) + j * bDim[0], &ione);
112 :     UNPROTECT(1);
113 :     return val;
114 :     }
115 :    
116 : bates 597 SEXP dtpMatrix_dgeMatrix_mm(SEXP x, SEXP y)
117 : bates 591 {
118 : bates 597 SEXP val = PROTECT(duplicate(y));
119 :     int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
120 :     *yDim = INTEGER(GET_SLOT(y, Matrix_DimSym));
121 :     int ione = 1, j;
122 :     char *uplo = CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0)),
123 :     *diag = CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0));
124 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym)),
125 :     *vx = REAL(GET_SLOT(val, Matrix_xSym));
126 : bates 591
127 : bates 597 if (yDim[0] != xDim[1])
128 :     error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
129 :     xDim[0], xDim[1], yDim[0], yDim[1]);
130 :     for (j = 0; j < yDim[1]; j++)
131 :     F77_CALL(dtpsv)(uplo, "N", diag, yDim, xx,
132 :     vx + j * yDim[0], &ione);
133 : bates 591 UNPROTECT(1);
134 :     return val;
135 :     }
136 :    
137 : bates 597 SEXP dtpMatrix_as_dtrMatrix(SEXP from)
138 : bates 591 {
139 : bates 597 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))),
140 : bates 591 uplo = GET_SLOT(from, Matrix_uploSym),
141 :     diag = GET_SLOT(from, Matrix_diagSym),
142 :     dimP = GET_SLOT(from, Matrix_DimSym);
143 :     int n = *INTEGER(dimP);
144 :    
145 :     SET_SLOT(val, Matrix_rcondSym,
146 :     duplicate(GET_SLOT(from, Matrix_rcondSym)));
147 :     SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
148 :     SET_SLOT(val, Matrix_diagSym, duplicate(diag));
149 :     SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
150 : bates 597 packed_to_full(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),
151 : bates 591 REAL(GET_SLOT(from, Matrix_xSym)), n,
152 : bates 597 *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW);
153 : bates 591 UNPROTECT(1);
154 :     return val;
155 :     }
156 : 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