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

revision 603, Fri Mar 4 12:53:01 2005 UTC revision 628, Thu Mar 10 17:57:47 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    
# Line 11  Line 13 
13                                             "LU", "uplo"))) return val;                                             "LU", "uplo"))) return val;
14      if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),      if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),
15                                             "NU", "diag"))) return val;                                             "NU", "diag"))) return val;
16      if (dims[0] != dims[1]) return mkString(_("Matrix in not square"));      if (dims[0] != dims[1]) return mkString(_("Matrix is not square"));
17      if (dims[0] != packed_ncol(length(GET_SLOT(obj, Matrix_xSym))))      if (dims[0] != packed_ncol(length(GET_SLOT(obj, Matrix_xSym))))
18          return(mkString(_("Incorrect length of 'x' slot")));          return(mkString(_("Incorrect length of 'x' slot")));
19      return ScalarLogical(1);      return ScalarLogical(1);
# Line 116  Line 118 
118  SEXP dtpMatrix_dgeMatrix_mm(SEXP x, SEXP y)  SEXP dtpMatrix_dgeMatrix_mm(SEXP x, SEXP y)
119  {  {
120      SEXP val = PROTECT(duplicate(y));      SEXP val = PROTECT(duplicate(y));
121        /* Since 'x' is square (n x n ),   dim(x %*% y) = dim(y) */
122      int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
123          *yDim = INTEGER(GET_SLOT(y, Matrix_DimSym));          *yDim = INTEGER(GET_SLOT(y, Matrix_DimSym));
124      int ione = 1, j;      int ione = 1, j;
# Line 127  Line 130 
130      if (yDim[0] != xDim[1])      if (yDim[0] != xDim[1])
131          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"),
132                xDim[0], xDim[1], yDim[0], yDim[1]);                xDim[0], xDim[1], yDim[0], yDim[1]);
133      for (j = 0; j < yDim[1]; j++)      for (j = 0; j < yDim[1]; j++) /* A %*% x  via BLAS 2 DTPMV(.) */
134          F77_CALL(dtpmv)(uplo, "N", diag, yDim, xx,          F77_CALL(dtpmv)(uplo, "N", diag, yDim, xx,
135                          vx + j * yDim[0], &ione);                          vx + j * yDim[0], &ione);
136      UNPROTECT(1);      UNPROTECT(1);
137      return val;      return val;
138  }  }
139    
140    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  SEXP dtpMatrix_matrix_mm(SEXP x, SEXP y)  SEXP dtpMatrix_matrix_mm(SEXP x, SEXP y)
163  {  {
164      SEXP val = PROTECT(duplicate(y));      SEXP val = PROTECT(duplicate(y));
165        /* Since 'x' is square (n x n ),   dim(x %*% y) = dim(y) */
166      int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
167          *yDim = INTEGER(getAttrib(y, R_DimSymbol));          *yDim = INTEGER(getAttrib(y, R_DimSymbol));
168      int ione = 1, j;      int ione = 1, j;
# Line 159  Line 185 
185      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))),
186          uplo = GET_SLOT(from, Matrix_uploSym),          uplo = GET_SLOT(from, Matrix_uploSym),
187          diag = GET_SLOT(from, Matrix_diagSym),          diag = GET_SLOT(from, Matrix_diagSym),
188          dimP = GET_SLOT(from, Matrix_DimSym);          dimP = GET_SLOT(from, Matrix_DimSym),
189            dmnP = GET_SLOT(from, Matrix_DimNamesSym);
190      int n = *INTEGER(dimP);      int n = *INTEGER(dimP);
191    
192      SET_SLOT(val, Matrix_rcondSym,      SET_SLOT(val, Matrix_rcondSym,
193               duplicate(GET_SLOT(from, Matrix_rcondSym)));               duplicate(GET_SLOT(from, Matrix_rcondSym)));
194      SET_SLOT(val, Matrix_DimSym, duplicate(dimP));      SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
195        SET_SLOT(val, Matrix_DimNamesSym, duplicate(dmnP));
196      SET_SLOT(val, Matrix_diagSym, duplicate(diag));      SET_SLOT(val, Matrix_diagSym, duplicate(diag));
197      SET_SLOT(val, Matrix_uploSym, duplicate(uplo));      SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
198      packed_to_full(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),      packed_to_full(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),

Legend:
Removed from v.603  
changed lines
  Added in v.628

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