SCM

SCM Repository

[matrix] Diff of /pkg/Matrix/src/dense.c
ViewVC logotype

Diff of /pkg/Matrix/src/dense.c

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

revision 2114, Mon Feb 18 17:15:41 2008 UTC revision 2115, Sat Feb 23 09:23:17 2008 UTC
# Line 262  Line 262 
262         argument for symmetry.         argument for symmetry.
263      */      */
264      CHM_SP chxs = cholmod_dense_to_sparse(chxd, 1, &c);      CHM_SP chxs = cholmod_dense_to_sparse(chxd, 1, &c);
265      int Rkind = (chxd->xtype == CHOLMOD_REAL) ? Real_KIND(x) : 0;      int Rkind = (chxd->xtype == CHOLMOD_REAL) ? Real_KIND2(x) : 0;
266        /* Note: when 'x' was integer Matrix, Real_KIND(x) = -1, but *_KIND2(.) = 0 */
267    
268      UNPROTECT(1);      UNPROTECT(1);
269      /* chm_sparse_to_SEXP() *could* deal with symmetric      /* chm_sparse_to_SEXP() *could* deal with symmetric
# Line 274  Line 275 
275  }  }
276    
277    
278  /* FIXME: generalize this to  dense_band() : */  SEXP dense_band(SEXP x, SEXP k1P, SEXP k2P)
   
 SEXP ddense_band(SEXP x, SEXP k1P, SEXP k2P)  
279  /* Always returns a full matrix with entries outside the band zeroed  /* Always returns a full matrix with entries outside the band zeroed
280   * Class of the value can be dtrMatrix or dgeMatrix   * Class of the value can be [dln]trMatrix or [dln]geMatrix
281   */   */
282  {  {
283      int k1 = asInteger(k1P), k2 = asInteger(k2P);      int k1 = asInteger(k1P), k2 = asInteger(k2P);
284    
285      if (k1 > k2) {      if (k1 > k2) {
286          error(_("Lower band %d > upper band %d"), k1, k2);          error(_("Lower band %d > upper band %d"), k1, k2);
287            return R_NilValue; /* -Wall */
288      }      }
289      else {      else {
290          SEXP aa, ans = PROTECT(dup_mMatrix_as_dgeMatrix(x));          SEXP ans = PROTECT(dup_mMatrix_as_geMatrix(x));
291          int *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),          int *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
292              i, j, m = adims[0], n = adims[1],              j, m = adims[0], n = adims[1],
293              sqr = (adims[0] == adims[1]),              sqr = (adims[0] == adims[1]),
294              tru = (k1 >= 0), trl = (k2 <= 0);              tru = (k1 >= 0), trl = (k2 <= 0);
295          double *ax = REAL(GET_SLOT(ans, Matrix_xSym));          const char *cl = class_P(ans);
296            enum dense_enum { ddense, ldense, ndense
297            } M_type = ( (cl[0] == 'd') ? ddense :
298                        ((cl[0] == 'l') ? ldense : ndense));
299    
300          for (j = 0; j < n; j++) {  
301              int i1 = j - k2, i2 = j + 1 - k1;  #define SET_ZERO_OUTSIDE                                \
302              for (i = 0; i < i1; i++) ax[i + j * m] = 0.;          for (j = 0; j < n; j++) {                       \
303              for (i = i2; i < m; i++) ax[i + j * m] = 0.;              int i, i1 = j - k2, i2 = j + 1 - k1;        \
304                for (i = 0; i < i1; i++) xx[i + j * m] = 0; \
305                for (i = i2; i < m; i++) xx[i + j * m] = 0; \
306          }          }
307          if (!sqr || (!tru && !trl)) { /* return the dgeMatrix */  
308            if(M_type == ddense) {
309                double *xx = REAL(GET_SLOT(ans, Matrix_xSym));
310                SET_ZERO_OUTSIDE
311            }
312            else { /* (M_type == ldense || M_type == ndense) */
313                int *xx = LOGICAL(GET_SLOT(ans, Matrix_xSym));
314                SET_ZERO_OUTSIDE
315            }
316    
317            if (!sqr || (!tru && !trl)) { /* return the *geMatrix */
318              UNPROTECT(1);              UNPROTECT(1);
319              return ans;              return ans;
320          }          }
321          /* Copy ans to a dtrMatrix object (must be square) */          else {
322                /* Copy ans to a *trMatrix object (must be square) */
323                SEXP aa= PROTECT(NEW_OBJECT(MAKE_CLASS(M_type == ddense? "dtrMatrix":
324                                                       (M_type== ldense? "ltrMatrix"
325                                                        : "ntrMatrix"))));
326          /* Because slots of ans are freshly allocated and ans will not be          /* Because slots of ans are freshly allocated and ans will not be
327           * used, we use the slots themselves and don't duplicate */           * used, we use the slots themselves and don't duplicate */
         aa = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix")));  
328          SET_SLOT(aa, Matrix_xSym, GET_SLOT(ans, Matrix_xSym));          SET_SLOT(aa, Matrix_xSym, GET_SLOT(ans, Matrix_xSym));
329          SET_SLOT(aa, Matrix_DimSym, GET_SLOT(ans, Matrix_DimSym));          SET_SLOT(aa, Matrix_DimSym, GET_SLOT(ans, Matrix_DimSym));
330          SET_SLOT(aa, Matrix_DimNamesSym, GET_SLOT(ans, Matrix_DimNamesSym));          SET_SLOT(aa, Matrix_DimNamesSym, GET_SLOT(ans, Matrix_DimNamesSym));
# Line 315  Line 333 
333          UNPROTECT(2);          UNPROTECT(2);
334          return aa;          return aa;
335      }      }
336        }
337      return R_NilValue;          /* -Wall */      return R_NilValue;          /* -Wall */
338  }  }
339    
340  SEXP ddense_to_symmetric(SEXP x, SEXP uplo)  SEXP dense_to_symmetric(SEXP x, SEXP uplo, SEXP symm_test)
341  /* Class of the value will be dsyMatrix  /* Class of result will be [dln]syMatrix */
  */  
342  {  {
343      SEXP xx = PROTECT(dup_mMatrix_as_dgeMatrix(x)),      int symm_tst = asLogical(symm_test);
344          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))), dns;      SEXP dx = PROTECT(dup_mMatrix_as_geMatrix(x));
345        SEXP ans, dns;
346        const char *cl = class_P(dx);
347        /* same as in ..._geMatrix() above:*/
348        enum dense_enum { ddense, ldense, ndense
349        } M_type = ( (cl[0] == 'd') ? ddense :
350                    ((cl[0] == 'l') ? ldense : ndense));
351    
352      /* Copy xx to ans;      if(symm_tst) {
353       * Because slots of xx are freshly allocated and xx will not be          int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0], i,j;
354       * used, we use the slots themselves and don't duplicate */          if(n != adims[1]) {
355      dns = GET_SLOT(xx, Matrix_DimNamesSym);              UNPROTECT(1);
356                error(_("ddense_to_symmetric(): matrix is not square!"));
357                return R_NilValue; /* -Wall */
358            }
359    #define CHECK_SYMMETRIC                                                 \
360            for (j = 0; j < n; j++)                                         \
361                for (i = 0; i < j; i++)                                     \
362                    if(xx[j * n + i] != xx[i * n + j]) {                    \
363                        UNPROTECT(1);                                       \
364                        error(_("matrix is not symmetric [%d,%d]"), i+1, j+1); \
365                        return R_NilValue; /* -Wall */                      \
366                    }
367            if(M_type == ddense) {
368    
369                double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
370                CHECK_SYMMETRIC
371            }
372            else { /* (M_type == ldense || M_type == ndense) */
373    
374                int *xx = LOGICAL(GET_SLOT(dx, Matrix_xSym));
375                CHECK_SYMMETRIC
376            }
377        }
378    
379        dns = GET_SLOT(dx, Matrix_DimNamesSym);
380      if(!equal_string_vectors(VECTOR_ELT(dns,0),      if(!equal_string_vectors(VECTOR_ELT(dns,0),
381                               VECTOR_ELT(dns,1))) {                               VECTOR_ELT(dns,1))) {
382          /* need _symmetric_ dimnames */          /* need _symmetric_ dimnames */
# Line 338  Line 386 
386              SET_VECTOR_ELT(dns,1, VECTOR_ELT(dns,0));              SET_VECTOR_ELT(dns,1, VECTOR_ELT(dns,0));
387      }      }
388    
389      SET_SLOT(ans, Matrix_xSym,        GET_SLOT(xx, Matrix_xSym));      /* Copy dx to ans;
390      SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(xx, Matrix_DimSym));       * Because slots of dx are freshly allocated and dx will not be
391         * used, we use the slots themselves and don't duplicate */
392        ans = PROTECT(NEW_OBJECT(MAKE_CLASS( M_type == ddense ? "dsyMatrix" :
393                                            (M_type == ldense ? "lsyMatrix" :
394                                             "nsyMatrix"))));
395        SET_SLOT(ans, Matrix_xSym,        GET_SLOT(dx, Matrix_xSym));
396        SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(dx, Matrix_DimSym));
397      SET_SLOT(ans, Matrix_DimNamesSym, dns);      SET_SLOT(ans, Matrix_DimNamesSym, dns);
398      SET_SLOT(ans, Matrix_uploSym,     ScalarString(asChar(uplo)));      SET_SLOT(ans, Matrix_uploSym,     ScalarString(asChar(uplo)));
399    
# Line 350  Line 404 
404  SEXP ddense_symmpart(SEXP x)  SEXP ddense_symmpart(SEXP x)
405  /* Class of the value will be dsyMatrix */  /* Class of the value will be dsyMatrix */
406  {  {
407      SEXP xx = PROTECT(dup_mMatrix_as_dgeMatrix(x));      SEXP dx = PROTECT(dup_mMatrix_as_dgeMatrix(x));
408      int *adims = INTEGER(GET_SLOT(xx, Matrix_DimSym)), n = adims[0];      int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
409    
410      if(n != adims[1]) {      if(n != adims[1]) {
411          UNPROTECT(1);          UNPROTECT(1);
412          error(_("matrix is not square! (symmetric part)"));          error(_("matrix is not square! (symmetric part)"));
413            return R_NilValue; /* -Wall */
414      } else {      } else {
415          SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))), dns;          SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))), dns;
416          double *ax = REAL(GET_SLOT(xx, Matrix_xSym));          double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
417          int i,j;          int i,j;
418    
419          /* only need to assign the *upper* triangle (uplo = "U");          /* only need to assign the *upper* triangle (uplo = "U");
420           * noting that diagonal remains unchanged */           * noting that diagonal remains unchanged */
421          for (j = 0; j < n; j++) {          for (j = 0; j < n; j++) {
422              for (i = 0; i < j; i++) {              for (i = 0; i < j; i++) {
423                  ax[j * n + i] = (ax[j * n + i] + ax[i * n + j]) / 2.;                  xx[j * n + i] = (xx[j * n + i] + xx[i * n + j]) / 2.;
424              }              }
425          }          }
426    
427          /* Copy xx to ans;          dns = GET_SLOT(dx, Matrix_DimNamesSym);
          * Because slots of xx are freshly allocated and xx will not be  
          * used, we use the slots themselves and don't duplicate */  
   
         dns = GET_SLOT(xx, Matrix_DimNamesSym);  
428          if(!equal_string_vectors(VECTOR_ELT(dns,0),          if(!equal_string_vectors(VECTOR_ELT(dns,0),
429                                   VECTOR_ELT(dns,1))) {                                   VECTOR_ELT(dns,1))) {
430              /* need _symmetric_ dimnames */              /* need _symmetric_ dimnames */
431              SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));/* ==> uplo = "U" */              SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));/* ==> uplo = "U" */
432          }          }
433    
434          SET_SLOT(ans, Matrix_xSym,        GET_SLOT(xx, Matrix_xSym));          /* Copy dx to ans;
435          SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(xx, Matrix_DimSym));           * Because slots of dx are freshly allocated and dx will not be
436             * used, we use the slots themselves and don't duplicate */
437    
438            SET_SLOT(ans, Matrix_xSym,        GET_SLOT(dx, Matrix_xSym));
439            SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(dx, Matrix_DimSym));
440          SET_SLOT(ans, Matrix_DimNamesSym, dns);          SET_SLOT(ans, Matrix_DimNamesSym, dns);
441          SET_SLOT(ans, Matrix_uploSym,     mkString("U"));          SET_SLOT(ans, Matrix_uploSym,     mkString("U"));
442    
# Line 394  Line 449 
449  SEXP ddense_skewpart(SEXP x)  SEXP ddense_skewpart(SEXP x)
450  /* Class of the value will be dgeMatrix */  /* Class of the value will be dgeMatrix */
451  {  {
452      SEXP xx = PROTECT(dup_mMatrix_as_dgeMatrix(x));      SEXP dx = PROTECT(dup_mMatrix_as_dgeMatrix(x));
453      int *adims = INTEGER(GET_SLOT(xx, Matrix_DimSym)), n = adims[0];      int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
454    
455      if(n != adims[1]) {      if(n != adims[1]) {
456          UNPROTECT(1);          UNPROTECT(1);
457          error(_("matrix is not square! (skew-symmetric part)"));          error(_("matrix is not square! (skew-symmetric part)"));
458            return R_NilValue; /* -Wall */
459      } else {      } else {
460          SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), dns;          SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), dns;
461          double *ax = REAL(GET_SLOT(xx, Matrix_xSym));          double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
462          int i,j;          int i,j;
463    
464          for (j = 0; j < n; j++) {          for (j = 0; j < n; j++) {
465              ax[j * n + j] = 0.;              xx[j * n + j] = 0.;
466              for (i = 0; i < j; i++) {              for (i = 0; i < j; i++) {
467                  double s = (ax[j * n + i] - ax[i * n + j]) / 2.;                  double s = (xx[j * n + i] - xx[i * n + j]) / 2.;
468                  ax[j * n + i] =  s;                  xx[j * n + i] =  s;
469                  ax[i * n + j] = -s;                  xx[i * n + j] = -s;
470              }              }
471          }          }
472    
473          /* Copy xx to ans;          dns = GET_SLOT(dx, Matrix_DimNamesSym);
          * Because slots of xx are freshly allocated and xx will not be  
          * used, we use the slots themselves and don't duplicate */  
   
         dns = GET_SLOT(xx, Matrix_DimNamesSym);  
474          if(!equal_string_vectors(VECTOR_ELT(dns,0),          if(!equal_string_vectors(VECTOR_ELT(dns,0),
475                                   VECTOR_ELT(dns,1))) {                                   VECTOR_ELT(dns,1))) {
476              /* need _symmetric_ dimnames */              /* need _symmetric_ dimnames */
477              SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));/* uplo = "U" */              SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));/* uplo = "U" */
478          }          }
479    
480          SET_SLOT(ans, Matrix_xSym,        GET_SLOT(xx, Matrix_xSym));          /* Copy dx to ans;
481          SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(xx, Matrix_DimSym));           * Because slots of dx are freshly allocated and dx will not be
482             * used, we use the slots themselves and don't duplicate */
483    
484            SET_SLOT(ans, Matrix_xSym,        GET_SLOT(dx, Matrix_xSym));
485            SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(dx, Matrix_DimSym));
486          SET_SLOT(ans, Matrix_DimNamesSym, dns);          SET_SLOT(ans, Matrix_DimNamesSym, dns);
487          SET_SLOT(ans, Matrix_uploSym,     mkString("U"));          SET_SLOT(ans, Matrix_uploSym,     mkString("U"));
488    

Legend:
Removed from v.2114  
changed lines
  Added in v.2115

R-Forge@R-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge