SCM

SCM Repository

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

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

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

revision 1736, Tue Jan 23 17:09:41 2007 UTC revision 1798, Sat Mar 24 14:52:47 2007 UTC
# Line 116  Line 116 
116      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
117          lo = uplo_P(a)[0] == 'L',          lo = uplo_P(a)[0] == 'L',
118          bnz = 10 * A->n;        /* initial estimate of nnz in b */          bnz = 10 * A->n;        /* initial estimate of nnz in b */
119      int *ti = Calloc(bnz, int), i, j, nz, pos = 0;      int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
120      double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);      double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);
121        cs *u = cs_spalloc(A->n, 1,1,1,0);  /* Sparse unit vector */
122        int top;                            /* top of stack */
123        int *xi = Calloc(2*A->n, int);      /* cs_reach uses this workspace */
124    
125      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
126      SET_DimNames(ans, a);      SET_DimNames(ans, a);
127      SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));      SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
128      SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));      SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
129        /* initialize the "constant part" of the sparse unit vector */
130        u->x[0] = 1.;
131        u->p[0] = 0; u->p[1] = 1;
132      bp[0] = 0;      bp[0] = 0;
133      for (j = 0; j < A->n; j++) {      for (j = 0; j < A->n; j++) {
134          AZERO(wrk, A->n);          u->i[0] = j;                    /* u := j'th unit vector */
135          wrk[j] = 1;          /* (wrk[top:n],xi[top:n]) :=  A^{-1} u  : */
136          lo ? cs_lsolve(A, wrk) : cs_usolve(A, wrk);          top = cs_spsolve (A, u, 0, xi, wrk, 0, lo);
137          for (i = 0, nz = 0; i < A->n; i++) if (wrk[i]) nz++;          nz = A->n - top;
138          bp[j + 1] = nz + bp[j];          bp[j + 1] = nz + bp[j];
139          if (bp[j + 1] > bnz) {          if (bp[j + 1] > bnz) {
140              while (bp[j + 1] > bnz) bnz *= 2;              while (bp[j + 1] > bnz) bnz *= 2;
141              ti = Realloc(ti, bnz, int);              ti = Realloc(ti, bnz, int);
142              tx = Realloc(tx, bnz, double);              tx = Realloc(tx, bnz, double);
143          }          }
144          for (i = 0; i < A->n; i++)          if (lo)
145              if (wrk[i]) {ti[pos] = i; tx[pos] = wrk[i]; pos++;}              for(p = top; p < A->n; p++, pos++) {
146                    ti[pos] = xi[p];
147                    tx[pos] = wrk[xi[p]];
148                }
149            else /* upper triagonal */
150                for(p = A->n - 1; p >= top; p--, pos++) {
151                    ti[pos] = xi[p];
152                    tx[pos] = wrk[xi[p]];
153                }
154      }      }
155      nz = bp[A->n];      nz = bp[A->n];
156      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
157      Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);      Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
158    
159      Free(A); Free(ti); Free(tx);      Free(A); Free(ti); Free(tx);
160        Free(wrk); cs_spfree(u); Free(xi);
161      UNPROTECT(1);      UNPROTECT(1);
162      return ans;      return ans;
163  }  }

Legend:
Removed from v.1736  
changed lines
  Added in v.1798

root@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