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 591, Thu Mar 3 05:16:09 2005 UTC revision 597, Thu Mar 3 19:58:59 2005 UTC
# Line 5  Line 5 
5  SEXP dtpMatrix_validate(SEXP obj)  SEXP dtpMatrix_validate(SEXP obj)
6  {  {
7      SEXP val;      SEXP val;
8        int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
9    
10      if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),      if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
11                                             "LU", "uplo"))) return val;                                             "LU", "uplo"))) return val;
12      if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),      if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),
13                                             "NU", "diag"))) return val;                                             "NU", "diag"))) return val;
14        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      return ScalarLogical(1);      return ScalarLogical(1);
18  }  }
19    
# Line 63  Line 67 
67      return ScalarReal(set_rcond(obj, CHAR(asChar(type))));      return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
68  }  }
69    
   
70  SEXP dtpMatrix_solve(SEXP a)  SEXP dtpMatrix_solve(SEXP a)
71  {  {
72      SEXP val = PROTECT(duplicate(a));      SEXP val = PROTECT(duplicate(a));
# Line 75  Line 78 
78      return val;      return val;
79  }  }
80    
81    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  SEXP dtpMatrix_matrix_solve(SEXP a, SEXP b)  SEXP dtpMatrix_matrix_solve(SEXP a, SEXP b)
97  {  {
98      SEXP val = PROTECT(duplicate(b));      SEXP val = PROTECT(duplicate(b));
99      int *Dim = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *Dim = INTEGER(GET_SLOT(a, Matrix_DimSym)),
100          *bDim = INTEGER(getAttrib(val, R_DimSymbol));          *bDim = INTEGER(getAttrib(val, R_DimSymbol));
101        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      int ione = 1, j;      int ione = 1, j;
105    
106      if (bDim[0] != Dim[1])      if (bDim[0] != Dim[1])
107          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"),
108                Dim[0], Dim[1], bDim[0], bDim[1]);                Dim[0], Dim[1], bDim[0], bDim[1]);
109      for (j = 0; j < bDim[1]; j++)      for (j = 0; j < bDim[1]; j++)
110          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)),  
111                          REAL(val) + j * bDim[0], &ione);                          REAL(val) + j * bDim[0], &ione);
112      UNPROTECT(1);      UNPROTECT(1);
113      return val;      return val;
114  }  }
115    
116  SEXP dtpMatrix_as_dtrMatrix(SEXP from)  SEXP dtpMatrix_dgeMatrix_mm(SEXP x, SEXP y)
117  {  {
118      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))),      SEXP val = PROTECT(duplicate(y));
119          uplo = GET_SLOT(from, Matrix_uploSym),      int *xDim = INTEGER(GET_SLOT(x, Matrix_DimSym)),
120          diag = GET_SLOT(from, Matrix_diagSym),          *yDim = INTEGER(GET_SLOT(y, Matrix_DimSym));
121          dimP = GET_SLOT(from, Matrix_DimSym);      int ione = 1, j;
122      int n = *INTEGER(dimP);      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    
127      SET_SLOT(val, Matrix_rcondSym,      if (yDim[0] != xDim[1])
128               duplicate(GET_SLOT(from, Matrix_rcondSym)));          error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
129      SET_SLOT(val, Matrix_DimSym, duplicate(dimP));                xDim[0], xDim[1], yDim[0], yDim[1]);
130      SET_SLOT(val, Matrix_diagSym, duplicate(diag));      for (j = 0; j < yDim[1]; j++)
131      SET_SLOT(val, Matrix_uploSym, duplicate(uplo));          F77_CALL(dtpsv)(uplo, "N", diag, yDim, xx,
132      packed_to_full(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),                          vx + j * yDim[0], &ione);
                    REAL(GET_SLOT(from, Matrix_xSym)), n,  
                    *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW);  
133      UNPROTECT(1);      UNPROTECT(1);
134      return val;      return val;
135  }  }
136    
137  SEXP dtrMatrix_as_dtpMatrix(SEXP from)  SEXP dtpMatrix_as_dtrMatrix(SEXP from)
138  {  {
139      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtpMatrix"))),      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))),
140          uplo = GET_SLOT(from, Matrix_uploSym),          uplo = GET_SLOT(from, Matrix_uploSym),
141          diag = GET_SLOT(from, Matrix_diagSym),          diag = GET_SLOT(from, Matrix_diagSym),
142          dimP = GET_SLOT(from, Matrix_DimSym);          dimP = GET_SLOT(from, Matrix_DimSym);
# Line 127  Line 147 
147      SET_SLOT(val, Matrix_DimSym, duplicate(dimP));      SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
148      SET_SLOT(val, Matrix_diagSym, duplicate(diag));      SET_SLOT(val, Matrix_diagSym, duplicate(diag));
149      SET_SLOT(val, Matrix_uploSym, duplicate(uplo));      SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
150      full_to_packed(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),      packed_to_full(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),
151                     REAL(GET_SLOT(from, Matrix_xSym)), n,                     REAL(GET_SLOT(from, Matrix_xSym)), n,
152                     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW,                     *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW);
                    *CHAR(STRING_ELT(diag, 0)) == 'U' ? UNT : NUN);  
153      UNPROTECT(1);      UNPROTECT(1);
154      return val;      return val;
155  }  }
156    

Legend:
Removed from v.591  
changed lines
  Added in v.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