SCM

SCM Repository

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

Diff of /pkg/src/dtCMatrix.c

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

revision 1386, Thu Aug 17 22:30:40 2006 UTC revision 1507, Mon Sep 4 18:41:34 2006 UTC
# Line 2  Line 2 
2  #include "dtCMatrix.h"  #include "dtCMatrix.h"
3  #include "cs_utils.h"  #include "cs_utils.h"
4    
5  SEXP tsc_validate(SEXP x)  SEXP dtCMatrix_validate(SEXP x)
6  {  {
7      return triangularMatrix_validate(x);      return triangularMatrix_validate(x);
8      /* see ./dsCMatrix.c or ./dtpMatrix.c  on how to do more testing here */      /* see ./dsCMatrix.c or ./dtpMatrix.c  on how to do more testing here */
9  }  }
10    
11    #if 0                           /* no longer used */
12  SEXP tsc_to_dgTMatrix(SEXP x)  SEXP tsc_to_dgTMatrix(SEXP x)
13  {  {
14      SEXP ans;      SEXP ans;
# Line 44  Line 44 
44      }      }
45      return ans;      return ans;
46  }  }
47    #endif
48    
49  /**  /**
50   * Derive the column pointer vector for the inverse of L from the parent array   * Derive the column pointer vector for the inverse of L from the parent array
# Line 114  Line 115 
115      return ans;      return ans;
116  }  }
117    
 #if 0  
 SEXP dtCMatrix_solve(SEXP a)  
 {  
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));  
     int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',  
         n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0], nnz,  
         *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),  
         *ap = INTEGER(GET_SLOT(a, Matrix_pSym)), *bi,  
         *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));  
     int bnz = 10 * ap[n];         /* initial estimate of nnz in b */  
     int *ri = Calloc(bnz, int), *ind = Calloc(n, int), j;  
     double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *bx;  
   
     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));  
     SET_SLOT(ans, Matrix_DimNamesSym,  
              duplicate(GET_SLOT(a, Matrix_DimNamesSym)));  
     SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));  
     SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));  
   
     if (!(lo && unit))  
         error("code for non-unit or upper triangular not yet written");  
     /* Initially bp will contain increasing negative values ending at zero. */  
     /* Later we add the negative of bp[0] to all values. */  
     bp[n] = 0;  
     for (j = n - 1; j >= 0; j--) { /* columns in reverse order */  
         int i, i1 = ap[j], i2 = ap[j + 1], k, nr;  
         if (i1 < i2) AZERO(ind, n);  
         for (i = i1; i < i2; i++) {  
             ind[ai[i]] = 1;  
             for (k = -bp[ai[i] + 1]; k < -bp[ai[i]]; k++) ind[ri[k]] = 1;  
         }  
         for (k = 0, nr = 0; k < n; k++) if (ind[k]) nr++;  
         if ((nr - bp[j + 1]) > bnz) {  
             while (nr > (bnz + bp[j + 1])) bnz *= 2;  
             ri = Realloc(ri, bnz, int);  
         }  
         bp[j] = bp[j + 1] - nr;  
         for (k = 0, i = -bp[j + 1]; k < n; k++) if (ind[k]) ri[i++] = k;  
     }  
     bnz = -bp[0];  
     bi = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, bnz));  
     bx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, bnz));  
     for (j = 0; j < n; j++) {  
         int bpnew = bp[j] + bnz;  
         Memcpy(bi + bpnew, ri - bp[j], bp[j + 1] - bp[j]);  
         bp[j] = bpnew;  
     }  
     /* insert code for calculating the actual values here */  
     for (j = 0; j < bnz; j++) bx[j] = 1;  
   
     Free(ind); Free(ri);  
     UNPROTECT(1);  
     return ans;  
 }  
 #else  
118  SEXP dtCMatrix_solve(SEXP a)  SEXP dtCMatrix_solve(SEXP a)
119  {  {
120      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
# Line 207  Line 153 
153      UNPROTECT(1);      UNPROTECT(1);
154      return ans;      return ans;
155  }  }
 #endif  
156    
157  SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)  SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
158  {  {

Legend:
Removed from v.1386  
changed lines
  Added in v.1507

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