SCM

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

pkg/src/dtpMatrix.c revision 597, Thu Mar 3 19:58:59 2005 UTC branches/Matrix-mer2/src/dtpMatrix.c revision 985, Fri Oct 21 19:33:37 2005 UTC
# Line 1  Line 1 
1  /* double (precision) Triangular Packed Matrices */  /* double (precision) Triangular Packed Matrices
2     * Note: this means *square* {n x n} matrices
3    */
4    
5  #include "dtpMatrix.h"  #include "dtpMatrix.h"
6    
7  SEXP dtpMatrix_validate(SEXP obj)  SEXP dtpMatrix_validate(SEXP obj)
8  {  {
9      SEXP val;      SEXP val = triangularMatrix_validate(obj);
10        if(isString(val))
11            return(val);
12        else {
13      int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));      int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
   
     if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),  
                                            "LU", "uplo"))) return val;  
     if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),  
                                            "NU", "diag"))) return val;  
     if (dims[0] != dims[1]) return mkString(_("Matrix in not square"));  
14      if (dims[0] != packed_ncol(length(GET_SLOT(obj, Matrix_xSym))))      if (dims[0] != packed_ncol(length(GET_SLOT(obj, Matrix_xSym))))
15          return(mkString(_("Incorrect length of 'x' slot")));          return(mkString(_("Incorrect length of 'x' slot")));
16      return ScalarLogical(1);      return ScalarLogical(1);
17    
18    
19        }
20  }  }
21    
22  static  static
# Line 28  Line 30 
30      if (*typnm == 'I') {      if (*typnm == 'I') {
31          work = (double *) R_alloc(dims[0], sizeof(double));          work = (double *) R_alloc(dims[0], sizeof(double));
32      }      }
33      return F77_CALL(dlantp)(typnm,      return F77_CALL(dlantp)(typnm, uplo_P(obj), diag_P(obj), dims,
34                              CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),                              REAL(GET_SLOT(obj, Matrix_xSym)), work);
                             CHAR(asChar(GET_SLOT(obj, Matrix_diagSym))),  
                             dims, REAL(GET_SLOT(obj, Matrix_xSym)), work);  
35  }  }
36    
37  SEXP dtpMatrix_norm(SEXP obj, SEXP type)  SEXP dtpMatrix_norm(SEXP obj, SEXP type)
# Line 49  Line 49 
49      typnm[0] = rcond_type(typstr);      typnm[0] = rcond_type(typstr);
50      if (R_IsNA(rcond)) {      if (R_IsNA(rcond)) {
51          int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;          int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
52          F77_CALL(dtpcon)(typnm,          F77_CALL(dtpcon)(typnm, uplo_P(obj), diag_P(obj), dims,
53                           CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))),                           REAL(GET_SLOT(obj, Matrix_xSym)), &rcond,
                          CHAR(asChar(GET_SLOT(obj, Matrix_diagSym))),  
                          dims, REAL(GET_SLOT(obj, Matrix_xSym)),  
                          &rcond,  
54                           (double *) R_alloc(3*dims[0], sizeof(double)),                           (double *) R_alloc(3*dims[0], sizeof(double)),
55                           (int *) R_alloc(dims[0], sizeof(int)), &info);                           (int *) R_alloc(dims[0], sizeof(int)), &info);
56          SET_SLOT(obj, Matrix_rcondSym,          SET_SLOT(obj, Matrix_rcondSym,
# Line 71  Line 68 
68  {  {
69      SEXP val = PROTECT(duplicate(a));      SEXP val = PROTECT(duplicate(a));
70      int info, *Dim = INTEGER(GET_SLOT(val, Matrix_DimSym));      int info, *Dim = INTEGER(GET_SLOT(val, Matrix_DimSym));
71      F77_CALL(dtptri)(CHAR(asChar(GET_SLOT(val, Matrix_uploSym))),      F77_CALL(dtptri)(uplo_P(val), diag_P(val), Dim,
72                       CHAR(asChar(GET_SLOT(val, Matrix_diagSym))),                       REAL(GET_SLOT(val, Matrix_xSym)), &info);
                      Dim, REAL(GET_SLOT(val, Matrix_xSym)), &info);  
73      UNPROTECT(1);      UNPROTECT(1);
74      return val;      return val;
75  }  }
# Line 83  Line 79 
79      int n = *INTEGER(GET_SLOT(x, Matrix_DimSym));      int n = *INTEGER(GET_SLOT(x, Matrix_DimSym));
80      SEXP val = PROTECT(allocVector(REALSXP, n));      SEXP val = PROTECT(allocVector(REALSXP, n));
81    
82      if (*CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0)) == 'U') {      if (*diag_P(x) == 'U') {
83          int j;          int j;
84          for (j = 0; j < n; j++) REAL(val)[j] = 1.;          for (j = 0; j < n; j++) REAL(val)[j] = 1.;
85      } else {      } else {
# Line 98  Line 94 
94      SEXP val = PROTECT(duplicate(b));      SEXP val = PROTECT(duplicate(b));
95      int *Dim = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *Dim = INTEGER(GET_SLOT(a, Matrix_DimSym)),
96          *bDim = INTEGER(getAttrib(val, R_DimSymbol));          *bDim = INTEGER(getAttrib(val, R_DimSymbol));
97      char *uplo = CHAR(STRING_ELT(GET_SLOT(a, Matrix_uploSym), 0)),      char *uplo = uplo_P(a), *diag = diag_P(a);
         *diag = CHAR(STRING_ELT(GET_SLOT(a, Matrix_diagSym), 0));  
98      double *ax = REAL(GET_SLOT(a, Matrix_xSym));      double *ax = REAL(GET_SLOT(a, Matrix_xSym));
99      int ione = 1, j;      int ione = 1, j;
100    
# Line 116  Line 111 
111  SEXP dtpMatrix_dgeMatrix_mm(SEXP x, SEXP y)  SEXP dtpMatrix_dgeMatrix_mm(SEXP x, SEXP y)
112  {  {
113      SEXP val = PROTECT(duplicate(y));      SEXP val = PROTECT(duplicate(y));
114        /* Since 'x' is square (n x n ),   dim(x %*% y) = dim(y) */
115      int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
116          *yDim = INTEGER(GET_SLOT(y, Matrix_DimSym));          *yDim = INTEGER(GET_SLOT(y, Matrix_DimSym));
117      int ione = 1, j;      int ione = 1, j;
118      char *uplo = CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0)),      char *uplo = uplo_P(x), *diag = diag_P(x);
         *diag = CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0));  
119      double *xx = REAL(GET_SLOT(x, Matrix_xSym)),      double *xx = REAL(GET_SLOT(x, Matrix_xSym)),
120          *vx = REAL(GET_SLOT(val, Matrix_xSym));          *vx = REAL(GET_SLOT(val, Matrix_xSym));
121    
122      if (yDim[0] != xDim[1])      if (yDim[0] != xDim[1])
123          error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),          error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
124                xDim[0], xDim[1], yDim[0], yDim[1]);                xDim[0], xDim[1], yDim[0], yDim[1]);
125      for (j = 0; j < yDim[1]; j++)      for (j = 0; j < yDim[1]; j++) /* X %*% y[,j]  via BLAS 2 DTPMV(.) */
126          F77_CALL(dtpsv)(uplo, "N", diag, yDim, xx,          F77_CALL(dtpmv)(uplo, "N", diag, yDim, xx,
127                          vx + j * yDim[0], &ione);                          vx + j * yDim[0], &ione);
128      UNPROTECT(1);      UNPROTECT(1);
129      return val;      return val;
130  }  }
131    
132    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        char *uplo = uplo_P(y), *diag = diag_P(y);
140        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                            vx + i, /* incr = */ xDim);
149        UNPROTECT(1);
150        return val;
151    }
152    
153    SEXP dtpMatrix_matrix_mm(SEXP x, SEXP y)
154    {
155        SEXP val = PROTECT(duplicate(y));
156        /* Since 'x' is square (n x n ),   dim(x %*% y) = dim(y) */
157        int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
158            *yDim = INTEGER(getAttrib(y, R_DimSymbol));
159        int ione = 1, j;
160        char *uplo = uplo_P(x), *diag = diag_P(x);
161        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  SEXP dtpMatrix_as_dtrMatrix(SEXP from)  SEXP dtpMatrix_as_dtrMatrix(SEXP from)
174  {  {
175      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))),
176          uplo = GET_SLOT(from, Matrix_uploSym),          uplo = GET_SLOT(from, Matrix_uploSym),
177          diag = GET_SLOT(from, Matrix_diagSym),          diag = GET_SLOT(from, Matrix_diagSym),
178          dimP = GET_SLOT(from, Matrix_DimSym);          dimP = GET_SLOT(from, Matrix_DimSym),
179            dmnP = GET_SLOT(from, Matrix_DimNamesSym);
180      int n = *INTEGER(dimP);      int n = *INTEGER(dimP);
181    
182      SET_SLOT(val, Matrix_rcondSym,      SET_SLOT(val, Matrix_rcondSym,
183               duplicate(GET_SLOT(from, Matrix_rcondSym)));               duplicate(GET_SLOT(from, Matrix_rcondSym)));
184      SET_SLOT(val, Matrix_DimSym, duplicate(dimP));      SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
185        SET_SLOT(val, Matrix_DimNamesSym, duplicate(dmnP));
186      SET_SLOT(val, Matrix_diagSym, duplicate(diag));      SET_SLOT(val, Matrix_diagSym, duplicate(diag));
187      SET_SLOT(val, Matrix_uploSym, duplicate(uplo));      SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
188      packed_to_full(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),      packed_to_full_double(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),
189                     REAL(GET_SLOT(from, Matrix_xSym)), n,                     REAL(GET_SLOT(from, Matrix_xSym)), n,
190                     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW);                     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW);
191      UNPROTECT(1);      UNPROTECT(1);

Legend:
Removed from v.597  
changed lines
  Added in v.985

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