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 591, Thu Mar 3 05:16:09 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      if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),          return(val);
12                                             "LU", "uplo"))) return val;      else {
13      if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),          int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
14                                             "NU", "diag"))) return val;          if (dims[0] != packed_ncol(length(GET_SLOT(obj, Matrix_xSym))))
15                return(mkString(_("Incorrect length of 'x' slot")));
16      return ScalarLogical(1);      return ScalarLogical(1);
17    
18    
19        }
20  }  }
21    
22  static  static
# Line 24  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 45  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 63  Line 64 
64      return ScalarReal(set_rcond(obj, CHAR(asChar(type))));      return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
65  }  }
66    
   
67  SEXP dtpMatrix_solve(SEXP a)  SEXP dtpMatrix_solve(SEXP a)
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);
73                       Dim, REAL(GET_SLOT(val, Matrix_xSym)), &info);      UNPROTECT(1);
74        return val;
75    }
76    
77    SEXP dtpMatrix_getDiag(SEXP x)
78    {
79        int n = *INTEGER(GET_SLOT(x, Matrix_DimSym));
80        SEXP val = PROTECT(allocVector(REALSXP, n));
81    
82        if (*diag_P(x) == 'U') {
83            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);      UNPROTECT(1);
89      return val;      return val;
90  }  }
# Line 80  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 = uplo_P(a), *diag = diag_P(a);
98        double *ax = REAL(GET_SLOT(a, Matrix_xSym));
99      int ione = 1, j;      int ione = 1, j;
100    
101      if (bDim[0] != Dim[1])      if (bDim[0] != Dim[1])
102          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"),
103                Dim[0], Dim[1], bDim[0], bDim[1]);                Dim[0], Dim[1], bDim[0], bDim[1]);
104      for (j = 0; j < bDim[1]; j++)      for (j = 0; j < bDim[1]; j++)
105          F77_CALL(dtpsv)(CHAR(asChar(GET_SLOT(val, Matrix_uploSym))),          F77_CALL(dtpsv)(uplo, "N", diag, bDim, ax,
                         "N", CHAR(asChar(GET_SLOT(val, Matrix_diagSym))),  
                         bDim, REAL(GET_SLOT(a, Matrix_xSym)),  
106                          REAL(val) + j * bDim[0], &ione);                          REAL(val) + j * bDim[0], &ione);
107      UNPROTECT(1);      UNPROTECT(1);
108      return val;      return val;
109  }  }
110    
111  SEXP dtpMatrix_as_dtrMatrix(SEXP from)  SEXP dtpMatrix_dgeMatrix_mm(SEXP x, SEXP y)
112  {  {
113      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))),      SEXP val = PROTECT(duplicate(y));
114          uplo = GET_SLOT(from, Matrix_uploSym),      /* Since 'x' is square (n x n ),   dim(x %*% y) = dim(y) */
115          diag = GET_SLOT(from, Matrix_diagSym),      int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
116          dimP = GET_SLOT(from, Matrix_DimSym);          *yDim = INTEGER(GET_SLOT(y, Matrix_DimSym));
117      int n = *INTEGER(dimP);      int ione = 1, j;
118        char *uplo = uplo_P(x), *diag = diag_P(x);
119        double *xx = REAL(GET_SLOT(x, Matrix_xSym)),
120            *vx = REAL(GET_SLOT(val, Matrix_xSym));
121    
122      SET_SLOT(val, Matrix_rcondSym,      if (yDim[0] != xDim[1])
123               duplicate(GET_SLOT(from, Matrix_rcondSym)));          error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
124      SET_SLOT(val, Matrix_DimSym, duplicate(dimP));                xDim[0], xDim[1], yDim[0], yDim[1]);
125      SET_SLOT(val, Matrix_diagSym, duplicate(diag));      for (j = 0; j < yDim[1]; j++) /* X %*% y[,j]  via BLAS 2 DTPMV(.) */
126      SET_SLOT(val, Matrix_uploSym, duplicate(uplo));          F77_CALL(dtpmv)(uplo, "N", diag, yDim, xx,
127      packed_to_full(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),                          vx + j * yDim[0], &ione);
128                     REAL(GET_SLOT(from, Matrix_xSym)), n,      UNPROTECT(1);
129                     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW);      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);      UNPROTECT(1);
150      return val;      return val;
151  }  }
152    
153  SEXP dtrMatrix_as_dtpMatrix(SEXP from)  SEXP dtpMatrix_matrix_mm(SEXP x, SEXP y)
154  {  {
155      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtpMatrix"))),      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)
174    {
175        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      full_to_packed(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),      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);
                    *CHAR(STRING_ELT(diag, 0)) == 'U' ? UNT : NUN);  
191      UNPROTECT(1);      UNPROTECT(1);
192      return val;      return val;
193  }  }
194    

Legend:
Removed from v.591  
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